Line data Source code
1 : //
2 : // This is a general purpose integrator that focuses on solving
3 : // elasticity/mechanics equations in the absence of other multiphysics
4 : // simulations.
5 : // It is enabled by :code:`alamo.program=mechanics` if used on its own, in
6 : // which case there is no prefix.
7 : // If it is being used by another integrator, see that integrator to determine
8 : // the value of :code:`[prefix]` (often equal to :code:`elastic`).
9 : //
10 : // This integrator inherits from :ref:`Integrator::Base::Mechanics`; see
11 : // documentation for that integrator for additional parameters.
12 : //
13 : // :bdg-primary-line:`Model setup`
14 : // There are two basic tools for setting up a mechanics problem in this
15 : // integrator.
16 : //
17 : // 1. The :code:`eta` field: this is used to mix models of different types.
18 : // Use :code:`nmodels` to specify how many material models to use, and then
19 : // specify each model as :code:`model1`, :code:`model2`, etc.
20 : // The type of moel is specified using the :code:`alamo.program.mechanics.model`
21 : // input.
22 : //
23 : // Once you have done this, you must determine the spatial distribution of each
24 : // model. This is done by choosing an IC for the eta field with :code:`ic.type`.
25 : // The :code:`IC::Expression` is the most general and recommended.
26 : // The models are then mixed linearly, i.e.
27 : //
28 : // .. math::
29 : //
30 : // W_{\textrm{eff}} = \sum_{i=1}^N W_i\,\eta_i(\mathbf{x})
31 : //
32 : // See the :ref:`Eshelby` test for an example of model mixing.
33 : //
34 : // 2. The :code:`psi` field: this is used specifically for cases where a "void" region
35 : // is desired. Its usage is similar to the :code:`eta` case, and is conceptually
36 : // similar in that it scales the model field to near-zero in order to mimic the
37 : // (lack of) mechanical behavior in an empty region.
38 : // It is important to use :code:`psi` here, for reasons that are discussed in detail
39 : // in
40 : // `this paper <https://doi.org/10.1007/s00466-023-02325-8>`_.
41 : // The initialization of :code:`psi` is similar to that for :code:`eta`.
42 : //
43 : // See the :ref:`PlateHole` and :ref:`RubberPlateHole` for canonical exmaples.
44 : // The :ref:`Integrator::Fracture` and :ref:`Integrator::TopOp` integrators are examples
45 : // of integrators that leverage the psi property.
46 : //
47 : // :bdg-primary-line:`Body forces`
48 : // currently have limited support due to the relatively low number of
49 : // times they are needed. See the :ref:`Integrator::Base::Mechanics` documentation for
50 : // detail.
51 : // See the :ref:`TrigTest` test for examples of current body force implementation.
52 : //
53 : // :bdg-primary-line:`Boundary conditions`
54 : // are implemented using the :ref:`BC::Operator::Elastic` classes.
55 : // See the documentation on these classes for more detail.
56 : // See any of the mechanics-based tests for examples of boundary condition application.
57 : //
58 :
59 : #ifndef INTEGRATOR_MECHANICS_H
60 : #define INTEGRATOR_MECHANICS_H
61 : #include <iostream>
62 : #include <fstream>
63 : #include <iomanip>
64 : #include <numeric>
65 :
66 : #include "AMReX.H"
67 : #include "AMReX_ParallelDescriptor.H"
68 : #include "AMReX_ParmParse.H"
69 :
70 : #include "IO/ParmParse.H"
71 : #include "Integrator/Base/Mechanics.H"
72 :
73 :
74 : #include "IC/IC.H"
75 : #include "BC/BC.H"
76 : #include "BC/Operator/Elastic/Constant.H"
77 : #include "BC/Operator/Elastic/TensionTest.H"
78 : #include "BC/Operator/Elastic/Expression.H"
79 :
80 : #include "IC/Ellipse.H"
81 : #include "IC/Voronoi.H"
82 : #include "IC/Constant.H"
83 : #include "IC/Expression.H"
84 : #include "IC/BMP.H"
85 : #include "IC/PNG.H"
86 : #include "IC/PSRead.H"
87 : #include "IC/Constant.H"
88 : #include "BC/Constant.H"
89 : #include "Numeric/Stencil.H"
90 :
91 : #include "Model/Solid/Solid.H"
92 : #include "Solver/Nonlocal/Linear.H"
93 : #include "Solver/Nonlocal/Newton.H"
94 :
95 : #include "Operator/Operator.H"
96 : #include "Base/Mechanics.H"
97 :
98 : namespace Integrator
99 : {
100 : template<class MODEL>
101 : class Mechanics : virtual public Base::Mechanics<MODEL>
102 : {
103 : public:
104 : static constexpr const char *name = "mechanics";
105 :
106 : Mechanics() : Base::Mechanics<MODEL>() { }
107 18 : Mechanics(IO::ParmParse& pp) : Base::Mechanics<MODEL>()
108 : {
109 18 : Parse(*this, pp);
110 18 : }
111 :
112 36 : ~Mechanics()
113 : {
114 3 : delete ic_eta;
115 18 : delete ic_psi;
116 18 : delete ic_trac_normal;
117 18 : delete bc_psi;
118 18 : delete bc_trac_normal;
119 54 : }
120 :
121 : // Mechanics inputs. See also :ref:`Integrator::Base::Mechanics`
122 18 : static void Parse(Mechanics& value, IO::ParmParse& pp)
123 : {
124 72 : pp.queryclass<Base::Mechanics<MODEL>>(value);
125 :
126 : int nmodels;
127 108 : pp_query_default("nmodels", nmodels, 1); // Number of elastic model varieties
128 :
129 90 : pp.queryclass_enumerate<MODEL>("model",value.models,nmodels);
130 :
131 126 : Util::Assert(INFO, TEST(value.models.size() > 0));
132 54 : value.RegisterNodalFab(value.eta_mf, value.models.size(), 2, "eta", true);
133 : // Refinement threshold for eta field
134 108 : pp_query_default("eta_ref_threshold", value.m_eta_ref_threshold, 0.01);
135 : // Refinement threshold for strain gradient
136 108 : pp_query_default("ref_threshold", value.m_elastic_ref_threshold, 0.01);
137 : // Explicity impose neumann condition on model at domain boundaries (2d only)
138 90 : pp_query_default("model_neumann_boundary",value.model_neumann_boundary,false);
139 :
140 : // Read in IC for eta
141 24 : if (nmodels > 1 && pp.contains("ic.type"))
142 : {
143 : // Select the initial condition for eta
144 6 : pp.select<IC::Constant,IC::Ellipse,IC::Voronoi,IC::BMP,IC::PNG,IC::Expression,IC::PSRead>("ic",value.ic_eta,value.geom);
145 :
146 3 : value.eta_reset_on_regrid = true;
147 : // Whether to re-initialize eta when re-gridding occurs.
148 : // Default is false unless eta ic is set, then default is.
149 : // true.
150 18 : pp_query_default("eta.reset_on_regrid", value.eta_reset_on_regrid, true);
151 : }
152 :
153 : // Read in IC for psi
154 36 : if (pp.contains("psi.ic.type"))
155 : {
156 : // Select initial condition for psi field
157 6 : pp.select<IC::Ellipse,IC::Constant,IC::Expression,IC::PSRead,IC::PNG>("psi.ic",value.ic_psi,value.geom);
158 :
159 3 : value.bc_psi = new BC::Nothing();
160 9 : value.RegisterNewFab(value.psi_mf, value.bc_psi, 1, 2, "psi", value.plot_psi);
161 3 : value.psi_on = true;
162 :
163 3 : value.psi_reset_on_regrid = true;
164 : // Whether to re-initialize psi when re-gridding occurs.
165 : // Default is false unless a psi ic is set, then default is
166 : // true.
167 18 : pp_query_default("psi.reset_on_regrid", value.psi_reset_on_regrid, true);
168 : }
169 :
170 36 : if (pp.contains("trac_normal.ic.type"))
171 : {
172 : // Read in IC for the "normal traction" field (applied at diffuse boundaries)
173 0 : pp.select<IC::Ellipse,IC::Constant,IC::Expression,IC::PSRead>("trac_normal.ic",value.ic_trac_normal,value.geom);
174 :
175 0 : if (value.ic_trac_normal)
176 : {
177 0 : value.bc_trac_normal = new BC::Nothing();
178 0 : value.RegisterNewFab(value.trac_normal_mf, value.bc_trac_normal, 1, 2, "trac_normal", true);
179 : }
180 : }
181 54 : Util::Message(INFO);
182 :
183 18 : }
184 :
185 73 : void Initialize(int lev) override
186 : {
187 73 : Base::Mechanics<MODEL>::Initialize(lev);
188 73 : eta_mf[lev]->setVal(0.0);
189 73 : if (models.size() > 1 && ic_eta) ic_eta->Initialize(lev, eta_mf);
190 45 : else eta_mf[lev]->setVal(1.0);
191 :
192 73 : if (psi_on) ic_psi->Initialize(lev, psi_mf);
193 73 : if (ic_trac_normal) ic_trac_normal->Initialize(lev, trac_normal_mf);
194 73 : }
195 :
196 424 : virtual void UpdateModel(int a_step, Set::Scalar time) override
197 : {
198 424 : if (m_type == Base::Mechanics<MODEL>::Type::Disable) return;
199 :
200 424 : if (ic_trac_normal)
201 : {
202 0 : for (int lev = 0; lev <= finest_level; ++lev)
203 : {
204 0 : ic_trac_normal->Initialize(lev, trac_normal_mf, time);
205 0 : psi_mf[lev]->FillBoundary();
206 0 : amrex::Box domain = this->geom[lev].Domain();
207 0 : domain.convert(amrex::IntVect::TheNodeVector());
208 0 : Set::Vector DX(geom[lev].CellSize());
209 :
210 0 : for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
211 : {
212 0 : amrex::Box bx = mfi.nodaltilebox();
213 0 : bx = bx & domain;
214 :
215 0 : amrex::Array4<Set::Vector> const& rhs = rhs_mf[lev]->array(mfi);
216 0 : amrex::Array4<const Set::Scalar> const& psi = psi_mf[lev]->array(mfi);
217 0 : amrex::Array4<const Set::Scalar> const& trac_normal = trac_normal_mf[lev]->array(mfi);
218 :
219 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
220 0 : Set::Vector grad = Numeric::CellGradientOnNode(psi, i, j, k, 0, DX.data());
221 0 : rhs(i,j,k) = trac_normal(i,j,k) * grad;
222 : });
223 : }
224 0 : Util::RealFillBoundary(*rhs_mf[lev], geom[lev]);
225 : }
226 : }
227 :
228 424 : if (a_step > 0) return;
229 :
230 64 : for (int lev = 0; lev <= finest_level; ++lev)
231 : {
232 46 : eta_mf[lev]->FillBoundary();
233 :
234 46 : amrex::Box domain = this->geom[lev].Domain();
235 46 : domain.convert(amrex::IntVect::TheNodeVector());
236 :
237 46 : Set::Vector DX(geom[lev].CellSize());
238 :
239 225 : for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
240 : {
241 179 : amrex::Box bx = mfi.grownnodaltilebox();
242 179 : bx = bx & domain;
243 :
244 179 : amrex::Array4<MODEL> const& model = model_mf[lev]->array(mfi);
245 179 : amrex::Array4<const Set::Scalar> const& eta = eta_mf[lev]->array(mfi);
246 :
247 105574 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
248 210790 : model(i, j, k) = MODEL::Zero();
249 253150 : for (unsigned int n = 0; n < models.size(); n++)
250 726235 : model(i, j, k) += eta(i, j, k, n) * models[n];
251 : });
252 : }
253 :
254 :
255 46 : if (model_neumann_boundary)
256 : {
257 0 : Util::AssertException(INFO,TEST(AMREX_SPACEDIM==2),"neumann boundary works in 2d only");
258 0 : for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
259 : {
260 0 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
261 0 : amrex::Array4<MODEL> const& model = this->model_mf[lev]->array(mfi);
262 0 : const Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
263 :
264 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
265 : {
266 0 : if (i==lo.x && j==lo.y)
267 0 : model(i,j,k) = 0.5*(model(i+1,j,k)+model(i,j+1,k));
268 0 : else if (i==lo.x && j==hi.y)
269 0 : model(i,j,k) = 0.5*(model(i+1,j,k)+model(i,j-1,k));
270 0 : else if (i==hi.x && j==lo.y)
271 0 : model(i,j,k) = 0.5*(model(i-1,j,k)+model(i,j+1,k));
272 0 : else if (i==hi.x && j==hi.y)
273 0 : model(i,j,k) = 0.5*(model(i-1,j,k)+model(i,j-1,k));
274 :
275 0 : else if (i==lo.x)
276 0 : model(i,j,k) = model(i+1,j,k);
277 0 : else if (i==hi.x)
278 0 : model(i,j,k) = model(i-1,j,k);
279 0 : else if (j==lo.y)
280 0 : model(i,j,k) = model(i,j+1,k);
281 0 : else if (j==hi.y)
282 0 : model(i,j,k) = model(i,j-1,k);
283 : });
284 : }
285 : }
286 :
287 46 : Util::RealFillBoundary(*model_mf[lev], geom[lev]);
288 : }
289 : }
290 :
291 211 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar a_time, int a_ngrow) override
292 : {
293 211 : if (m_type == Base::Mechanics<MODEL>::Type::Disable) return;
294 211 : Base::Mechanics<MODEL>::TagCellsForRefinement(lev, a_tags, a_time, a_ngrow);
295 :
296 211 : Set::Vector DX(geom[lev].CellSize());
297 211 : Set::Scalar DXnorm = DX.lpNorm<2>();
298 1192 : for (amrex::MFIter mfi(*eta_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
299 : {
300 981 : amrex::Box bx = mfi.nodaltilebox();
301 981 : amrex::Array4<char> const& tags = a_tags.array(mfi);
302 981 : amrex::Array4<Set::Scalar> const& eta = eta_mf[lev]->array(mfi);
303 402254 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
304 : {
305 401273 : auto sten = Numeric::GetStencil(i, j, k, bx);
306 1087610 : for (int n = 0; n < eta_mf[lev]->nComp(); n++)
307 : {
308 686337 : Set::Vector grad = Numeric::Gradient(eta, i, j, k, n, DX.data(), sten);
309 686337 : if (grad.lpNorm<2>() * DXnorm > m_eta_ref_threshold)
310 225368 : tags(i, j, k) = amrex::TagBox::SET;
311 : }
312 : });
313 981 : if (psi_on)
314 : {
315 296 : amrex::Array4<Set::Scalar> const& psi = psi_mf[lev]->array(mfi);
316 114712 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
317 : {
318 114416 : auto sten = Numeric::GetStencil(i, j, k, bx);
319 : {
320 114416 : Set::Vector gradpsi = Numeric::Gradient(psi, i, j, k, 0, DX.data(), sten);
321 114416 : if (gradpsi.lpNorm<2>() * DXnorm > m_eta_ref_threshold)
322 42996 : tags(i, j, k) = amrex::TagBox::SET;
323 : }
324 : });
325 : }
326 : }
327 : }
328 :
329 4 : void Regrid(int lev, Set::Scalar time) override
330 : {
331 4 : if (eta_reset_on_regrid && models.size() > 1 && ic_eta) ic_eta->Initialize(lev, eta_mf, time);
332 4 : if (psi_reset_on_regrid) ic_psi->Initialize(lev, psi_mf, time);
333 4 : }
334 :
335 : protected:
336 : Set::Field<Set::Scalar> eta_mf;
337 : Set::Scalar m_eta_ref_threshold = NAN;
338 : std::vector<MODEL> models;
339 : IC::IC<Set::Scalar>* ic_eta = nullptr;
340 : IC::IC<Set::Scalar>* ic_psi = nullptr;
341 : IC::IC<Set::Scalar>* ic_trac_normal = nullptr;
342 : BC::BC<Set::Scalar>* bc_psi = nullptr;
343 : BC::BC<Set::Scalar>* bc_trac_normal = nullptr;
344 : bool psi_reset_on_regrid = false;
345 : bool eta_reset_on_regrid = false;
346 :
347 : bool model_neumann_boundary = false;
348 :
349 : Set::Field<Set::Scalar> trac_normal_mf;
350 :
351 : using Base::Mechanics<MODEL>::m_type;
352 : using Base::Mechanics<MODEL>::finest_level;
353 : using Base::Mechanics<MODEL>::geom;
354 : using Base::Mechanics<MODEL>::model_mf;
355 : using Base::Mechanics<MODEL>::psi_mf;
356 : using Base::Mechanics<MODEL>::psi_on;
357 : using Base::Mechanics<MODEL>::rhs_mf;
358 : };
359 :
360 :
361 :
362 :
363 :
364 :
365 :
366 :
367 :
368 :
369 : } // namespace Integrator
370 : #endif
|
