diff --git a/include/aspect/utilities.h b/include/aspect/utilities.h index 779b018ca25..c5862ac01af 100644 --- a/include/aspect/utilities.h +++ b/include/aspect/utilities.h @@ -1060,6 +1060,13 @@ namespace aspect * entry in the list must match one of the allowed operations. */ std::vector create_model_operator_list(const std::vector &operator_names); + + /** + * Create matrix with unit at independent indices + */ + template + SymmetricTensor<2,dim> symmetric_independent_component_matrix (const unsigned int k); + } } diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc index 975bd8c8909..61f68ccb4a4 100644 --- a/source/material_model/visco_plastic.cc +++ b/source/material_model/visco_plastic.cc @@ -22,6 +22,7 @@ #include #include #include +#include namespace aspect { @@ -354,6 +355,16 @@ namespace aspect evaluate(const MaterialModel::MaterialModelInputs &in, MaterialModel::MaterialModelOutputs &out) const { + // set up additional output for the derivatives + MaterialModel::MaterialModelDerivatives *derivatives; + derivatives = out.template get_additional_output >(); + + // If we use the Newton solver, get the derivative scaling factor, + // otherwise set it to zero. + double derivative_scaling_factor = 0.0; + if (this->get_parameters().nonlinear_solver == Parameters::NonlinearSolver::Newton_Stokes) + derivative_scaling_factor = this->get_newton_handler().get_newton_derivative_scaling_factor(); + // Loop through points for (unsigned int i=0; i < in.temperature.size(); ++i) { @@ -361,6 +372,7 @@ namespace aspect const double pressure = in.pressure[i]; const std::vector &composition = in.composition[i]; const std::vector volume_fractions = compute_volume_fractions(composition); + const SymmetricTensor<2,dim> strain_rate = in.strain_rate[i]; // Averaging composition-field dependent properties @@ -397,7 +409,10 @@ namespace aspect // TODO: This is only consistent with viscosity averaging if the arithmetic averaging // scheme is chosen. It would be useful to have a function to calculate isostress viscosities. const std::vector composition_viscosities = - calculate_isostrain_viscosities(volume_fractions, pressure, temperature, composition, in.strain_rate[i],viscous_flow_law,yield_mechanism); + calculate_isostrain_viscosities(volume_fractions, pressure, temperature, composition, strain_rate,viscous_flow_law,yield_mechanism); + + std::vector > composition_viscosities_derivatives(volume_fractions.size()); + std::vector composition_dviscosities_dpressure(volume_fractions.size()); // The isostrain condition implies that the viscosity averaging should be arithmetic (see above). // We have given the user freedom to apply alternative bounds, because in diffusion-dominated @@ -405,6 +420,71 @@ namespace aspect // of compositional field viscosities is consistent with any averaging scheme. out.viscosities[i] = average_value(composition, composition_viscosities, viscosity_averaging); + // compute derivatives if nessesary + if (derivative_scaling_factor != 0 && derivatives != NULL) + { + const double finite_difference_accuracy = 1e-7; + + // For each independent component, compute the derivative. + for (unsigned int independent_component_k = 0; independent_component_k < SymmetricTensor<2,dim>::n_independent_components; independent_component_k++) + { + const TableIndices<2> indices_ij = SymmetricTensor<2,dim>::unrolled_to_component_indices (independent_component_k); + SymmetricTensor<2,dim> strain_rate_component = strain_rate + std::max(std::fabs(strain_rate[indices_ij]), min_strain_rate) * finite_difference_accuracy + * Utilities::symmetric_independent_component_matrix(independent_component_k); + std::vector eta_component = calculate_isostrain_viscosities(volume_fractions, pressure, temperature, composition, strain_rate_component,viscous_flow_law,yield_mechanism); + + // For each composition of the independent component, compute the derivative. + for (unsigned int composition_i = 0; composition_i < eta_component.size(); composition_i++) + { + // compute the difference between the viscosity with and without the strain-rate difference. + double derivative_component = eta_component[composition_i] - composition_viscosities[composition_i]; + if (derivative_component != 0) + { + // when the difference is non-zero, divide by the difference. + derivative_component /= std::max(std::fabs(strain_rate_component[indices_ij]), min_strain_rate) * finite_difference_accuracy; + } + composition_viscosities_derivatives[composition_i][indices_ij] = derivative_component; + } + } + + /** + * Now compute the derivative of the viscosity to the pressure + */ + double pressure_difference = in.pressure[i] + (std::fabs(in.pressure[i]) * finite_difference_accuracy); + + std::vector pressure_difference_eta = calculate_isostrain_viscosities(volume_fractions, pressure_difference, temperature, composition, strain_rate,viscous_flow_law,yield_mechanism); + + + for (unsigned int composition_i = 0; composition_i < pressure_difference_eta.size(); composition_i++) + { + double deriv_pressure = pressure_difference_eta[composition_i] - composition_viscosities[composition_i]; + if (pressure_difference_eta[composition_i] != 0) + { + if (in.pressure[i] != 0) + { + deriv_pressure /= std::fabs(in.pressure[i]) * finite_difference_accuracy; + } + else + { + deriv_pressure = 0; + } + } + composition_dviscosities_dpressure[composition_i] = deriv_pressure; + } + + double viscosity_averaging_p = 0; // Geometric + if (viscosity_averaging == harmonic) + viscosity_averaging_p = -1; + if (viscosity_averaging == arithmetic) + viscosity_averaging_p = 1; + if (viscosity_averaging == maximum_composition) + viscosity_averaging_p = 1000; + + + derivatives->viscosity_derivative_wrt_strain_rate[i] = Utilities::derivative_of_weighted_p_norm_average(out.viscosities[i],volume_fractions, composition_viscosities, composition_viscosities_derivatives, viscosity_averaging_p); + derivatives->viscosity_derivative_wrt_pressure[i] = Utilities::derivative_of_weighted_p_norm_average(out.viscosities[i],volume_fractions, composition_viscosities, composition_dviscosities_dpressure, viscosity_averaging_p); + + } } out.densities[i] = density; @@ -432,7 +512,7 @@ namespace aspect double e_ii = 0.; if (use_strain_weakening == true && use_finite_strain_tensor == false && this->get_timestep_number() > 0) { - edot_ii = std::max(sqrt(std::fabs(second_invariant(deviator(in.strain_rate[i])))),min_strain_rate); + edot_ii = std::max(sqrt(std::fabs(second_invariant(deviator(strain_rate)))),min_strain_rate); e_ii = edot_ii*this->get_timestep(); // Update reaction term out.reaction_terms[i][0] = e_ii; @@ -739,7 +819,6 @@ namespace aspect { prm.enter_subsection ("Visco Plastic"); { - // Reference and minimum/maximum values reference_T = prm.get_double("Reference temperature"); min_strain_rate = prm.get_double("Minimum strain rate"); @@ -770,6 +849,7 @@ namespace aspect if (use_strain_weakening) AssertThrow(this->n_compositional_fields() >= 1, ExcMessage("There must be at least one compositional field. ")); + use_finite_strain_tensor = prm.get_bool ("Use finite strain tensor"); if (use_finite_strain_tensor) AssertThrow(this->n_compositional_fields() >= s, diff --git a/source/utilities.cc b/source/utilities.cc index 49e27e111f5..86b8b2f09d9 100644 --- a/source/utilities.cc +++ b/source/utilities.cc @@ -2536,6 +2536,21 @@ namespace aspect } + template + SymmetricTensor<2,dim> + symmetric_independent_component_matrix (const unsigned int k) + { + // somehow I need this, because otherwise it complains that I am passing 3 arguments and it only takes 2... + const bool boolian = k < SymmetricTensor<2,dim>::n_independent_components; + Assert(boolian, ExcMessage("The component is larger then the amount of independent components in the matrix.") ); + + const TableIndices<2> indices_ij = SymmetricTensor<2,dim>::unrolled_to_component_indices (k); + + Tensor<2,dim> result; + result[indices_ij] = 1; + + return symmetrize(result); + } // Explicit instantiations @@ -2622,5 +2637,8 @@ namespace aspect template std_cxx11::array convert_point_to_array<2>(const Point<2> &point); template std_cxx11::array convert_point_to_array<3>(const Point<3> &point); + template SymmetricTensor<2,2> symmetric_independent_component_matrix (const unsigned int k); + template SymmetricTensor<2,3> symmetric_independent_component_matrix (const unsigned int k); + } } diff --git a/tests/visco_plastic_derivatives.cc b/tests/visco_plastic_derivatives.cc new file mode 100644 index 00000000000..a9305171aac --- /dev/null +++ b/tests/visco_plastic_derivatives.cc @@ -0,0 +1,270 @@ +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +template +void f(const aspect::SimulatorAccess &simulator_access, + aspect::Assemblers::Manager &, + unsigned int parameter) +{ + // Prepare NewtonHandler for test + const_cast &> (simulator_access.get_newton_handler()).set_newton_derivative_scaling_factor(1.0); + + std::cout << std::endl << "Testing ViscoPlastic derivatives against analytical derivatives" << std::endl; + + using namespace aspect::MaterialModel; + MaterialModelInputs in_base(5,3); + in_base.composition[0][0] = 0; + in_base.composition[0][1] = 0; + in_base.composition[0][2] = 0; + in_base.composition[1][0] = 0.75; + in_base.composition[1][1] = 0.15; + in_base.composition[1][2] = 0.10; + in_base.composition[2][0] = 0; + in_base.composition[2][1] = 0.2; + in_base.composition[2][2] = 0.4; + in_base.composition[3][0] = 0; + in_base.composition[3][1] = 0.2; + in_base.composition[3][2] = 0.4; + in_base.composition[4][0] = 1; + in_base.composition[4][1] = 0; + in_base.composition[4][2] = 0; + + in_base.pressure[0] = 1e9; + in_base.pressure[1] = 5e9; + in_base.pressure[2] = 2e10; + in_base.pressure[3] = 2e11; + in_base.pressure[4] = 5e8; + + /** + * We can't take too small strain-rates, because then the difference in the + * viscosity will be too small for the double accuracy which stores + * the viscosity solutions and the finite difference solution. + */ + in_base.strain_rate[0] = SymmetricTensor<2,dim>(); + in_base.strain_rate[0][0][0] = 1e-12; + in_base.strain_rate[0][0][1] = 1e-12; + in_base.strain_rate[0][1][1] = 1e-11; + in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[1][0][0] = -1.71266e-13; + in_base.strain_rate[1][0][1] = -5.82647e-12; + in_base.strain_rate[1][1][1] = 4.21668e-14; + in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[2][1][1] = 1e-13; + in_base.strain_rate[2][0][1] = 1e-11; + in_base.strain_rate[2][0][0] = -1e-12; + in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[3][1][1] = 4.9e-21; + in_base.strain_rate[3][0][1] = 4.9e-21; + in_base.strain_rate[3][0][0] = 4.9e-21; + in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[4][1][1] = 1e-11; + in_base.strain_rate[4][0][1] = 1e-11; + in_base.strain_rate[4][0][0] = -1e-11; + + in_base.temperature[0] = 293; + in_base.temperature[1] = 1600; + in_base.temperature[2] = 2000; + in_base.temperature[3] = 2100; + in_base.temperature[4] = 600; + + SymmetricTensor<2,dim> zerozero = SymmetricTensor<2,dim>(); + SymmetricTensor<2,dim> onezero = SymmetricTensor<2,dim>(); + SymmetricTensor<2,dim> oneone = SymmetricTensor<2,dim>(); + + zerozero[0][0] = 1; + onezero[1][0] = 0.5; // because symmetry doubles this entry + oneone[1][1] = 1; + + double finite_difference_accuracy = 1e-7; + double finite_difference_factor = 1+finite_difference_accuracy; + + bool Error = false; + + MaterialModelInputs in_dviscositydpressure(in_base); + in_dviscositydpressure.pressure[0] *= finite_difference_factor; + in_dviscositydpressure.pressure[1] *= finite_difference_factor; + in_dviscositydpressure.pressure[2] *= finite_difference_factor; + in_dviscositydpressure.pressure[3] *= finite_difference_factor; + in_dviscositydpressure.pressure[4] *= finite_difference_factor; + + MaterialModelInputs in_dviscositydstrainrate_zerozero(in_base); + MaterialModelInputs in_dviscositydstrainrate_onezero(in_base); + MaterialModelInputs in_dviscositydstrainrate_oneone(in_base); + + in_dviscositydstrainrate_zerozero.strain_rate[0] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[0][0][0]) * finite_difference_accuracy * zerozero; + in_dviscositydstrainrate_zerozero.strain_rate[1] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[1][0][0]) * finite_difference_accuracy * zerozero; + in_dviscositydstrainrate_zerozero.strain_rate[2] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[2][0][0]) * finite_difference_accuracy * zerozero; + in_dviscositydstrainrate_zerozero.strain_rate[3] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[3][0][0]) * finite_difference_accuracy * zerozero; + in_dviscositydstrainrate_zerozero.strain_rate[4] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[4][0][0]) * finite_difference_accuracy * zerozero; + in_dviscositydstrainrate_onezero.strain_rate[0] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[0][1][0]) * finite_difference_accuracy * onezero; + in_dviscositydstrainrate_onezero.strain_rate[1] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[1][1][0]) * finite_difference_accuracy * onezero; + in_dviscositydstrainrate_onezero.strain_rate[2] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[2][1][0]) * finite_difference_accuracy * onezero; + in_dviscositydstrainrate_onezero.strain_rate[3] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[3][1][0]) * finite_difference_accuracy * onezero; + in_dviscositydstrainrate_onezero.strain_rate[4] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[4][1][0]) * finite_difference_accuracy * onezero; + in_dviscositydstrainrate_oneone.strain_rate[0] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[0][1][1]) * finite_difference_accuracy * oneone; + in_dviscositydstrainrate_oneone.strain_rate[1] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[1][1][1]) * finite_difference_accuracy * oneone; + in_dviscositydstrainrate_oneone.strain_rate[2] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[2][1][1]) * finite_difference_accuracy * oneone; + in_dviscositydstrainrate_oneone.strain_rate[3] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[3][1][1]) * finite_difference_accuracy * oneone; + in_dviscositydstrainrate_oneone.strain_rate[4] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[4][1][1]) * finite_difference_accuracy * oneone; + + MaterialModelInputs in_dviscositydtemperature(in_base); + in_dviscositydtemperature.temperature[0] *= 1.0000000001; + in_dviscositydtemperature.temperature[1] *= 1.0000000001; + in_dviscositydtemperature.temperature[2] *= 1.0000000001; + in_dviscositydtemperature.temperature[3] *= 1.0000000001; + in_dviscositydtemperature.temperature[4] *= 1.0000000001; + + + MaterialModelOutputs out_base(5,3); + + MaterialModelOutputs out_dviscositydpressure(5,3); + MaterialModelOutputs out_dviscositydstrainrate_zerozero(5,3); + MaterialModelOutputs out_dviscositydstrainrate_onezero(5,3); + MaterialModelOutputs out_dviscositydstrainrate_oneone(5,3); + MaterialModelOutputs out_dviscositydtemperature(5,3); + + out_base.additional_outputs.push_back(std_cxx11::make_shared > (5)); + + simulator_access.get_material_model().evaluate(in_base, out_base); + simulator_access.get_material_model().evaluate(in_dviscositydpressure, out_dviscositydpressure); + simulator_access.get_material_model().evaluate(in_dviscositydstrainrate_zerozero, out_dviscositydstrainrate_zerozero); + simulator_access.get_material_model().evaluate(in_dviscositydstrainrate_onezero, out_dviscositydstrainrate_onezero); + simulator_access.get_material_model().evaluate(in_dviscositydstrainrate_oneone, out_dviscositydstrainrate_oneone); + simulator_access.get_material_model().evaluate(in_dviscositydtemperature, out_dviscositydtemperature); + + // set up additional output for the derivatives + MaterialModelDerivatives *derivatives; + derivatives = out_base.template get_additional_output >(); + + double temp; + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = (out_dviscositydpressure.viscosities[i] - out_base.viscosities[i]); + if (in_base.pressure[i] != 0) + { + temp /= (in_base.pressure[i] * finite_difference_accuracy); + } + std::cout << "pressure at point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_pressure[i] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_pressure[i]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_pressure[i]))) + // if (temp > derivatives->viscosity_derivative_wrt_pressure[i] * finite_difference_factor || temp < derivatives->viscosity_derivative_wrt_pressure[i] * (2-finite_difference_factor)) + { + std::cout << " Error: The derivative of the viscosity to the pressure is too different from the analitical value." << std::endl; + Error = true; + } + + } + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = out_dviscositydstrainrate_zerozero.viscosities[i] - out_base.viscosities[i]; + if (temp != 0) + { + temp /= std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[i][0][0]) * finite_difference_accuracy; + } + std::cout << "zerozero at point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][0][0] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][0][0]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][0][0]))) + { + std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analitical value." << std::endl; + Error = true; + } + + + + } + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = out_dviscositydstrainrate_onezero.viscosities[i] - out_base.viscosities[i]; + if (temp != 0) + { + temp /= std::fabs(in_dviscositydstrainrate_onezero.strain_rate[i][1][0]) * finite_difference_accuracy; + } + std::cout << "onezero at point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][1][0] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][1][0]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][1][0])) ) + { + std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; + Error = true; + } + } + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = out_dviscositydstrainrate_oneone.viscosities[i] - out_base.viscosities[i]; + if (temp != 0) + { + temp /= std::fabs(in_dviscositydstrainrate_oneone.strain_rate[i][1][1]) * finite_difference_accuracy; + } + std::cout << "oneone at point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][1][1] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][1][1]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][1][1])) ) + { + std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; + Error = true; + } + + } + + if (Error) + { + std::cout << "Some parts of the test where not successful." << std::endl; + } + else + { + std::cout << "OK" << std::endl; + } + + // Restore NewtonHandler state after test + const_cast &> (simulator_access.get_newton_handler()).set_newton_derivative_scaling_factor(0.0); + +} + +template <> +void f(const aspect::SimulatorAccess<3> &simulator_access, + aspect::Assemblers::Manager<3> &, + unsigned int parameter) +{ + AssertThrow(false,dealii::ExcInternalError()); +} + +template +void signal_connector (aspect::SimulatorSignals &signals) +{ + using namespace dealii; + std::cout << "* Connecting signals" << std::endl; + signals.set_assemblers.connect (std_cxx11::bind(&f, + std_cxx11::_1, + std_cxx11::_2, + 1)); + + signals.set_assemblers.connect (std_cxx11::bind(&f, + std_cxx11::_1, + std_cxx11::_2, + 2)); +} + +ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, + signal_connector<3>) + diff --git a/tests/visco_plastic_derivatives.prm b/tests/visco_plastic_derivatives.prm new file mode 100644 index 00000000000..90cd29400e6 --- /dev/null +++ b/tests/visco_plastic_derivatives.prm @@ -0,0 +1,107 @@ +# This test checks wether the derivatives of the drucker prager +# material model are the same as the finite difference derivatives +# created by the drucker prager material model. +set Additional shared libraries = tests/libvisco_plastic_derivatives.so + +set Dimension = 2 +set End time = 0 +set Use years in output instead of seconds = true +set Nonlinear solver scheme = Newton Stokes +set Max nonlinear iterations = 1 +set Number of cheap Stokes solver steps = 0 + +# Model geometry (100x100 km, 10 km spacing) +subsection Geometry model + set Model name = box + subsection Box + set X repetitions = 10 + set Y repetitions = 10 + set X extent = 100e3 + set Y extent = 100e3 + end +end + +# Mesh refinement specifications +subsection Mesh refinement + set Initial adaptive refinement = 0 + set Initial global refinement = 0 + set Time steps between mesh refinement = 0 +end + + +# Boundary classifications (fixed T boundaries, prescribed velocity) +subsection Model settings + set Fixed temperature boundary indicators = bottom, top, left, right + set Prescribed velocity boundary indicators = bottom y: function, top y: function, left x: function, right x: function +end + +# Velocity on boundaries characterized by functions +subsection Boundary velocity model + subsection Function + set Variable names = x,y + set Function constants = m=0.0005, year=1 + set Function expression = if (x<50e3 , -1*m/year, 1*m/year); if (y<50e3 , 1*m/year, -1*m/year); + end +end + +# Temperature boundary and initial conditions +subsection Boundary temperature model + set List of model names = box + subsection Box + set Bottom temperature = 273 + set Left temperature = 273 + set Right temperature = 273 + set Top temperature = 273 + end +end +subsection Initial temperature model + set Model name = function + subsection Function + set Function expression = 273 + end +end + +# Compositional fields used to track finite strain invariant +subsection Compositional fields + set Number of fields = 3 +end + +# We prescribe some initial strain at the center of the domain +subsection Initial composition model + set Model name = function + subsection Function + set Variable names = x,y + set Function expression = if(x>=45e3&x<=55e3&y>=45.3e3&y<=55.e3,0.2,0);0;0 + end +end + +# Boundary composition specification +subsection Boundary composition model + set Model name = initial composition +end + +# Material model (values for background material) +subsection Material model + set Model name = visco plastic + subsection Visco Plastic + set Angles of internal friction = 30. + end +end + +# Gravity model +subsection Gravity model + set Model name = vertical + subsection Vertical + set Magnitude = 10.0 + end +end + +# Post processing +# named additional outputs includes the weakened cohesions and friction angles +subsection Postprocess + set List of postprocessors = velocity statistics, mass flux statistics, visualization + subsection Visualization + set List of output variables = viscosity, strain rate, named additional outputs + set Output format = gnuplot + end +end diff --git a/tests/visco_plastic_derivatives/screen-output b/tests/visco_plastic_derivatives/screen-output new file mode 100644 index 00000000000..b6d38a72207 --- /dev/null +++ b/tests/visco_plastic_derivatives/screen-output @@ -0,0 +1,44 @@ + +Loading shared library <./libvisco_plastic_derivatives.so> + +* Connecting signals + +Testing ViscoPlastic derivatives against analytical derivatives +pressure at point 0: Finite difference = 0. Analytical derivative = 0 +pressure at point 1: Finite difference = 2.12594e+08. Analytical derivative = 2.12594e+08 +pressure at point 2: Finite difference = 2.68615e+08. Analytical derivative = 2.68615e+08 +pressure at point 3: Finite difference = 0. Analytical derivative = 0 +pressure at point 4: Finite difference = 1.8585e+16. Analytical derivative = 1.8585e+16 +zerozero at point 0: Finite difference = 0. Analytical derivative = 0 +zerozero at point 1: Finite difference = 6.58655e+26. Analytical derivative = 6.58655e+26 +zerozero at point 2: Finite difference = 1.23761e+27. Analytical derivative = 1.23761e+27 +zerozero at point 3: Finite difference = 0. Analytical derivative = 0 +zerozero at point 4: Finite difference = 4.13904e+35. Analytical derivative = 4.13904e+35 +onezero at point 0: Finite difference = 0. Analytical derivative = 0 +onezero at point 1: Finite difference = 3.5961e+28. Analytical derivative = 3.5961e+28 +onezero at point 2: Finite difference = -2.2502e+28. Analytical derivative = -2.2502e+28 +onezero at point 3: Finite difference = 0. Analytical derivative = 0 +onezero at point 4: Finite difference = -4.13904e+35. Analytical derivative = -4.13904e+35 +oneone at point 0: Finite difference = 0. Analytical derivative = 0 +oneone at point 1: Finite difference = -6.58657e+26. Analytical derivative = -6.58657e+26 +oneone at point 2: Finite difference = -1.23761e+27. Analytical derivative = -1.23761e+27 +oneone at point 3: Finite difference = 0. Analytical derivative = 0 +oneone at point 4: Finite difference = -4.13904e+35. Analytical derivative = -4.13904e+35 +OK +Number of active cells: 100 (on 1 levels) +Number of degrees of freedom: 2,767 (882+121+441+441+441+441) + +*** Timestep 0: t=0 years + + Postprocessing: + RMS, max velocity: 0 m/year, 0 m/year + Mass fluxes through boundary parts: 0 kg/yr, 0 kg/yr, 0 kg/yr, 0 kg/yr + Writing graphical output: output-visco_plastic_derivatives/solution/solution-00000 + +Termination requested by criterion: end time + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ +