Skip to content

Commit 3df9c3a

Browse files
MFratersMFraters
committed
Add tests for derivatives of visco-plastic material model derivatives.
1 parent 8793a39 commit 3df9c3a

File tree

3 files changed

+374
-0
lines changed

3 files changed

+374
-0
lines changed

‎tests/visco_plastic_derivatives.cc‎

Lines changed: 274 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,274 @@
1+
#include <aspect/simulator.h>
2+
#include <deal.II/grid/tria.h>
3+
#include <aspect/material_model/interface.h>
4+
#include <aspect/material_model/visco_plastic.h>
5+
#include <aspect/simulator_access.h>
6+
#include <aspect/newton.h>
7+
8+
#include <iostream>
9+
10+
int f(double parameter)
11+
{
12+
13+
std::cout << std::endl << "Test for p = " << parameter << std::endl;
14+
15+
const int dim=2;
16+
using namespace aspect::MaterialModel;
17+
MaterialModelInputs<dim> in_base(5,3);
18+
in_base.composition[0][0] = 0;
19+
in_base.composition[0][1] = 0;
20+
in_base.composition[0][2] = 0;
21+
in_base.composition[1][0] = 0.75;
22+
in_base.composition[1][1] = 0.15;
23+
in_base.composition[1][2] = 0.10;
24+
in_base.composition[2][0] = 0;
25+
in_base.composition[2][1] = 0.2;
26+
in_base.composition[2][2] = 0.4;
27+
in_base.composition[3][0] = 0;
28+
in_base.composition[3][1] = 0.2;
29+
in_base.composition[3][2] = 0.4;
30+
in_base.composition[4][0] = 1;
31+
in_base.composition[4][1] = 0;
32+
in_base.composition[4][2] = 0;
33+
34+
in_base.pressure[0] = 1e9;
35+
in_base.pressure[1] = 5e9;
36+
in_base.pressure[2] = 2e10;
37+
in_base.pressure[3] = 2e11;
38+
in_base.pressure[4] = 5e8;
39+
40+
/**
41+
* We can't take to small strain-rates, because then the difference in the
42+
* visocisty will be too small for the double accuracy which stores
43+
* the visocity solutions and the finite diference solution.
44+
*/
45+
in_base.strain_rate[0] = SymmetricTensor<2,dim>();
46+
in_base.strain_rate[0][0][0] = 1e-12;
47+
in_base.strain_rate[0][0][1] = 1e-12;
48+
in_base.strain_rate[0][1][1] = 1e-11;
49+
in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]);
50+
in_base.strain_rate[1][0][0] = -1.71266e-13;
51+
in_base.strain_rate[1][0][1] = -5.82647e-12;
52+
in_base.strain_rate[1][1][1] = 4.21668e-14;
53+
in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]);
54+
in_base.strain_rate[2][1][1] = 1e-13;
55+
in_base.strain_rate[2][0][1] = 1e-11;
56+
in_base.strain_rate[2][0][0] = -1e-12;
57+
in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]);
58+
in_base.strain_rate[3][1][1] = 4.9e-21;
59+
in_base.strain_rate[3][0][1] = 4.9e-21;
60+
in_base.strain_rate[3][0][0] = 4.9e-21;
61+
in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]);
62+
in_base.strain_rate[4][1][1] = 1e-11;
63+
in_base.strain_rate[4][0][1] = 1e-11;
64+
in_base.strain_rate[4][0][0] = -1e-11;
65+
66+
in_base.temperature[0] = 293;
67+
in_base.temperature[1] = 1600;
68+
in_base.temperature[2] = 2000;
69+
in_base.temperature[3] = 2100;
70+
in_base.temperature[4] = 600;
71+
72+
SymmetricTensor<2,dim> zerozero = SymmetricTensor<2,dim>();
73+
SymmetricTensor<2,dim> onezero = SymmetricTensor<2,dim>();
74+
SymmetricTensor<2,dim> oneone = SymmetricTensor<2,dim>();
75+
76+
zerozero[0][0] = 1;
77+
onezero[1][0] = 0.5; // because symmetry doubles this entry
78+
oneone[1][1] = 1;
79+
80+
double finite_difference_accuracy = 1e-7;
81+
double finite_difference_factor = 1+finite_difference_accuracy;
82+
83+
bool Error = false;
84+
85+
MaterialModelInputs<dim> in_dviscositydpressure(in_base);
86+
in_dviscositydpressure.pressure[0] *= finite_difference_factor;
87+
in_dviscositydpressure.pressure[1] *= finite_difference_factor;
88+
in_dviscositydpressure.pressure[2] *= finite_difference_factor;
89+
in_dviscositydpressure.pressure[3] *= finite_difference_factor;
90+
in_dviscositydpressure.pressure[4] *= finite_difference_factor;
91+
92+
MaterialModelInputs<dim> in_dviscositydstrainrate_zerozero(in_base);
93+
MaterialModelInputs<dim> in_dviscositydstrainrate_onezero(in_base);
94+
MaterialModelInputs<dim> in_dviscositydstrainrate_oneone(in_base);
95+
96+
in_dviscositydstrainrate_zerozero.strain_rate[0] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[0][0][0]) * finite_difference_accuracy * zerozero;
97+
in_dviscositydstrainrate_zerozero.strain_rate[1] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[1][0][0]) * finite_difference_accuracy * zerozero;
98+
in_dviscositydstrainrate_zerozero.strain_rate[2] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[2][0][0]) * finite_difference_accuracy * zerozero;
99+
in_dviscositydstrainrate_zerozero.strain_rate[3] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[3][0][0]) * finite_difference_accuracy * zerozero;
100+
in_dviscositydstrainrate_zerozero.strain_rate[4] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[4][0][0]) * finite_difference_accuracy * zerozero;
101+
in_dviscositydstrainrate_onezero.strain_rate[0] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[0][1][0]) * finite_difference_accuracy * onezero;
102+
in_dviscositydstrainrate_onezero.strain_rate[1] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[1][1][0]) * finite_difference_accuracy * onezero;
103+
in_dviscositydstrainrate_onezero.strain_rate[2] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[2][1][0]) * finite_difference_accuracy * onezero;
104+
in_dviscositydstrainrate_onezero.strain_rate[3] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[3][1][0]) * finite_difference_accuracy * onezero;
105+
in_dviscositydstrainrate_onezero.strain_rate[4] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[4][1][0]) * finite_difference_accuracy * onezero;
106+
in_dviscositydstrainrate_oneone.strain_rate[0] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[0][1][1]) * finite_difference_accuracy * oneone;
107+
in_dviscositydstrainrate_oneone.strain_rate[1] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[1][1][1]) * finite_difference_accuracy * oneone;
108+
in_dviscositydstrainrate_oneone.strain_rate[2] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[2][1][1]) * finite_difference_accuracy * oneone;
109+
in_dviscositydstrainrate_oneone.strain_rate[3] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[3][1][1]) * finite_difference_accuracy * oneone;
110+
in_dviscositydstrainrate_oneone.strain_rate[4] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[4][1][1]) * finite_difference_accuracy * oneone;
111+
112+
MaterialModelInputs<dim> in_dviscositydtemperature(in_base);
113+
in_dviscositydtemperature.temperature[0] *= 1.0000000001;
114+
in_dviscositydtemperature.temperature[1] *= 1.0000000001;
115+
in_dviscositydtemperature.temperature[2] *= 1.0000000001;
116+
in_dviscositydtemperature.temperature[3] *= 1.0000000001;
117+
in_dviscositydtemperature.temperature[4] *= 1.0000000001;
118+
119+
120+
MaterialModelOutputs<dim> out_base(5,3);
121+
122+
MaterialModelOutputs<dim> out_dviscositydpressure(5,3);
123+
MaterialModelOutputs<dim> out_dviscositydstrainrate_zerozero(5,3);
124+
MaterialModelOutputs<dim> out_dviscositydstrainrate_onezero(5,3);
125+
MaterialModelOutputs<dim> out_dviscositydstrainrate_oneone(5,3);
126+
MaterialModelOutputs<dim> out_dviscositydtemperature(5,3);
127+
128+
if (out_base.get_additional_output<MaterialModelDerivatives<dim> >() != NULL)
129+
throw "error";
130+
131+
out_base.additional_outputs.push_back(std::make_shared<MaterialModelDerivatives<dim> > (5));
132+
133+
ViscoPlastic<dim> mat;
134+
ParameterHandler prm;
135+
mat.declare_parameters(prm);
136+
137+
prm.enter_subsection("Compositional fields");
138+
{
139+
prm.set ("Number of fields", "3");
140+
}
141+
prm.leave_subsection();
142+
143+
prm.enter_subsection("Material model");
144+
{
145+
prm.enter_subsection ("Visco Plastic");
146+
{
147+
// prm.enter_subsection ("Viscosity");
148+
// {
149+
//prm.set ("Reference strain rate", "1e-20");
150+
prm.set ("Angles of internal friction", "30");
151+
// }
152+
// prm.leave_subsection();
153+
}
154+
prm.leave_subsection();
155+
}
156+
prm.leave_subsection();
157+
158+
mat.parse_parameters(prm);
159+
160+
mat.evaluate(in_base, out_base);
161+
mat.evaluate(in_dviscositydpressure, out_dviscositydpressure);
162+
mat.evaluate(in_dviscositydstrainrate_zerozero, out_dviscositydstrainrate_zerozero);
163+
mat.evaluate(in_dviscositydstrainrate_onezero, out_dviscositydstrainrate_onezero);
164+
mat.evaluate(in_dviscositydstrainrate_oneone, out_dviscositydstrainrate_oneone);
165+
mat.evaluate(in_dviscositydtemperature, out_dviscositydtemperature);
166+
167+
// set up additional output for the derivatives
168+
MaterialModelDerivatives<dim> *derivatives;
169+
derivatives = out_base.get_additional_output<MaterialModelDerivatives<dim> >();
170+
171+
double temp;
172+
for (unsigned int i = 0; i < 5; i++)
173+
{
174+
// prevent division by zero. If it is zero, the test has passed, because or
175+
// the finite difference and the analytical result match perfectly, or (more
176+
// likely) the material model in independent of this variable.
177+
temp = (out_dviscositydpressure.viscosities[i] - out_base.viscosities[i]);
178+
if (in_base.pressure[i] != 0)
179+
{
180+
temp /= (in_base.pressure[i] * finite_difference_accuracy);
181+
}
182+
std::cout << "pressure at point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_pressure[i] << std::endl;
183+
if (std::fabs(temp - derivatives->viscosity_derivative_wrt_pressure[i]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_pressure[i])))
184+
// if (temp > derivatives->viscosity_derivative_wrt_pressure[i] * finite_difference_factor || temp < derivatives->viscosity_derivative_wrt_pressure[i] * (2-finite_difference_factor))
185+
{
186+
std::cout << " Error: The derivative of the viscosity to the pressure is too different from the analitical value." << std::endl;
187+
Error = true;
188+
}
189+
190+
}
191+
192+
for (unsigned int i = 0; i < 5; i++)
193+
{
194+
// prevent division by zero. If it is zero, the test has passed, because or
195+
// the finite difference and the analytical result match perfectly, or (more
196+
// likely) the material model in independent of this variable.
197+
temp = out_dviscositydstrainrate_zerozero.viscosities[i] - out_base.viscosities[i];
198+
if (temp != 0)
199+
{
200+
temp /= std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[i][0][0]) * finite_difference_accuracy;
201+
}
202+
std::cout << "zerozero at point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][0][0] << std::endl;
203+
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])))
204+
{
205+
std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analitical value." << std::endl;
206+
Error = true;
207+
}
208+
209+
210+
211+
}
212+
213+
for (unsigned int i = 0; i < 5; i++)
214+
{
215+
// prevent division by zero. If it is zero, the test has passed, because or
216+
// the finite difference and the analytical result match perfectly, or (more
217+
// likely) the material model in independent of this variable.
218+
temp = out_dviscositydstrainrate_onezero.viscosities[i] - out_base.viscosities[i];
219+
if (temp != 0)
220+
{
221+
temp /= std::fabs(in_dviscositydstrainrate_onezero.strain_rate[i][1][0]) * finite_difference_accuracy;
222+
}
223+
std::cout << "onezero at point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][1][0] << std::endl;
224+
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])) )
225+
{
226+
std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analitical value." << std::endl;
227+
Error = true;
228+
}
229+
}
230+
231+
for (unsigned int i = 0; i < 5; i++)
232+
{
233+
// prevent division by zero. If it is zero, the test has passed, because or
234+
// the finite difference and the analytical result match perfectly, or (more
235+
// likely) the material model in independent of this variable.
236+
temp = out_dviscositydstrainrate_oneone.viscosities[i] - out_base.viscosities[i];
237+
if (temp != 0)
238+
{
239+
temp /= std::fabs(in_dviscositydstrainrate_oneone.strain_rate[i][1][1]) * finite_difference_accuracy;
240+
}
241+
std::cout << "oneone at point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][1][1] << std::endl;
242+
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])) )
243+
{
244+
std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analitical value." << std::endl;
245+
Error = true;
246+
}
247+
248+
}
249+
250+
if (Error)
251+
{
252+
std::cout << "Some parts of the test where not succesful." << std::endl;
253+
}
254+
else
255+
{
256+
std::cout << "OK" << std::endl;
257+
}
258+
259+
return 42;
260+
}
261+
262+
int exit_function()
263+
{
264+
exit(0);
265+
return 42;
266+
}
267+
// run this function by initializing a global variable by it
268+
int ik = f(-1); // Testing harmonic mean
269+
int ji = f(0); // Testing geometric mean
270+
int jj = f(1); // Testing arithmetic mean
271+
int kj = f(1000); // Testing max function
272+
int kl = exit_function();
273+
274+

‎tests/visco_plastic_derivatives.prm‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
# This test checks wether the derivatives of the drucker prager
2+
# material model are the same as the finite difference derivatives
3+
# created by the drucker prager material model.
4+
set Additional shared libraries = tests/libvisco_plastic_derivatives.so

0 commit comments

Comments
 (0)