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 192 : Mechanics() : Integrator()
31 36 : {}
32 :
33 36 : ~Mechanics()
34 : {
35 9 : delete ic_rhs;
36 36 : delete velocity_ic;
37 36 : delete bc;
38 72 : }
39 :
40 : // The mechanics integrator manages the solution of an elastic
41 : // solve using the MLMG solver.
42 33 : static void Parse(Mechanics& value, IO::ParmParse& pp)
43 : {
44 : BL_PROFILE("Integrator::Base::Mechanics::Parse");
45 66 : if (pp.contains("type"))
46 : {
47 12 : 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 48 : pp_query_validate("type", type_str, {"disable","static","dynamic"});
53 12 : 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 12 : }
58 33 : 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 60 : pp_query_default("time_evolving", value.m_time_evolving,false);
64 :
65 60 : pp_query_default("plot_disp", value.plot_disp, true); // Include displacement field in output
66 60 : pp_query_default("plot_rhs", value.plot_rhs, true); // Include right-hand side in output
67 60 : pp_query_default("plot_psi", value.plot_psi, true); // Include :math:`\psi` field in output
68 60 : pp_query_default("plot_stress", value.plot_stress, true); // Include stress in output
69 60 : pp_query_default("plot_strain", value.plot_strain, true); // Include strain in output
70 :
71 60 : value.RegisterGeneralFab(value.disp_mf, 1, 2, value.plot_disp, "disp", value.m_time_evolving);
72 60 : value.RegisterGeneralFab(value.rhs_mf, 1, 2, value.plot_rhs, "rhs", value.m_time_evolving);
73 60 : value.RegisterGeneralFab(value.stress_mf, 1, 2, value.plot_stress, "stress", value.m_time_evolving);
74 60 : value.RegisterGeneralFab(value.strain_mf, 1, 2, value.plot_strain, "strain", value.m_time_evolving);
75 :
76 30 : if (value.m_type == Type::Static)
77 : {
78 : // Read parameters for :ref:`Solver::Nonlocal::Newton` solver
79 90 : pp_queryclass("solver", value.solver);
80 : }
81 30 : 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 60 : pp.select<BC::Operator::Elastic::Constant,BC::Operator::Elastic::TensionTest,BC::Operator::Elastic::Expression>("bc",value.bc);
99 :
100 60 : pp_query_default("print_model", value.m_print_model, false); // Print out model variables (if enabled by model)
101 46 : if (value.m_print_model) value.RegisterGeneralFab(value.model_mf, 1, 2, "model", value.m_time_evolving);
102 22 : 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 60 : if (pp.contains("rhs.type"))
108 : // initial condition for right hand side (body force)
109 27 : 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 60 : pp_query_default("interval", value.m_interval, 0);
113 :
114 60 : value.RegisterIntegratedVariable(&(value.disp_hi[0].data()[0]), "disp_xhi_x");
115 60 : value.RegisterIntegratedVariable(&(value.disp_hi[0].data()[1]), "disp_xhi_y");
116 60 : value.RegisterIntegratedVariable(&(value.disp_hi[1].data()[0]), "disp_yhi_x");
117 60 : value.RegisterIntegratedVariable(&(value.disp_hi[1].data()[1]), "disp_yhi_y");
118 60 : value.RegisterIntegratedVariable(&(value.trac_hi[0].data()[0]), "trac_xhi_x");
119 60 : value.RegisterIntegratedVariable(&(value.trac_hi[0].data()[1]), "trac_xhi_y");
120 60 : value.RegisterIntegratedVariable(&(value.trac_hi[1].data()[0]), "trac_yhi_x");
121 60 : value.RegisterIntegratedVariable(&(value.trac_hi[1].data()[1]), "trac_yhi_y");
122 :
123 : // Maximum multigrid coarsening level (default - none, maximum coarsening)
124 60 : pp_query_default("max_coarsening_level", value.m_max_coarsening_level,-1);
125 :
126 : // Whether to include residual output field
127 60 : pp_query_default("print_residual", value.m_print_residual,false);
128 52 : if (value.m_print_residual) value.RegisterGeneralFab(value.res_mf, 1, 2, "res", false);
129 :
130 : // Whether to refine based on elastic solution
131 60 : 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 60 : 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 60 : pp_query_default("tstart", value.tstart,-1.0);
139 :
140 :
141 : // Relative tolerance in mechanics solve
142 60 : pp.query_default("tol_rel",value.tol_rel,1E-8);
143 :
144 : // Absolute tolerance in mechanics solve
145 90 : pp.query_default("tol_abs",value.tol_abs,1E-8);
146 : }
147 :
148 : protected:
149 : /// \brief Use the #ic object to initialize#Temp
150 130 : void Initialize(int lev) override
151 : {
152 : BL_PROFILE("Integrator::Base::Mechanics::Initialize");
153 130 : if (m_type == Mechanics<MODEL>::Type::Disable) return;
154 :
155 93 : disp_mf[lev]->setVal(Set::Vector::Zero());
156 : //disp_old_mf[lev]->setVal(Set::Vector::Zero());
157 :
158 93 : if (m_type == Type::Dynamic)
159 0 : if (velocity_ic)
160 : {
161 0 : velocity_ic->Initialize(lev,vel_mf);
162 : }
163 :
164 93 : if (ic_rhs) ic_rhs->Initialize(lev, rhs_mf);
165 76 : else rhs_mf[lev]->setVal(Set::Vector::Zero());
166 : }
167 :
168 : virtual void UpdateModel(int a_step, Set::Scalar a_time) = 0;
169 :
170 1596 : virtual void TimeStepBegin(Set::Scalar a_time, int a_step) override
171 : {
172 : BL_PROFILE("Integrator::Base::Mechanics::TimeStepBegin");
173 1596 : if (m_type == Mechanics<MODEL>::Type::Disable) return;
174 :
175 1522 : for (int lev = 0; lev <= finest_level; ++lev)
176 : {
177 788 : rhs_mf[lev]->setVal(Set::Vector::Zero());
178 788 : if (ic_rhs) ic_rhs->Initialize(lev, rhs_mf);
179 : }
180 :
181 734 : UpdateModel(a_step, a_time);
182 :
183 734 : bc->SetTime(a_time);
184 734 : bc->Init(rhs_mf, geom);
185 :
186 734 : if (m_type != Mechanics<MODEL>::Type::Static) return;
187 734 : if (a_time < tstart) return;
188 734 : if (m_interval && a_step % m_interval) return;
189 :
190 734 : amrex::LPInfo info;
191 734 : if (m_max_coarsening_level >= 0)
192 1 : info.setMaxCoarseningLevel(m_max_coarsening_level);
193 734 : Operator::Elastic<MODEL::sym> elastic_op(Geom(0, finest_level), grids, DistributionMap(0, finest_level), info);
194 734 : elastic_op.SetUniform(false);
195 734 : elastic_op.SetHomogeneous(false);
196 734 : elastic_op.SetBC(bc);
197 734 : IO::ParmParse pp("elasticop");
198 734 : pp_queryclass(elastic_op);
199 :
200 734 : solver.Define(elastic_op);
201 734 : if (psi_on) solver.setPsi(psi_mf);
202 :
203 1522 : for (int lev = 0; lev <= finest_level; ++lev)
204 788 : if (m_zero_out_displacement) disp_mf[lev]->setVal(Set::Vector::Zero());
205 :
206 734 : solver.solve(disp_mf, rhs_mf, model_mf, tol_rel, tol_abs);
207 734 : if (m_print_residual) solver.compLinearSolverResidual(res_mf, disp_mf, rhs_mf);
208 734 : solver.Clear();
209 :
210 1522 : for (int lev = 0; lev <= disp_mf.finest_level; lev++)
211 : {
212 788 : amrex::Box domain = geom[lev].Domain();
213 788 : domain.convert(amrex::IntVect::TheNodeVector());
214 :
215 788 : const amrex::Real* DX = geom[lev].CellSize();
216 1709 : for (MFIter mfi(*disp_mf[lev], false); mfi.isValid(); ++mfi)
217 : {
218 921 : amrex::Box bx = mfi.nodaltilebox();
219 921 : bx.grow(2);
220 921 : bx = bx & domain;
221 921 : amrex::Array4<MODEL> const& model = model_mf[lev]->array(mfi);
222 921 : amrex::Array4<Set::Matrix> const& stress = stress_mf[lev]->array(mfi);
223 921 : amrex::Array4<Set::Matrix> const& strain = strain_mf[lev]->array(mfi);
224 921 : amrex::Array4<const Set::Vector> const& disp = disp_mf[lev]->array(mfi);
225 :
226 :
227 253158 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
228 : {
229 252237 : auto sten = Numeric::GetStencil(i, j, k, bx);
230 504474 : if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
231 : {
232 97068 : Set::Matrix F = Set::Matrix::Identity() + Numeric::Gradient(disp, i, j, k, DX, sten);
233 145602 : stress(i, j, k) = model(i, j, k).DW(F);
234 97068 : strain(i, j, k) = F;
235 : }
236 : else
237 : {
238 203703 : Set::Matrix gradu = Numeric::Gradient(disp, i, j, k, DX, sten);
239 611109 : stress(i, j, k) = model(i, j, k).DW(gradu);
240 407406 : strain(i, j, k) = 0.5 * (gradu + gradu.transpose());
241 : }
242 : });
243 : }
244 788 : stress_mf[lev]->setMultiGhost(true);
245 788 : stress_mf[lev]->FillBoundaryAndSync(geom[lev].periodicity());
246 788 : strain_mf[lev]->setMultiGhost(true);
247 788 : strain_mf[lev]->FillBoundaryAndSync(geom[lev].periodicity());
248 : }
249 734 : }
250 :
251 9578 : void Advance(int lev, Set::Scalar time, Set::Scalar dt) override
252 : {
253 : BL_PROFILE("Integrator::Base::Mechanics::Advance");
254 9578 : if (m_type == Mechanics<MODEL>::Type::Disable) return;
255 1048 : const amrex::Real* DX = geom[lev].CellSize();
256 :
257 1048 : amrex::Box domain = geom[lev].Domain();
258 2096 : domain.convert(amrex::IntVect::TheNodeVector());
259 1048 : const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
260 :
261 1048 : if (m_type == Type::Dynamic)
262 : {
263 :
264 0 : std::swap(*disp_mf[lev], *disp_old_mf[lev]);
265 0 : std::swap(*vel_mf[lev], *vel_old_mf[lev]);
266 :
267 :
268 0 : for (amrex::MFIter mfi(*disp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
269 : {
270 0 : Set::Patch<const Set::Vector> u = disp_old_mf.Patch(lev,mfi);
271 0 : Set::Patch<const Set::Vector> v = vel_old_mf.Patch(lev,mfi);
272 0 : Set::Patch<const Set::Vector> b = rhs_mf.Patch(lev,mfi);
273 :
274 0 : Set::Patch<Set::Vector> unew = disp_mf.Patch(lev,mfi);
275 0 : Set::Patch<Set::Vector> vnew = vel_mf.Patch(lev,mfi);
276 0 : Set::Patch<Set::Matrix> eps = strain_mf.Patch(lev,mfi);
277 0 : Set::Patch<Set::Matrix> sig = stress_mf.Patch(lev,mfi);
278 : //Set::Patch<MATRIX4> ddw = ddw_mf.Patch(lev,mfi);
279 0 : Set::Patch<MODEL> model = model_mf.Patch(lev,mfi);
280 :
281 0 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
282 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
283 : {
284 0 : auto sten = Numeric::GetStencil(i,j,k,bx);
285 0 : eps(i, j, k) = Numeric::Gradient(u, i, j, k, DX,sten);
286 0 : sig(i, j, k) = model(i, j, k).DW(eps(i, j, k));
287 : //ddw(i,j,k) = model(i,j,k).DDW(eps(i,j,k));
288 : });
289 :
290 0 : bx = mfi.nodaltilebox() & domain;
291 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
292 : {
293 :
294 0 : bool AMREX_D_DECL(xmin = (i == lo.x), ymin = (j == lo.y), zmin = (k == lo.z));
295 0 : bool AMREX_D_DECL(xmax = (i == hi.x), ymax = (j == hi.y), zmax = (k == hi.z));
296 :
297 0 : if (AMREX_D_TERM(xmax || xmin, || ymax || ymin, || zmax || zmin))
298 : {
299 0 : auto sten = Numeric::GetStencil(i,j,k,domain);
300 :
301 0 : auto bctype = bc->getType(i,j,k,domain);
302 :
303 0 : for (int d = 0; d < AMREX_SPACEDIM; d++)
304 : {
305 0 : if (bctype[d] == BC::Operator::Elastic::Elastic::Type::Displacement)
306 : {
307 0 : unew(i,j,k)(d) = b(i,j,k)(d);
308 : }
309 0 : else if (bctype[d] == BC::Operator::Elastic::Elastic::Type::Traction)
310 : {
311 :
312 0 : Set::Vector N = Set::Normal(AMREX_D_DECL(xmin,ymin,zmin),
313 : AMREX_D_DECL(xmax,ymax,zmax));
314 :
315 0 : auto [phi, offdiag] = Numeric::GradientSplit(u,i,j,k,DX,sten);
316 :
317 0 : Set::Matrix A = Set::Matrix::Zero();
318 0 : Set::Vector rhs = b(i,j,k);
319 0 : Set::Matrix DW_F0 = model(i,j,k).DW(Set::Matrix::Zero());
320 0 : MATRIX4 ddw = model(i,j,k).DDW(eps(i,j,k));
321 : //Util::Message(INFO,b(i,j,k).transpose());
322 :
323 0 : for (int p = 0; p < AMREX_SPACEDIM; p++)
324 0 : for (int q = 0; q < AMREX_SPACEDIM; q++)
325 : {
326 0 : for (int r = 0; r < AMREX_SPACEDIM; r++)
327 0 : for (int s = 0; s < AMREX_SPACEDIM; s++)
328 : {
329 0 : A(p,r) += ddw(p,q,r,s) * phi(s) * N(q);
330 :
331 0 : rhs(p) -= ddw(p,q,r,s) * eps(i,j,k)(r,s) * N(q);
332 : }
333 :
334 0 : rhs(p) -= DW_F0(p,q)*N(q);
335 : }
336 0 : Set::Vector delta_u = (A.inverse() * rhs);
337 :
338 0 : unew(i,j,k)(d) = u(i,j,k)(d) + delta_u(d);
339 0 : vnew(i,j,k)(d) = (unew(i,j,k)(d) - u(i,j,k)(d))/dt;
340 : }
341 : else
342 : {
343 0 : Util::Abort(INFO,"Elastic dynamics not supported for other BCs yet");
344 : }
345 : }
346 0 : }
347 : else
348 : {
349 : //Set::Matrix gradu = Numeric::Gradient(u,i,j,k,DX);
350 : //Set::Matrix3 gradgradu = Numeric::Hessian(u,i,j,k,DX);
351 :
352 : //Set::Vector f = ddw(i,j,k)*gradgradu;
353 0 : Set::Vector f = Numeric::Divergence(sig, i, j, k, DX);
354 :
355 : //MATRIX4 AMREX_D_DECL(
356 : // Cgrad1 = (Numeric::Stencil<MATRIX4, 1, 0, 0>::D(ddw, i, j, k, 0, DX)),
357 : // Cgrad2 = (Numeric::Stencil<MATRIX4, 0, 1, 0>::D(ddw, i, j, k, 0, DX)),
358 : // Cgrad3 = (Numeric::Stencil<MATRIX4, 0, 0, 1>::D(ddw, i, j, k, 0, DX)));
359 :
360 : //f += AMREX_D_TERM( ( Cgrad1*gradu).col(0),
361 : // +(Cgrad2*gradu).col(1),
362 : // +(Cgrad3*gradu).col(2));
363 :
364 0 : Set::Vector lapv = Numeric::Laplacian(v, i, j, k, DX);
365 0 : f += mu_newton * lapv + b(i,j,k) - mu_dashpot*v(i,j,k);
366 :
367 0 : Set::Vector udotdot = f / rho;
368 0 : vnew(i, j, k) = v(i, j, k) + dt * udotdot;
369 0 : unew(i, j, k) = u(i, j, k) + dt * v(i, j, k);
370 : }
371 : });
372 : }
373 : }
374 :
375 5084 : for (amrex::MFIter mfi(*disp_mf[lev], false); mfi.isValid(); ++mfi)
376 : {
377 4036 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
378 4036 : amrex::Array4<Set::Matrix> const& eps = (*strain_mf[lev]).array(mfi);
379 4036 : amrex::Array4<Set::Matrix> const& sig = (*stress_mf[lev]).array(mfi);
380 4036 : amrex::Array4<MODEL> const& model = (*model_mf[lev]).array(mfi);
381 3699108 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
382 : {
383 7440624 : model(i, j, k).Advance(dt, eps(i, j, k), sig(i, j, k),time);
384 : });
385 : }
386 :
387 1048 : model_mf[lev]->setMultiGhost(true);
388 1048 : model_mf[lev]->FillBoundaryAndSync(geom[lev].periodicity());
389 : }
390 :
391 13084 : void Integrate(int amrlev, Set::Scalar /*time*/, int /*step*/,
392 : const amrex::MFIter& mfi, const amrex::Box& a_box) override
393 : {
394 : BL_PROFILE("Integrator::Base::Mechanics::Integrate");
395 13084 : if (m_type == Type::Disable) return;
396 :
397 1186 : if (amrex::ParallelDescriptor::NProcs() > 1 && a_box.contains(amrex::IntVect::TheZeroVector()) &&
398 : amrlev == 0)
399 : {
400 0 : Util::Warning(INFO,"There is a known bug when calculating trac/disp in Base::Mechanics in parallel.");
401 0 : Util::Warning(INFO,"The thermo.dat values likely will not be correct; use the boxlib output instead.");
402 : }
403 :
404 1186 : const amrex::Real* DX = geom[amrlev].CellSize();
405 1186 : amrex::Box domain = geom[amrlev].Domain();
406 1186 : domain.convert(amrex::IntVect::TheNodeVector());
407 :
408 1186 : amrex::Box box = a_box;
409 1186 : box.convert(amrex::IntVect::TheNodeVector());
410 :
411 :
412 : //Set::Scalar dv = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
413 : #if AMREX_SPACEDIM == 2
414 586 : Set::Vector da0(DX[1], 0);
415 586 : Set::Vector da1(0, DX[0]);
416 : #elif AMREX_SPACEDIM == 3
417 600 : Set::Vector da(DX[1] * DX[2], 0, 0);
418 : #endif
419 :
420 1186 : const Dim3 /*lo= amrex::lbound(domain),*/ hi = amrex::ubound(domain);
421 1186 : const Dim3 /*boxlo= amrex::lbound(box),*/ boxhi = amrex::ubound(box);
422 :
423 1186 : amrex::Array4<const Set::Matrix> const& stress = (*stress_mf[amrlev]).array(mfi);
424 1186 : amrex::Array4<const Set::Vector> const& disp = (*disp_mf[amrlev]).array(mfi);
425 195899 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
426 : {
427 : #if AMREX_SPACEDIM == 2
428 119713 : if (i == hi.x && j < boxhi.y)
429 : {
430 3840 : trac_hi[0] += (0.5 * (stress(i, j, k) + stress(i, j + 1, k)) * da0);
431 2560 : disp_hi[0] = disp(i, j, k);
432 : }
433 119713 : if (j == hi.y && i < boxhi.x)
434 : {
435 4032 : trac_hi[1] += (0.5 * (stress(i, j, k) + stress(i + 1, j, k)) * da1);
436 2688 : disp_hi[1] = disp(i, j, k);
437 : }
438 : #elif AMREX_SPACEDIM == 3
439 75000 : if (i == hi.x && (j < boxhi.y && k < boxhi.z))
440 : {
441 19200 : trac_hi[0] += (0.25 * (stress(i, j, k) + stress(i, j + 1, k)
442 28800 : + stress(i, j, k + 1) + stress(i, j + 1, k + 1)) * da);
443 19200 : disp_hi[0] = disp(i, j, k);
444 : }
445 : #endif
446 : });
447 :
448 : }
449 :
450 455 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/) override
451 : {
452 : BL_PROFILE("Integrator::Base::Mechanics::TagCellsForRefinement");
453 455 : if (m_type == Type::Disable) return;
454 :
455 215 : Set::Vector DX(geom[lev].CellSize());
456 215 : Set::Scalar DXnorm = DX.lpNorm<2>();
457 1202 : for (amrex::MFIter mfi(*strain_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
458 : {
459 987 : amrex::Box bx = mfi.tilebox();
460 987 : bx.convert(amrex::IntVect::TheCellVector());
461 987 : amrex::Array4<char> const& tags = a_tags.array(mfi);
462 987 : amrex::Array4<Set::Matrix> const& eps = strain_mf[lev]->array(mfi);
463 371675 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
464 : {
465 741376 : Set::Matrix3 grad = Numeric::NodeGradientOnCell(eps, i, j, k, DX.data());
466 370688 : if (grad.norm() * DXnorm > m_elastic_ref_threshold)
467 412 : tags(i, j, k) = amrex::TagBox::SET;
468 : });
469 : }
470 : }
471 :
472 : protected:
473 : typedef Set::Matrix4<AMREX_SPACEDIM,MODEL::sym> MATRIX4;
474 : Set::Field<MODEL> model_mf;
475 : Set::Field<MATRIX4> ddw_mf;
476 : Set::Field<Set::Scalar> psi_mf;
477 : bool psi_on = false;
478 :
479 : int m_interval = 0;
480 : Type m_type = Type::Static;
481 :
482 : Set::Field<Set::Vector> disp_mf;
483 : Set::Field<Set::Vector> rhs_mf;
484 : Set::Field<Set::Vector> res_mf;
485 : Set::Field<Set::Matrix> stress_mf;
486 : Set::Field<Set::Matrix> strain_mf;
487 :
488 : // Only use these if using the "dynamics" option
489 : Set::Field<Set::Vector> disp_old_mf;
490 : Set::Field<Set::Vector> vel_mf;
491 : Set::Field<Set::Vector> vel_old_mf;
492 : //Set::Field<Set::Matrix4<AMREX_SPACEDIM,MODEL::sym>> ddw_mf;
493 : Set::Scalar rho = 1.0;
494 : Set::Scalar mu_dashpot = NAN;
495 : Set::Scalar mu_newton = NAN;
496 :
497 : //Set::Vector trac_lo[AMREX_SPACEDIM];
498 : Set::Vector trac_hi[AMREX_SPACEDIM];
499 : Set::Vector disp_hi[AMREX_SPACEDIM];
500 :
501 :
502 : IC::IC<Set::Vector>* ic_rhs = nullptr;
503 :
504 : IC::IC<Set::Vector>* velocity_ic = nullptr;
505 :
506 : Solver::Nonlocal::Newton<MODEL> solver;//(elastic.op);
507 : BC::Operator::Elastic::Elastic* bc = nullptr;
508 :
509 : Set::Scalar m_elastic_ref_threshold = 0.01;
510 : bool m_print_model = false;
511 : bool m_print_residual = false;
512 : bool m_time_evolving = false;
513 : int m_max_coarsening_level = -1;
514 :
515 : bool m_zero_out_displacement = false;
516 :
517 : bool plot_disp = true;
518 : bool plot_stress = true;
519 : bool plot_strain = true;
520 : bool plot_psi = true;
521 : bool plot_rhs = true;
522 :
523 :
524 : Set::Scalar tstart = -1;
525 :
526 : Set::Scalar tol_rel=NAN, tol_abs=NAN;
527 : };
528 : }
529 : }
530 : #endif
|
