Line data Source code
1 : #ifndef INTEGRATOR_BASE_MECHANICS_H
2 : #define INTEGRATOR_BASE_MECHANICS_H
3 :
4 : #include "AMReX.H"
5 : #include "BC/Operator/Elastic/Elastic.H"
6 : #include "BC/Operator/Elastic/Constant.H"
7 : #include "BC/Operator/Elastic/TensionTest.H"
8 : #include "BC/Operator/Elastic/Expression.H"
9 : #include "Integrator/Integrator.H"
10 : #include "Numeric/Stencil.H"
11 : #include "Model/Solid/Solid.H"
12 : #include "Solver/Nonlocal/Linear.H"
13 : #include "Solver/Nonlocal/Newton.H"
14 : #include "Operator/Operator.H"
15 : #include "IC/Constant.H"
16 : #include "IC/Expression.H"
17 :
18 : namespace Integrator
19 : {
20 : namespace Base
21 : {
22 : template<class MODEL>
23 : class Mechanics: virtual public Integrator
24 : {
25 : public:
26 :
27 : enum Type { Static, Dynamic, Disable };
28 :
29 131 : Mechanics() : Integrator()
30 25 : {}
31 :
32 : // The mechanics integrator manages the solution of an elastic
33 : // solve using the MLMG solver.
34 25 : static void Parse(Mechanics& value, IO::ParmParse& pp)
35 : {
36 : BL_PROFILE("Integrator::Base::Mechanics::Parse");
37 25 : if (pp.contains("type"))
38 : {
39 18 : std::string type_str;
40 : // Type of mecahnics to use.
41 : // Static: do full implicit solve.
42 : // Dynamic: evolve dynamic equations with explicit dynamics
43 : // Disable: do nothing.
44 9 : pp_query_validate("type", type_str, {"disable","static","dynamic"});
45 9 : if (type_str == "static") value.m_type = Type::Static;
46 3 : else if (type_str == "dynamic") value.m_type = Type::Dynamic;
47 3 : else if (type_str == "disable") value.m_type = Type::Disable;
48 0 : else Util::Abort(INFO, "Invalid type ", type_str, " specified");
49 : }
50 25 : if (value.m_type == Type::Disable) return;
51 :
52 : // Treat mechanics fields as changing in time. [false]
53 : // You should use this if you care about other physics driven by
54 : // the output of this integrator.
55 19 : pp_query_default("time_evolving", value.m_time_evolving,false);
56 :
57 19 : pp_query_default("plot_disp", value.plot_disp, true); // Include displacement field in output
58 19 : pp_query_default("plot_rhs", value.plot_rhs, true); // Include right-hand side in output
59 19 : pp_query_default("plot_psi", value.plot_psi, true); // Include :math:`\psi` field in output
60 19 : pp_query_default("plot_stress", value.plot_stress, true); // Include stress in output
61 19 : pp_query_default("plot_strain", value.plot_strain, true); // Include strain in output
62 :
63 19 : value.RegisterGeneralFab(value.disp_mf, 1, 2, value.plot_disp, "disp", value.m_time_evolving);
64 19 : value.RegisterGeneralFab(value.rhs_mf, 1, 2, value.plot_rhs, "rhs", value.m_time_evolving);
65 19 : value.RegisterGeneralFab(value.stress_mf, 1, 2, value.plot_stress, "stress", value.m_time_evolving);
66 19 : value.RegisterGeneralFab(value.strain_mf, 1, 2, value.plot_strain, "strain", value.m_time_evolving);
67 :
68 19 : if (value.m_type == Type::Static)
69 : {
70 : // Read parameters for :ref:`Solver::Nonlocal::Newton` solver
71 19 : pp_queryclass("solver", value.solver);
72 : }
73 19 : if (value.m_type == Type::Dynamic)
74 : {
75 0 : value.RegisterGeneralFab(value.vel_mf, 1, 2, "vel",true);
76 0 : value.RegisterGeneralFab(value.disp_old_mf, 1, 2, "dispold");
77 0 : value.RegisterGeneralFab(value.vel_old_mf, 1, 2, "velold");
78 0 : value.RegisterGeneralFab(value.ddw_mf, 1, 2);
79 0 : pp_forbid("viscous.mu","replaced with viscous.mu_dashpot");
80 0 : pp_query_default("viscous.mu_dashpot", value.mu_dashpot,0.0); // Dashpot damping (damps velocity)
81 0 : pp_forbid("viscous.mu2","replaced with viscous.mu_newton");
82 0 : pp_query_default("viscous.mu_newton", value.mu_newton,0.0); // Newtonian viscous damping (damps velocity gradient)
83 :
84 0 : std::string velocity_ic_str;
85 0 : pp_query_validate("velocity.ic.type",velocity_ic_str,{"none","expression"}); // Initializer for RHS
86 0 : if (velocity_ic_str == "expression") value.velocity_ic = new IC::Expression(value.geom, pp,"velocity.ic.expression");
87 : }
88 :
89 19 : pp.select<BC::Operator::Elastic::Constant,BC::Operator::Elastic::TensionTest,BC::Operator::Elastic::Expression>("bc",value.bc);
90 :
91 19 : pp_query_default("print_model", value.m_print_model, false); // Print out model variables (if enabled by model)
92 19 : if (value.m_print_model) value.RegisterGeneralFab(value.model_mf, 1, 2, "model", value.m_time_evolving);
93 19 : else value.RegisterGeneralFab(value.model_mf, 1, 2, value.m_time_evolving);
94 :
95 : // Read in IC for RHS
96 19 : std::string rhstype;
97 19 : pp_query_validate("rhs.type", rhstype,{"none","trig","constant"}); // Initializer for RHS
98 19 : if (rhstype == "trig") value.ic_rhs = new IC::Trig(value.geom, pp, "rhs.trig");
99 15 : else if (rhstype == "constant") value.ic_rhs = new IC::Constant(value.geom, pp, "rhs.constant");
100 :
101 : // Timestep interval for elastic solves (default - solve every time)
102 19 : pp_query_default("interval", value.m_interval, 0);
103 :
104 19 : value.RegisterIntegratedVariable(&(value.disp_hi[0].data()[0]), "disp_xhi_x");
105 19 : value.RegisterIntegratedVariable(&(value.disp_hi[0].data()[1]), "disp_xhi_y");
106 19 : value.RegisterIntegratedVariable(&(value.disp_hi[1].data()[0]), "disp_yhi_x");
107 19 : value.RegisterIntegratedVariable(&(value.disp_hi[1].data()[1]), "disp_yhi_y");
108 19 : value.RegisterIntegratedVariable(&(value.trac_hi[0].data()[0]), "trac_xhi_x");
109 19 : value.RegisterIntegratedVariable(&(value.trac_hi[0].data()[1]), "trac_xhi_y");
110 19 : value.RegisterIntegratedVariable(&(value.trac_hi[1].data()[0]), "trac_yhi_x");
111 19 : value.RegisterIntegratedVariable(&(value.trac_hi[1].data()[1]), "trac_yhi_y");
112 :
113 : // Maximum multigrid coarsening level (default - none, maximum coarsening)
114 19 : pp_query_default("max_coarsening_level", value.m_max_coarsening_level,-1);
115 :
116 : // Whether to include residual output field
117 19 : pp_query_default("print_residual", value.m_print_residual,false);
118 19 : if (value.m_print_residual) value.RegisterGeneralFab(value.res_mf, 1, 2, "res", false);
119 :
120 : // Whether to refine based on elastic solution
121 19 : pp_query_default("elastic_ref_threshold", value.m_elastic_ref_threshold,0.01);
122 :
123 : // Set this to true to zero out the displacement before each solve.
124 : // (This is a temporary fix - we need to figure out why this is needed.)
125 19 : pp_query_default("zero_out_displacement", value.m_zero_out_displacement,false);
126 :
127 : // Time to start doing the elastic solve (by default, start immediately)
128 19 : pp_query_default("tstart", value.tstart,-1.0);
129 : }
130 :
131 : protected:
132 : /// \brief Use the #ic object to initialize#Temp
133 121 : void Initialize(int lev) override
134 : {
135 : BL_PROFILE("Integrator::Base::Mechanics::Initialize");
136 121 : if (m_type == Mechanics<MODEL>::Type::Disable) return;
137 :
138 74 : disp_mf[lev]->setVal(Set::Vector::Zero());
139 : //disp_old_mf[lev]->setVal(Set::Vector::Zero());
140 :
141 74 : if (m_type == Type::Dynamic)
142 0 : if (velocity_ic)
143 : {
144 0 : velocity_ic->Initialize(lev,vel_mf);
145 : }
146 :
147 74 : if (ic_rhs) ic_rhs->Initialize(lev, rhs_mf);
148 67 : else rhs_mf[lev]->setVal(Set::Vector::Zero());
149 : }
150 :
151 : virtual void UpdateModel(int a_step, Set::Scalar a_time) = 0;
152 :
153 682 : virtual void TimeStepBegin(Set::Scalar a_time, int a_step) override
154 : {
155 : BL_PROFILE("Integrator::Base::Mechanics::TimeStepBegin");
156 682 : if (m_type == Mechanics<MODEL>::Type::Disable) return;
157 :
158 898 : for (int lev = 0; lev <= finest_level; ++lev)
159 : {
160 472 : if (m_zero_out_displacement) disp_mf[lev]->setVal(Set::Vector::Zero());
161 472 : rhs_mf[lev]->setVal(Set::Vector::Zero());
162 472 : if (ic_rhs) ic_rhs->Initialize(lev, rhs_mf);
163 : }
164 :
165 426 : UpdateModel(a_step, a_time);
166 :
167 426 : bc->SetTime(a_time);
168 426 : bc->Init(rhs_mf, geom);
169 :
170 426 : if (m_type != Mechanics<MODEL>::Type::Static) return;
171 426 : if (a_time < tstart) return;
172 426 : if (m_interval && a_step % m_interval) return;
173 :
174 426 : amrex::LPInfo info;
175 426 : if (m_max_coarsening_level >= 0)
176 0 : info.setMaxCoarseningLevel(m_max_coarsening_level);
177 1278 : Operator::Elastic<MODEL::sym> elastic_op(Geom(0, finest_level), grids, DistributionMap(0, finest_level), info);
178 426 : elastic_op.SetUniform(false);
179 426 : elastic_op.SetHomogeneous(false);
180 426 : elastic_op.SetBC(bc);
181 1278 : IO::ParmParse pp("elasticop");
182 426 : pp_queryclass(elastic_op);
183 :
184 426 : Set::Scalar tol_rel = 1E-8, tol_abs = 1E-8;
185 :
186 426 : solver.Define(elastic_op);
187 426 : if (psi_on) solver.setPsi(psi_mf);
188 426 : solver.solve(disp_mf, rhs_mf, model_mf, tol_rel, tol_abs);
189 426 : if (m_print_residual) solver.compLinearSolverResidual(res_mf, disp_mf, rhs_mf);
190 426 : solver.Clear();
191 :
192 898 : for (int lev = 0; lev <= disp_mf.finest_level; lev++)
193 : {
194 472 : amrex::Box domain = geom[lev].Domain();
195 472 : domain.convert(amrex::IntVect::TheNodeVector());
196 :
197 472 : const amrex::Real* DX = geom[lev].CellSize();
198 1077 : for (MFIter mfi(*disp_mf[lev], false); mfi.isValid(); ++mfi)
199 : {
200 605 : amrex::Box bx = mfi.nodaltilebox();
201 605 : bx.grow(2);
202 605 : bx = bx & domain;
203 605 : amrex::Array4<MODEL> const& model = model_mf[lev]->array(mfi);
204 605 : amrex::Array4<Set::Matrix> const& stress = stress_mf[lev]->array(mfi);
205 605 : amrex::Array4<Set::Matrix> const& strain = strain_mf[lev]->array(mfi);
206 605 : amrex::Array4<const Set::Vector> const& disp = disp_mf[lev]->array(mfi);
207 :
208 :
209 187862 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
210 : {
211 187257 : auto sten = Numeric::GetStencil(i, j, k, bx);
212 374514 : if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
213 : {
214 97068 : Set::Matrix F = Set::Matrix::Identity() + Numeric::Gradient(disp, i, j, k, DX, sten);
215 145602 : stress(i, j, k) = model(i, j, k).DW(F);
216 97068 : strain(i, j, k) = F;
217 : }
218 : else
219 : {
220 138723 : Set::Matrix gradu = Numeric::Gradient(disp, i, j, k, DX, sten);
221 416169 : stress(i, j, k) = model(i, j, k).DW(gradu);
222 277446 : strain(i, j, k) = 0.5 * (gradu + gradu.transpose());
223 : }
224 : });
225 : }
226 472 : Util::RealFillBoundary(*stress_mf[lev], geom[lev]);
227 472 : Util::RealFillBoundary(*strain_mf[lev], geom[lev]);
228 : }
229 : }
230 :
231 55102 : void Advance(int lev, Set::Scalar time, Set::Scalar dt) override
232 : {
233 : BL_PROFILE("Integrator::Base::Mechanics::Advance");
234 55102 : if (m_type == Mechanics<MODEL>::Type::Disable) return;
235 718 : const amrex::Real* DX = geom[lev].CellSize();
236 :
237 718 : amrex::Box domain = geom[lev].Domain();
238 1436 : domain.convert(amrex::IntVect::TheNodeVector());
239 718 : const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
240 :
241 718 : if (m_type == Type::Dynamic)
242 : {
243 :
244 0 : std::swap(*disp_mf[lev], *disp_old_mf[lev]);
245 0 : std::swap(*vel_mf[lev], *vel_old_mf[lev]);
246 :
247 :
248 0 : for (amrex::MFIter mfi(*disp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
249 : {
250 0 : Set::Patch<const Set::Vector> u = disp_old_mf.Patch(lev,mfi);
251 0 : Set::Patch<const Set::Vector> v = vel_old_mf.Patch(lev,mfi);
252 0 : Set::Patch<const Set::Vector> b = rhs_mf.Patch(lev,mfi);
253 :
254 0 : Set::Patch<Set::Vector> unew = disp_mf.Patch(lev,mfi);
255 0 : Set::Patch<Set::Vector> vnew = vel_mf.Patch(lev,mfi);
256 0 : Set::Patch<Set::Matrix> eps = strain_mf.Patch(lev,mfi);
257 0 : Set::Patch<Set::Matrix> sig = stress_mf.Patch(lev,mfi);
258 : //Set::Patch<MATRIX4> ddw = ddw_mf.Patch(lev,mfi);
259 0 : Set::Patch<MODEL> model = model_mf.Patch(lev,mfi);
260 :
261 0 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
262 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
263 : {
264 0 : auto sten = Numeric::GetStencil(i,j,k,bx);
265 0 : eps(i, j, k) = Numeric::Gradient(u, i, j, k, DX,sten);
266 0 : sig(i, j, k) = model(i, j, k).DW(eps(i, j, k));
267 : //ddw(i,j,k) = model(i,j,k).DDW(eps(i,j,k));
268 : });
269 :
270 0 : bx = mfi.nodaltilebox() & domain;
271 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
272 : {
273 :
274 0 : bool AMREX_D_DECL(xmin = (i == lo.x), ymin = (j == lo.y), zmin = (k == lo.z));
275 0 : bool AMREX_D_DECL(xmax = (i == hi.x), ymax = (j == hi.y), zmax = (k == hi.z));
276 :
277 0 : if (AMREX_D_TERM(xmax || xmin, || ymax || ymin, || zmax || zmin))
278 : {
279 0 : auto sten = Numeric::GetStencil(i,j,k,domain);
280 :
281 0 : auto bctype = bc->getType(i,j,k,domain);
282 :
283 0 : for (int d = 0; d < AMREX_SPACEDIM; d++)
284 : {
285 0 : if (bctype[d] == BC::Operator::Elastic::Elastic::Type::Displacement)
286 : {
287 0 : unew(i,j,k)(d) = b(i,j,k)(d);
288 : }
289 0 : else if (bctype[d] == BC::Operator::Elastic::Elastic::Type::Traction)
290 : {
291 :
292 0 : Set::Vector N = Set::Normal(AMREX_D_DECL(xmin,ymin,zmin),
293 : AMREX_D_DECL(xmax,ymax,zmax));
294 :
295 0 : auto [phi, offdiag] = Numeric::GradientSplit(u,i,j,k,DX,sten);
296 :
297 0 : Set::Matrix A = Set::Matrix::Zero();
298 0 : Set::Vector rhs = b(i,j,k);
299 0 : Set::Matrix DW_F0 = model(i,j,k).DW(Set::Matrix::Zero());
300 0 : MATRIX4 ddw = model(i,j,k).DDW(eps(i,j,k));
301 : //Util::Message(INFO,b(i,j,k).transpose());
302 :
303 0 : for (int p = 0; p < AMREX_SPACEDIM; p++)
304 0 : for (int q = 0; q < AMREX_SPACEDIM; q++)
305 : {
306 0 : for (int r = 0; r < AMREX_SPACEDIM; r++)
307 0 : for (int s = 0; s < AMREX_SPACEDIM; s++)
308 : {
309 0 : A(p,r) += ddw(p,q,r,s) * phi(s) * N(q);
310 :
311 0 : rhs(p) -= ddw(p,q,r,s) * eps(i,j,k)(r,s) * N(q);
312 : }
313 :
314 0 : rhs(p) -= DW_F0(p,q)*N(q);
315 : }
316 0 : Set::Vector delta_u = (A.inverse() * rhs);
317 :
318 0 : unew(i,j,k)(d) = u(i,j,k)(d) + delta_u(d);
319 0 : vnew(i,j,k)(d) = (unew(i,j,k)(d) - u(i,j,k)(d))/dt;
320 : }
321 : else
322 : {
323 0 : Util::Abort(INFO,"Elastic dynamics not supported for other BCs yet");
324 : }
325 0 : }
326 : }
327 : else
328 : {
329 : //Set::Matrix gradu = Numeric::Gradient(u,i,j,k,DX);
330 : //Set::Matrix3 gradgradu = Numeric::Hessian(u,i,j,k,DX);
331 :
332 : //Set::Vector f = ddw(i,j,k)*gradgradu;
333 0 : Set::Vector f = Numeric::Divergence(sig, i, j, k, DX);
334 :
335 : //MATRIX4 AMREX_D_DECL(
336 : // Cgrad1 = (Numeric::Stencil<MATRIX4, 1, 0, 0>::D(ddw, i, j, k, 0, DX)),
337 : // Cgrad2 = (Numeric::Stencil<MATRIX4, 0, 1, 0>::D(ddw, i, j, k, 0, DX)),
338 : // Cgrad3 = (Numeric::Stencil<MATRIX4, 0, 0, 1>::D(ddw, i, j, k, 0, DX)));
339 :
340 : //f += AMREX_D_TERM( ( Cgrad1*gradu).col(0),
341 : // +(Cgrad2*gradu).col(1),
342 : // +(Cgrad3*gradu).col(2));
343 :
344 0 : Set::Vector lapv = Numeric::Laplacian(v, i, j, k, DX);
345 0 : f += mu_newton * lapv + b(i,j,k) - mu_dashpot*v(i,j,k);
346 :
347 0 : Set::Vector udotdot = f / rho;
348 0 : vnew(i, j, k) = v(i, j, k) + dt * udotdot;
349 0 : unew(i, j, k) = u(i, j, k) + dt * v(i, j, k);
350 : }
351 : });
352 : }
353 : }
354 :
355 4462 : for (amrex::MFIter mfi(*disp_mf[lev], false); mfi.isValid(); ++mfi)
356 : {
357 3744 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
358 3744 : amrex::Array4<Set::Matrix> const& eps = (*strain_mf[lev]).array(mfi);
359 3744 : amrex::Array4<Set::Matrix> const& sig = (*stress_mf[lev]).array(mfi);
360 3744 : amrex::Array4<MODEL> const& model = (*model_mf[lev]).array(mfi);
361 7012790 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
362 : {
363 7078642 : model(i, j, k).Advance(dt, eps(i, j, k), sig(i, j, k),time);
364 : });
365 : }
366 718 : Util::RealFillBoundary(*model_mf[lev], geom[lev]);
367 :
368 : }
369 :
370 12754 : void Integrate(int amrlev, Set::Scalar /*time*/, int /*step*/,
371 : const amrex::MFIter& mfi, const amrex::Box& a_box) override
372 : {
373 : BL_PROFILE("Integrator::Base::Mechanics::Integrate");
374 12754 : if (m_type == Type::Disable) return;
375 :
376 856 : if (amrex::ParallelDescriptor::NProcs() > 1)
377 : {
378 0 : Util::Warning(INFO,"There is a known bug when calculating trac/disp in Base::Mechanics in parallel.");
379 0 : Util::Warning(INFO,"The thermo.dat values likely will not be correct; use the boxlib output instead.");
380 : }
381 :
382 856 : const amrex::Real* DX = geom[amrlev].CellSize();
383 856 : amrex::Box domain = geom[amrlev].Domain();
384 856 : domain.convert(amrex::IntVect::TheNodeVector());
385 :
386 856 : amrex::Box box = a_box;
387 856 : box.convert(amrex::IntVect::TheNodeVector());
388 :
389 :
390 : //Set::Scalar dv = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
391 : #if AMREX_SPACEDIM == 2
392 556 : Set::Vector da0(DX[1], 0);
393 556 : Set::Vector da1(0, DX[0]);
394 : #elif AMREX_SPACEDIM == 3
395 300 : Set::Vector da(DX[1] * DX[2], 0, 0);
396 : #endif
397 :
398 856 : const Dim3 /*lo= amrex::lbound(domain),*/ hi = amrex::ubound(domain);
399 856 : const Dim3 /*boxlo= amrex::lbound(box),*/ boxhi = amrex::ubound(box);
400 :
401 856 : amrex::Array4<const Set::Matrix> const& stress = (*stress_mf[amrlev]).array(mfi);
402 856 : amrex::Array4<const Set::Vector> const& disp = (*disp_mf[amrlev]).array(mfi);
403 136248 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
404 : {
405 : #if AMREX_SPACEDIM == 2
406 97892 : if (i == hi.x && j < boxhi.y)
407 : {
408 3072 : trac_hi[0] += (0.5 * (stress(i, j, k) + stress(i, j + 1, k)) * da0);
409 2048 : disp_hi[0] = disp(i, j, k);
410 : }
411 97892 : if (j == hi.y && i < boxhi.x)
412 : {
413 3264 : trac_hi[1] += (0.5 * (stress(i, j, k) + stress(i + 1, j, k)) * da1);
414 2176 : disp_hi[1] = disp(i, j, k);
415 : }
416 : #elif AMREX_SPACEDIM == 3
417 37500 : if (i == hi.x && (j < boxhi.y && k < boxhi.z))
418 : {
419 9600 : trac_hi[0] += (0.25 * (stress(i, j, k) + stress(i, j + 1, k)
420 14400 : + stress(i, j, k + 1) + stress(i, j + 1, k + 1)) * da);
421 9600 : disp_hi[0] = disp(i, j, k);
422 : }
423 : #endif
424 : });
425 :
426 : }
427 :
428 528 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/) override
429 : {
430 : BL_PROFILE("Integrator::Base::Mechanics::TagCellsForRefinement");
431 528 : if (m_type == Type::Disable) return;
432 :
433 215 : Set::Vector DX(geom[lev].CellSize());
434 215 : Set::Scalar DXnorm = DX.lpNorm<2>();
435 1209 : for (amrex::MFIter mfi(*strain_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
436 : {
437 994 : amrex::Box bx = mfi.tilebox();
438 994 : bx.convert(amrex::IntVect::TheCellVector());
439 994 : amrex::Array4<char> const& tags = a_tags.array(mfi);
440 994 : amrex::Array4<Set::Matrix> const& eps = strain_mf[lev]->array(mfi);
441 372402 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
442 : {
443 742816 : Set::Matrix3 grad = Numeric::NodeGradientOnCell(eps, i, j, k, DX.data());
444 371408 : if (grad.norm() * DXnorm > m_elastic_ref_threshold)
445 412 : tags(i, j, k) = amrex::TagBox::SET;
446 : });
447 : }
448 : }
449 :
450 : protected:
451 : typedef Set::Matrix4<AMREX_SPACEDIM,MODEL::sym> MATRIX4;
452 : Set::Field<MODEL> model_mf;
453 : Set::Field<MATRIX4> ddw_mf;
454 : Set::Field<Set::Scalar> psi_mf;
455 : bool psi_on = false;
456 :
457 : int m_interval = 0;
458 : Type m_type = Type::Static;
459 :
460 : Set::Field<Set::Vector> disp_mf;
461 : Set::Field<Set::Vector> rhs_mf;
462 : Set::Field<Set::Vector> res_mf;
463 : Set::Field<Set::Matrix> stress_mf;
464 : Set::Field<Set::Matrix> strain_mf;
465 :
466 : // Only use these if using the "dynamics" option
467 : Set::Field<Set::Vector> disp_old_mf;
468 : Set::Field<Set::Vector> vel_mf;
469 : Set::Field<Set::Vector> vel_old_mf;
470 : //Set::Field<Set::Matrix4<AMREX_SPACEDIM,MODEL::sym>> ddw_mf;
471 : Set::Scalar rho = 1.0;
472 : Set::Scalar mu_dashpot = NAN;
473 : Set::Scalar mu_newton = NAN;
474 :
475 : //Set::Vector trac_lo[AMREX_SPACEDIM];
476 : Set::Vector trac_hi[AMREX_SPACEDIM];
477 : Set::Vector disp_hi[AMREX_SPACEDIM];
478 :
479 :
480 : IC::IC* ic_rhs = nullptr;
481 : BC::BC<Set::Scalar>* mybc;
482 :
483 : IC::IC* velocity_ic = nullptr;
484 :
485 : Solver::Nonlocal::Newton<MODEL> solver;//(elastic.op);
486 : BC::Operator::Elastic::Elastic* bc;
487 :
488 : Set::Scalar m_elastic_ref_threshold = 0.01;
489 : bool m_print_model = false;
490 : bool m_print_residual = false;
491 : bool m_time_evolving = false;
492 : int m_max_coarsening_level = -1;
493 :
494 : bool m_zero_out_displacement = false;
495 :
496 : bool plot_disp = true;
497 : bool plot_stress = true;
498 : bool plot_strain = true;
499 : bool plot_psi = true;
500 : bool plot_rhs = true;
501 :
502 :
503 : Set::Scalar tstart = -1;
504 : };
505 : }
506 : }
507 : #endif
|
