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