Skip to content

Commit cbd582f

Browse files
SachinJalanSachinJalan
committed
Test File Added
1 parent 68caa89 commit cbd582f

2 files changed

Lines changed: 202 additions & 8 deletions

File tree

‎diffmpm/material.py‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -203,18 +203,18 @@ def thermodynamic_pressure(self, volumetric_strain):
203203
def compute_stress(self, dstrain, particles, state_vars):
204204
shear_rate_threshold = 1e-15
205205
dirac_delta = jnp.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0]).reshape((6, 1))
206-
lax.cond(
206+
dirac_delta=lax.cond(
207207
self.ndim == 1,
208208
lambda x: x.at[0, 0].set(1.0),
209209
lambda x: x.at[0:2, 0].set(1.0),
210210
dirac_delta,
211211
)
212-
lax.select(
212+
self.properties["critical_shear_rate"] = lax.select(
213213
self.properties["critical_shear_rate"] < shear_rate_threshold,
214214
shear_rate_threshold,
215215
self.properties["critical_shear_rate"],
216216
)
217-
217+
@jit
218218
def compute_stress_per_particle(
219219
particle_strain_rate,
220220
self,
@@ -227,17 +227,17 @@ def compute_stress_per_particle(
227227
shear_rate = jnp.sqrt(
228228
2.0 * (strain_r.T @ (strain_r) + strain_r[-3:].T @ strain_r[-3:])
229229
)
230-
apparent_velocity_true = 2.0 * (
230+
apparent_viscosity_true = 2.0 * (
231231
(self.properties["tau0"] / shear_rate[0, 0]) + self.properties["mu"]
232232
)
233233
condition = (shear_rate[0, 0] * shear_rate[0, 0]) > (
234234
self.properties["critical_shear_rate"]
235235
* self.properties["critical_shear_rate"]
236236
)
237-
apparent_viscosity = lax.select(condition, apparent_velocity_true, 0.0)
237+
apparent_viscosity = lax.select(condition, apparent_viscosity_true, 0.0)
238238
tau = apparent_viscosity * strain_r
239239
trace_invariant2 = 0.5 * jnp.dot(tau[:3, 0], tau[:3, 0])
240-
lax.cond(
240+
tau=lax.cond(
241241
trace_invariant2 < (self.properties["tau0"] * self.properties["tau0"]),
242242
lambda x: x.at[:].set(0),
243243
lambda x: x,
@@ -254,8 +254,8 @@ def compute_stress_per_particle(
254254
)
255255
return updated_stress_per_particle, state_vars_pressure
256256

257-
updated_stress, state_vars["pressure"] = jit(
258-
vmap(compute_stress_per_particle, in_axes=(0, None, 0, 0, None))
257+
updated_stress, state_vars["pressure"] = vmap(
258+
compute_stress_per_particle, in_axes=(0, None, 0, 0, None)
259259
)(
260260
particles.strain_rate,
261261
self,

‎tests/test_bingham.py‎

Lines changed: 194 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,194 @@
1+
import pytest
2+
import jax.numpy as jnp
3+
from diffmpm.material import Bingham
4+
from diffmpm.particle import Particles
5+
from diffmpm.element import Quadrilateral4Node
6+
from diffmpm.constraint import Constraint
7+
from diffmpm.node import Nodes
8+
9+
material_dstrain_particles_state_vars_element_targets_dt = [
10+
(
11+
Bingham(
12+
{
13+
"density": 1000,
14+
"youngs_modulus": 1.0e7,
15+
"poisson_ratio": 0.3,
16+
"tau0": 771.8,
17+
"mu": 0.0451,
18+
"critical_shear_rate": 0.2,
19+
"ndim": 2,
20+
}
21+
),
22+
None,
23+
Particles(
24+
jnp.array([[0.5, 0.5]]).reshape(1, 1, 2),
25+
(
26+
Bingham(
27+
{
28+
"density": 1000,
29+
"youngs_modulus": 1.0e7,
30+
"poisson_ratio": 0.3,
31+
"tau0": 771.8,
32+
"mu": 0.0451,
33+
"critical_shear_rate": 0.2,
34+
"ndim": 2,
35+
}
36+
)
37+
),
38+
jnp.array([0]),
39+
),
40+
{"pressure": jnp.zeros(1)},
41+
Quadrilateral4Node(
42+
(1, 1),
43+
1,
44+
(4.0, 4.0),
45+
[],
46+
Nodes(4, jnp.array([-2, -2, 2, -2, -2, 2, 2, 2]).reshape((4, 1, 2))),
47+
),
48+
jnp.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0]).reshape((6, 1)),
49+
1.0,
50+
),
51+
(
52+
Bingham(
53+
{
54+
"density": 1000,
55+
"youngs_modulus": 1.0e7,
56+
"poisson_ratio": 0.3,
57+
"tau0": 771.8,
58+
"mu": 0.0451,
59+
"critical_shear_rate": 0.2,
60+
"ndim": 2,
61+
}
62+
),
63+
None,
64+
Particles(
65+
jnp.array([[0.5, 0.5]]).reshape(1, 1, 2),
66+
(
67+
Bingham(
68+
{
69+
"density": 1000,
70+
"youngs_modulus": 1.0e7,
71+
"poisson_ratio": 0.3,
72+
"tau0": 771.8,
73+
"mu": 0.0451,
74+
"critical_shear_rate": 0.2,
75+
"ndim": 2,
76+
}
77+
)
78+
),
79+
jnp.array([0]),
80+
),
81+
{"pressure": jnp.zeros(1)},
82+
Quadrilateral4Node(
83+
(1, 1),
84+
1,
85+
(4.0, 4.0),
86+
[(0, Constraint(0, 0.02)), (0, Constraint(1, 0.03))],
87+
Nodes(4, jnp.array([-2, -2, 2, -2, -2, 2, 2, 2]).reshape((4, 1, 2))),
88+
),
89+
jnp.array([-52083.3333333333, -52083.3333333333, 0.0, 0.0, 0.0, 0.0]).reshape(
90+
(6, 1)
91+
),
92+
1.0,
93+
),
94+
(
95+
Bingham(
96+
{
97+
"density": 1000,
98+
"youngs_modulus": 1.0e7,
99+
"poisson_ratio": 0.3,
100+
"tau0": 200.0,
101+
"mu": 200.0,
102+
"critical_shear_rate": 0.2,
103+
"ndim": 2,
104+
}
105+
),
106+
None,
107+
Particles(
108+
jnp.array([[0.5, 0.5]]).reshape(1, 1, 2),
109+
(
110+
Bingham(
111+
{
112+
"density": 1000,
113+
"youngs_modulus": 1.0e7,
114+
"poisson_ratio": 0.3,
115+
"tau0": 200.0,
116+
"mu": 200.0,
117+
"critical_shear_rate": 0.2,
118+
"ndim": 2,
119+
}
120+
)
121+
),
122+
jnp.array([0]),
123+
),
124+
{"pressure": jnp.zeros(1)},
125+
Quadrilateral4Node(
126+
(1, 1),
127+
1,
128+
(4.0, 4.0),
129+
[(0, Constraint(0, 2.0)), (0, Constraint(1, 3.0))],
130+
Nodes(4, jnp.array([-2, -2, 2, -2, -2, 2, 2, 2]).reshape((4, 1, 2))),
131+
),
132+
jnp.array(
133+
[-5208520.35574006, -5208613.86694342, 0.0, -233.778008402801, 0.0, 0.0]
134+
).reshape((6, 1)),
135+
1.0,
136+
),
137+
(
138+
Bingham(
139+
{
140+
"density": 1000,
141+
"youngs_modulus": 1.0e7,
142+
"poisson_ratio": 0.3,
143+
"tau0": 200.0,
144+
"mu": 200.0,
145+
"critical_shear_rate": 0.2,
146+
"ndim": 2,
147+
"incompressible": True,
148+
}
149+
),
150+
None,
151+
Particles(
152+
jnp.array([[0.5, 0.5]]).reshape(1, 1, 2),
153+
(
154+
Bingham(
155+
{
156+
"density": 1000,
157+
"youngs_modulus": 1.0e7,
158+
"poisson_ratio": 0.3,
159+
"tau0": 200.0,
160+
"mu": 200.0,
161+
"critical_shear_rate": 0.2,
162+
"ndim": 2,
163+
"incompressible": True,
164+
}
165+
)
166+
),
167+
jnp.array([0]),
168+
),
169+
{"pressure": jnp.zeros(1)},
170+
Quadrilateral4Node(
171+
(1, 1),
172+
1,
173+
(4.0, 4.0),
174+
[(0, Constraint(0, 2.0)), (0, Constraint(1, 3.0))],
175+
Nodes(4, jnp.array([-2, -2, 2, -2, -2, 2, 2, 2]).reshape((4, 1, 2))),
176+
),
177+
jnp.array(
178+
[-187.0224067222, -280.5336100834, 0.0, -233.778008402801, 0.0, 0.0]
179+
).reshape((6, 1)),
180+
1.0,
181+
),
182+
]
183+
184+
185+
@pytest.mark.parametrize(
186+
"material, dstrain, particles, state_vars, element, target, dt",
187+
material_dstrain_particles_state_vars_element_targets_dt,
188+
)
189+
def test_compute_stress(material, dstrain, particles, state_vars, element, target, dt):
190+
if element.constraints != []:
191+
element.apply_boundary_constraints()
192+
particles.compute_strain(element, dt)
193+
stress = material.compute_stress(dstrain, particles, state_vars)
194+
assert jnp.allclose(stress, target)

0 commit comments

Comments
 (0)