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 :
199 : // Elastic operator
200 734 : pp.queryclass(elastic_op);
201 :
202 734 : solver.Define(elastic_op);
203 734 : if (psi_on) solver.setPsi(psi_mf);
204 :
205 1522 : for (int lev = 0; lev <= finest_level; ++lev)
206 788 : if (m_zero_out_displacement) disp_mf[lev]->setVal(Set::Vector::Zero());
207 :
208 734 : solver.solve(disp_mf, rhs_mf, model_mf, tol_rel, tol_abs);
209 734 : if (m_print_residual) solver.compLinearSolverResidual(res_mf, disp_mf, rhs_mf);
210 734 : solver.Clear();
211 :
212 1522 : for (int lev = 0; lev <= disp_mf.finest_level; lev++)
213 : {
214 788 : amrex::Box domain = geom[lev].Domain();
215 788 : domain.convert(amrex::IntVect::TheNodeVector());
216 :
217 788 : const amrex::Real* DX = geom[lev].CellSize();
218 1709 : for (MFIter mfi(*disp_mf[lev], false); mfi.isValid(); ++mfi)
219 : {
220 921 : amrex::Box bx = mfi.nodaltilebox();
221 921 : bx.grow(2);
222 921 : bx = bx & domain;
223 921 : amrex::Array4<MODEL> const& model = model_mf[lev]->array(mfi);
224 921 : amrex::Array4<Set::Matrix> const& stress = stress_mf[lev]->array(mfi);
225 921 : amrex::Array4<Set::Matrix> const& strain = strain_mf[lev]->array(mfi);
226 921 : amrex::Array4<const Set::Vector> const& disp = disp_mf[lev]->array(mfi);
227 :
228 :
229 253158 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
230 : {
231 252237 : auto sten = Numeric::GetStencil(i, j, k, bx);
232 504474 : if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
233 : {
234 97068 : Set::Matrix F = Set::Matrix::Identity() + Numeric::Gradient(disp, i, j, k, DX, sten);
235 145602 : stress(i, j, k) = model(i, j, k).DW(F);
236 97068 : strain(i, j, k) = F;
237 : }
238 : else
239 : {
240 203703 : Set::Matrix gradu = Numeric::Gradient(disp, i, j, k, DX, sten);
241 611109 : stress(i, j, k) = model(i, j, k).DW(gradu);
242 407406 : strain(i, j, k) = 0.5 * (gradu + gradu.transpose());
243 : }
244 : });
245 : }
246 788 : stress_mf[lev]->setMultiGhost(true);
247 788 : stress_mf[lev]->FillBoundaryAndSync(geom[lev].periodicity());
248 788 : strain_mf[lev]->setMultiGhost(true);
249 788 : strain_mf[lev]->FillBoundaryAndSync(geom[lev].periodicity());
250 : }
251 734 : }
252 :
253 9578 : void Advance(int lev, Set::Scalar time, Set::Scalar dt) override
254 : {
255 : BL_PROFILE("Integrator::Base::Mechanics::Advance");
256 9578 : if (m_type == Mechanics<MODEL>::Type::Disable) return;
257 1048 : const amrex::Real* DX = geom[lev].CellSize();
258 :
259 1048 : amrex::Box domain = geom[lev].Domain();
260 2096 : domain.convert(amrex::IntVect::TheNodeVector());
261 1048 : const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
262 :
263 1048 : if (m_type == Type::Dynamic)
264 : {
265 :
266 0 : std::swap(*disp_mf[lev], *disp_old_mf[lev]);
267 0 : std::swap(*vel_mf[lev], *vel_old_mf[lev]);
268 :
269 :
270 0 : for (amrex::MFIter mfi(*disp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
271 : {
272 0 : Set::Patch<const Set::Vector> u = disp_old_mf.Patch(lev,mfi);
273 0 : Set::Patch<const Set::Vector> v = vel_old_mf.Patch(lev,mfi);
274 0 : Set::Patch<const Set::Vector> b = rhs_mf.Patch(lev,mfi);
275 :
276 0 : Set::Patch<Set::Vector> unew = disp_mf.Patch(lev,mfi);
277 0 : Set::Patch<Set::Vector> vnew = vel_mf.Patch(lev,mfi);
278 0 : Set::Patch<Set::Matrix> eps = strain_mf.Patch(lev,mfi);
279 0 : Set::Patch<Set::Matrix> sig = stress_mf.Patch(lev,mfi);
280 : //Set::Patch<MATRIX4> ddw = ddw_mf.Patch(lev,mfi);
281 0 : Set::Patch<MODEL> model = model_mf.Patch(lev,mfi);
282 :
283 0 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
284 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
285 : {
286 0 : auto sten = Numeric::GetStencil(i,j,k,bx);
287 0 : eps(i, j, k) = Numeric::Gradient(u, i, j, k, DX,sten);
288 0 : sig(i, j, k) = model(i, j, k).DW(eps(i, j, k));
289 : //ddw(i,j,k) = model(i,j,k).DDW(eps(i,j,k));
290 : });
291 :
292 0 : bx = mfi.nodaltilebox() & domain;
293 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
294 : {
295 :
296 0 : bool AMREX_D_DECL(xmin = (i == lo.x), ymin = (j == lo.y), zmin = (k == lo.z));
297 0 : bool AMREX_D_DECL(xmax = (i == hi.x), ymax = (j == hi.y), zmax = (k == hi.z));
298 :
299 0 : if (AMREX_D_TERM(xmax || xmin, || ymax || ymin, || zmax || zmin))
300 : {
301 0 : auto sten = Numeric::GetStencil(i,j,k,domain);
302 :
303 0 : auto bctype = bc->getType(i,j,k,domain);
304 :
305 0 : for (int d = 0; d < AMREX_SPACEDIM; d++)
306 : {
307 0 : if (bctype[d] == BC::Operator::Elastic::Elastic::Type::Displacement)
308 : {
309 0 : unew(i,j,k)(d) = b(i,j,k)(d);
310 : }
311 0 : else if (bctype[d] == BC::Operator::Elastic::Elastic::Type::Traction)
312 : {
313 :
314 0 : Set::Vector N = Set::Normal(AMREX_D_DECL(xmin,ymin,zmin),
315 : AMREX_D_DECL(xmax,ymax,zmax));
316 :
317 0 : auto [phi, offdiag] = Numeric::GradientSplit(u,i,j,k,DX,sten);
318 :
319 0 : Set::Matrix A = Set::Matrix::Zero();
320 0 : Set::Vector rhs = b(i,j,k);
321 0 : Set::Matrix DW_F0 = model(i,j,k).DW(Set::Matrix::Zero());
322 0 : MATRIX4 ddw = model(i,j,k).DDW(eps(i,j,k));
323 : //Util::Message(INFO,b(i,j,k).transpose());
324 :
325 0 : for (int p = 0; p < AMREX_SPACEDIM; p++)
326 0 : for (int q = 0; q < AMREX_SPACEDIM; q++)
327 : {
328 0 : for (int r = 0; r < AMREX_SPACEDIM; r++)
329 0 : for (int s = 0; s < AMREX_SPACEDIM; s++)
330 : {
331 0 : A(p,r) += ddw(p,q,r,s) * phi(s) * N(q);
332 :
333 0 : rhs(p) -= ddw(p,q,r,s) * eps(i,j,k)(r,s) * N(q);
334 : }
335 :
336 0 : rhs(p) -= DW_F0(p,q)*N(q);
337 : }
338 0 : Set::Vector delta_u = (A.inverse() * rhs);
339 :
340 0 : unew(i,j,k)(d) = u(i,j,k)(d) + delta_u(d);
341 0 : vnew(i,j,k)(d) = (unew(i,j,k)(d) - u(i,j,k)(d))/dt;
342 : }
343 : else
344 : {
345 0 : Util::Abort(INFO,"Elastic dynamics not supported for other BCs yet");
346 : }
347 : }
348 0 : }
349 : else
350 : {
351 : //Set::Matrix gradu = Numeric::Gradient(u,i,j,k,DX);
352 : //Set::Matrix3 gradgradu = Numeric::Hessian(u,i,j,k,DX);
353 :
354 : //Set::Vector f = ddw(i,j,k)*gradgradu;
355 0 : Set::Vector f = Numeric::Divergence(sig, i, j, k, DX);
356 :
357 : //MATRIX4 AMREX_D_DECL(
358 : // Cgrad1 = (Numeric::Stencil<MATRIX4, 1, 0, 0>::D(ddw, i, j, k, 0, DX)),
359 : // Cgrad2 = (Numeric::Stencil<MATRIX4, 0, 1, 0>::D(ddw, i, j, k, 0, DX)),
360 : // Cgrad3 = (Numeric::Stencil<MATRIX4, 0, 0, 1>::D(ddw, i, j, k, 0, DX)));
361 :
362 : //f += AMREX_D_TERM( ( Cgrad1*gradu).col(0),
363 : // +(Cgrad2*gradu).col(1),
364 : // +(Cgrad3*gradu).col(2));
365 :
366 0 : Set::Vector lapv = Numeric::Laplacian(v, i, j, k, DX);
367 0 : f += mu_newton * lapv + b(i,j,k) - mu_dashpot*v(i,j,k);
368 :
369 0 : Set::Vector udotdot = f / rho;
370 0 : vnew(i, j, k) = v(i, j, k) + dt * udotdot;
371 0 : unew(i, j, k) = u(i, j, k) + dt * v(i, j, k);
372 : }
373 : });
374 : }
375 : }
376 :
377 5084 : for (amrex::MFIter mfi(*disp_mf[lev], false); mfi.isValid(); ++mfi)
378 : {
379 4036 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
380 4036 : amrex::Array4<Set::Matrix> const& eps = (*strain_mf[lev]).array(mfi);
381 4036 : amrex::Array4<Set::Matrix> const& sig = (*stress_mf[lev]).array(mfi);
382 4036 : amrex::Array4<MODEL> const& model = (*model_mf[lev]).array(mfi);
383 3699108 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
384 : {
385 7440624 : model(i, j, k).Advance(dt, eps(i, j, k), sig(i, j, k),time);
386 : });
387 : }
388 :
389 1048 : model_mf[lev]->setMultiGhost(true);
390 1048 : model_mf[lev]->FillBoundaryAndSync(geom[lev].periodicity());
391 : }
392 :
393 13304 : void Integrate(int amrlev, Set::Scalar /*time*/, int /*step*/,
394 : const amrex::MFIter& mfi, const amrex::Box& a_box) override
395 : {
396 : BL_PROFILE("Integrator::Base::Mechanics::Integrate");
397 13304 : if (m_type == Type::Disable) return;
398 :
399 1186 : if (amrex::ParallelDescriptor::NProcs() > 1 && a_box.contains(amrex::IntVect::TheZeroVector()) &&
400 : amrlev == 0)
401 : {
402 0 : Util::Warning(INFO,"There is a known bug when calculating trac/disp in Base::Mechanics in parallel.");
403 0 : Util::Warning(INFO,"The thermo.dat values likely will not be correct; use the boxlib output instead.");
404 : }
405 :
406 1186 : const amrex::Real* DX = geom[amrlev].CellSize();
407 1186 : amrex::Box domain = geom[amrlev].Domain();
408 1186 : domain.convert(amrex::IntVect::TheNodeVector());
409 :
410 1186 : amrex::Box box = a_box;
411 1186 : box.convert(amrex::IntVect::TheNodeVector());
412 :
413 :
414 : //Set::Scalar dv = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
415 : #if AMREX_SPACEDIM == 2
416 586 : Set::Vector da0(DX[1], 0);
417 586 : Set::Vector da1(0, DX[0]);
418 : #elif AMREX_SPACEDIM == 3
419 600 : Set::Vector da(DX[1] * DX[2], 0, 0);
420 : #endif
421 :
422 1186 : const Dim3 /*lo= amrex::lbound(domain),*/ hi = amrex::ubound(domain);
423 1186 : const Dim3 /*boxlo= amrex::lbound(box),*/ boxhi = amrex::ubound(box);
424 :
425 1186 : amrex::Array4<const Set::Matrix> const& stress = (*stress_mf[amrlev]).array(mfi);
426 1186 : amrex::Array4<const Set::Vector> const& disp = (*disp_mf[amrlev]).array(mfi);
427 195899 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
428 : {
429 : #if AMREX_SPACEDIM == 2
430 119713 : if (i == hi.x && j < boxhi.y)
431 : {
432 3840 : trac_hi[0] += (0.5 * (stress(i, j, k) + stress(i, j + 1, k)) * da0);
433 2560 : disp_hi[0] = disp(i, j, k);
434 : }
435 119713 : if (j == hi.y && i < boxhi.x)
436 : {
437 4032 : trac_hi[1] += (0.5 * (stress(i, j, k) + stress(i + 1, j, k)) * da1);
438 2688 : disp_hi[1] = disp(i, j, k);
439 : }
440 : #elif AMREX_SPACEDIM == 3
441 75000 : if (i == hi.x && (j < boxhi.y && k < boxhi.z))
442 : {
443 19200 : trac_hi[0] += (0.25 * (stress(i, j, k) + stress(i, j + 1, k)
444 28800 : + stress(i, j, k + 1) + stress(i, j + 1, k + 1)) * da);
445 19200 : disp_hi[0] = disp(i, j, k);
446 : }
447 : #endif
448 : });
449 :
450 : }
451 :
452 455 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/) override
453 : {
454 : BL_PROFILE("Integrator::Base::Mechanics::TagCellsForRefinement");
455 455 : if (m_type == Type::Disable) return;
456 :
457 215 : Set::Vector DX(geom[lev].CellSize());
458 215 : Set::Scalar DXnorm = DX.lpNorm<2>();
459 1202 : for (amrex::MFIter mfi(*strain_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
460 : {
461 987 : amrex::Box bx = mfi.tilebox();
462 987 : bx.convert(amrex::IntVect::TheCellVector());
463 987 : amrex::Array4<char> const& tags = a_tags.array(mfi);
464 987 : amrex::Array4<Set::Matrix> const& eps = strain_mf[lev]->array(mfi);
465 371675 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
466 : {
467 741376 : Set::Matrix3 grad = Numeric::NodeGradientOnCell(eps, i, j, k, DX.data());
468 370688 : if (grad.norm() * DXnorm > m_elastic_ref_threshold)
469 412 : tags(i, j, k) = amrex::TagBox::SET;
470 : });
471 : }
472 : }
473 :
474 : protected:
475 : typedef Set::Matrix4<AMREX_SPACEDIM,MODEL::sym> MATRIX4;
476 : Set::Field<MODEL> model_mf;
477 : Set::Field<MATRIX4> ddw_mf;
478 : Set::Field<Set::Scalar> psi_mf;
479 : bool psi_on = false;
480 :
481 : int m_interval = 0;
482 : Type m_type = Type::Static;
483 :
484 : Set::Field<Set::Vector> disp_mf;
485 : Set::Field<Set::Vector> rhs_mf;
486 : Set::Field<Set::Vector> res_mf;
487 : Set::Field<Set::Matrix> stress_mf;
488 : Set::Field<Set::Matrix> strain_mf;
489 :
490 : // Only use these if using the "dynamics" option
491 : Set::Field<Set::Vector> disp_old_mf;
492 : Set::Field<Set::Vector> vel_mf;
493 : Set::Field<Set::Vector> vel_old_mf;
494 : //Set::Field<Set::Matrix4<AMREX_SPACEDIM,MODEL::sym>> ddw_mf;
495 : Set::Scalar rho = 1.0;
496 : Set::Scalar mu_dashpot = NAN;
497 : Set::Scalar mu_newton = NAN;
498 :
499 : //Set::Vector trac_lo[AMREX_SPACEDIM];
500 : Set::Vector trac_hi[AMREX_SPACEDIM];
501 : Set::Vector disp_hi[AMREX_SPACEDIM];
502 :
503 :
504 : IC::IC<Set::Vector>* ic_rhs = nullptr;
505 :
506 : IC::IC<Set::Vector>* velocity_ic = nullptr;
507 :
508 : Solver::Nonlocal::Newton<MODEL> solver;//(elastic.op);
509 : BC::Operator::Elastic::Elastic* bc = nullptr;
510 :
511 : Set::Scalar m_elastic_ref_threshold = 0.01;
512 : bool m_print_model = false;
513 : bool m_print_residual = false;
514 : bool m_time_evolving = false;
515 : int m_max_coarsening_level = -1;
516 :
517 : bool m_zero_out_displacement = false;
518 :
519 : bool plot_disp = true;
520 : bool plot_stress = true;
521 : bool plot_strain = true;
522 : bool plot_psi = true;
523 : bool plot_rhs = true;
524 :
525 :
526 : Set::Scalar tstart = -1;
527 :
528 : Set::Scalar tol_rel=NAN, tol_abs=NAN;
529 : };
530 : }
531 : }
532 : #endif
|
