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 29 : Mechanics(IO::ParmParse& pp) : Base::Mechanics<MODEL>()
108 : {
109 29 : Parse(*this, pp);
110 29 : }
111 :
112 58 : ~Mechanics()
113 : {
114 3 : delete ic_eta;
115 29 : delete ic_psi;
116 29 : delete ic_trac_normal;
117 29 : delete bc_psi;
118 29 : delete bc_trac_normal;
119 87 : }
120 :
121 : // Mechanics inputs. See also :ref:`Integrator::Base::Mechanics`
122 29 : static void Parse(Mechanics& value, IO::ParmParse& pp)
123 : {
124 29 : pp.queryclass<Base::Mechanics<MODEL>>(value);
125 :
126 : int nmodels;
127 29 : pp_query_default("nmodels", nmodels, 1); // Number of elastic model varieties
128 :
129 58 : pp.queryclass_enumerate<MODEL>("model",value.models,nmodels);
130 :
131 203 : Util::Assert(INFO, TEST((int)nmodels == (int)value.models.size()),
132 : "nmodels does not match the # of specified models");
133 :
134 203 : Util::Assert(INFO, TEST(value.models.size() > 0));
135 87 : value.RegisterNodalFab(value.eta_mf, value.models.size(), 2, "eta", true);
136 : // Refinement threshold for eta field
137 58 : pp_query_default("eta_ref_threshold", value.m_eta_ref_threshold, 0.01);
138 : // Refinement threshold for strain gradient
139 58 : pp_query_default("ref_threshold", value.m_elastic_ref_threshold, 0.01);
140 : // Explicity impose neumann condition on model at domain boundaries (2d only)
141 58 : pp_query_default("model_neumann_boundary",value.model_neumann_boundary,false);
142 :
143 : // Read in IC for eta
144 35 : if (nmodels > 1 && pp.contains("ic.type"))
145 : {
146 : // Select the initial condition for eta
147 6 : pp.select<IC::Constant,IC::Ellipse,IC::Voronoi,IC::BMP,IC::PNG,IC::Expression,IC::PSRead>("ic",value.ic_eta,value.geom);
148 :
149 3 : value.eta_reset_on_regrid = true;
150 : // Whether to re-initialize eta when re-gridding occurs.
151 : // Default is false unless eta ic is set, then default is.
152 : // true.
153 9 : pp_query_default("eta.reset_on_regrid", value.eta_reset_on_regrid, true);
154 : }
155 :
156 : // Read in IC for psi
157 58 : if (pp.contains("psi.ic.type"))
158 : {
159 : // Select initial condition for psi field
160 6 : pp.select<IC::Ellipse,IC::Constant,IC::Expression,IC::PSRead,IC::PNG>("psi.ic",value.ic_psi,value.geom);
161 :
162 3 : value.bc_psi = new BC::Nothing();
163 9 : value.RegisterNewFab(value.psi_mf, value.bc_psi, 1, 2, "psi", value.plot_psi);
164 3 : value.psi_on = true;
165 :
166 3 : value.psi_reset_on_regrid = true;
167 : // Whether to re-initialize psi when re-gridding occurs.
168 : // Default is false unless a psi ic is set, then default is
169 : // true.
170 9 : pp_query_default("psi.reset_on_regrid", value.psi_reset_on_regrid, true);
171 : }
172 :
173 58 : if (pp.contains("trac_normal.ic.type"))
174 : {
175 : // Read in IC for the "normal traction" field (applied at diffuse boundaries)
176 0 : pp.select<IC::Ellipse,IC::Constant,IC::Expression,IC::PSRead>("trac_normal.ic",value.ic_trac_normal,value.geom);
177 :
178 0 : if (value.ic_trac_normal)
179 : {
180 0 : value.bc_trac_normal = new BC::Nothing();
181 0 : value.RegisterNewFab(value.trac_normal_mf, value.bc_trac_normal, 1, 2, "trac_normal", true);
182 : }
183 : }
184 87 : Util::Message(INFO);
185 :
186 29 : }
187 :
188 92 : void Initialize(int lev) override
189 : {
190 92 : Base::Mechanics<MODEL>::Initialize(lev);
191 92 : eta_mf[lev]->setVal(0.0);
192 92 : if (models.size() > 1 && ic_eta) ic_eta->Initialize(lev, eta_mf);
193 64 : else eta_mf[lev]->setVal(1.0);
194 :
195 92 : if (psi_on) ic_psi->Initialize(lev, psi_mf);
196 92 : if (ic_trac_normal) ic_trac_normal->Initialize(lev, trac_normal_mf);
197 92 : }
198 :
199 732 : virtual void UpdateModel(int a_step, Set::Scalar time) override
200 : {
201 732 : if (m_type == Base::Mechanics<MODEL>::Type::Disable) return;
202 :
203 732 : if (ic_trac_normal)
204 : {
205 0 : for (int lev = 0; lev <= finest_level; ++lev)
206 : {
207 0 : ic_trac_normal->Initialize(lev, trac_normal_mf, time);
208 0 : psi_mf[lev]->FillBoundary();
209 0 : amrex::Box domain = this->geom[lev].Domain();
210 0 : domain.convert(amrex::IntVect::TheNodeVector());
211 0 : Set::Vector DX(geom[lev].CellSize());
212 :
213 0 : for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
214 : {
215 0 : amrex::Box bx = mfi.nodaltilebox();
216 0 : bx = bx & domain;
217 :
218 0 : amrex::Array4<Set::Vector> const& rhs = rhs_mf[lev]->array(mfi);
219 0 : amrex::Array4<const Set::Scalar> const& psi = psi_mf[lev]->array(mfi);
220 0 : amrex::Array4<const Set::Scalar> const& trac_normal = trac_normal_mf[lev]->array(mfi);
221 :
222 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
223 0 : Set::Vector grad = Numeric::CellGradientOnNode(psi, i, j, k, 0, DX.data());
224 0 : rhs(i,j,k) = trac_normal(i,j,k) * grad;
225 : });
226 : }
227 0 : Util::RealFillBoundary(*rhs_mf[lev], geom[lev]);
228 : }
229 : }
230 :
231 732 : if (a_step > 0) return;
232 :
233 94 : for (int lev = 0; lev <= finest_level; ++lev)
234 : {
235 65 : eta_mf[lev]->FillBoundary();
236 :
237 65 : amrex::Box domain = this->geom[lev].Domain();
238 65 : domain.convert(amrex::IntVect::TheNodeVector());
239 :
240 65 : Set::Vector DX(geom[lev].CellSize());
241 :
242 263 : for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
243 : {
244 198 : amrex::Box bx = mfi.grownnodaltilebox();
245 :
246 198 : amrex::Array4<MODEL> const& model = model_mf[lev]->array(mfi);
247 198 : amrex::Array4<const Set::Scalar> const& eta = eta_mf[lev]->array(mfi);
248 :
249 246306 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
250 185130 : model(i, j, k) = MODEL::Zero();
251 330564 : for (unsigned int n = 0; n < models.size(); n++)
252 921038 : model(i, j, k) += eta(i, j, k, n) * models[n];
253 : });
254 : }
255 :
256 :
257 65 : if (model_neumann_boundary)
258 : {
259 0 : Util::AssertException(INFO,TEST(AMREX_SPACEDIM==2),"neumann boundary works in 2d only");
260 0 : for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
261 : {
262 0 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
263 0 : amrex::Array4<MODEL> const& model = this->model_mf[lev]->array(mfi);
264 0 : const Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
265 :
266 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
267 : {
268 0 : if (i==lo.x && j==lo.y)
269 0 : model(i,j,k) = 0.5*(model(i+1,j,k)+model(i,j+1,k));
270 0 : else if (i==lo.x && j==hi.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==lo.y)
273 0 : model(i,j,k) = 0.5*(model(i-1,j,k)+model(i,j+1,k));
274 0 : else if (i==hi.x && j==hi.y)
275 0 : model(i,j,k) = 0.5*(model(i-1,j,k)+model(i,j-1,k));
276 :
277 0 : else if (i==lo.x)
278 0 : model(i,j,k) = model(i+1,j,k);
279 0 : else if (i==hi.x)
280 0 : model(i,j,k) = model(i-1,j,k);
281 0 : else if (j==lo.y)
282 0 : model(i,j,k) = model(i,j+1,k);
283 0 : else if (j==hi.y)
284 0 : model(i,j,k) = model(i,j-1,k);
285 : });
286 : }
287 : }
288 :
289 65 : Util::RealFillBoundary(*model_mf[lev], geom[lev]);
290 : }
291 : }
292 :
293 211 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar a_time, int a_ngrow) override
294 : {
295 211 : if (m_type == Base::Mechanics<MODEL>::Type::Disable) return;
296 211 : Base::Mechanics<MODEL>::TagCellsForRefinement(lev, a_tags, a_time, a_ngrow);
297 :
298 211 : Set::Vector DX(geom[lev].CellSize());
299 211 : Set::Scalar DXnorm = DX.lpNorm<2>();
300 1192 : for (amrex::MFIter mfi(*eta_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
301 : {
302 981 : amrex::Box bx = mfi.nodaltilebox();
303 981 : amrex::Array4<char> const& tags = a_tags.array(mfi);
304 981 : amrex::Array4<Set::Scalar> const& eta = eta_mf[lev]->array(mfi);
305 402254 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
306 : {
307 401273 : auto sten = Numeric::GetStencil(i, j, k, bx);
308 1087610 : for (int n = 0; n < eta_mf[lev]->nComp(); n++)
309 : {
310 686337 : Set::Vector grad = Numeric::Gradient(eta, i, j, k, n, DX.data(), sten);
311 686337 : if (grad.lpNorm<2>() * DXnorm > m_eta_ref_threshold)
312 225368 : tags(i, j, k) = amrex::TagBox::SET;
313 : }
314 : });
315 981 : if (psi_on)
316 : {
317 296 : amrex::Array4<Set::Scalar> const& psi = psi_mf[lev]->array(mfi);
318 114712 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
319 : {
320 114416 : auto sten = Numeric::GetStencil(i, j, k, bx);
321 : {
322 114416 : Set::Vector gradpsi = Numeric::Gradient(psi, i, j, k, 0, DX.data(), sten);
323 114416 : if (gradpsi.lpNorm<2>() * DXnorm > m_eta_ref_threshold)
324 42996 : tags(i, j, k) = amrex::TagBox::SET;
325 : }
326 : });
327 : }
328 : }
329 : }
330 :
331 4 : void Regrid(int lev, Set::Scalar time) override
332 : {
333 4 : if (eta_reset_on_regrid && models.size() > 1 && ic_eta) ic_eta->Initialize(lev, eta_mf, time);
334 4 : if (psi_reset_on_regrid) ic_psi->Initialize(lev, psi_mf, time);
335 4 : }
336 :
337 : protected:
338 : Set::Field<Set::Scalar> eta_mf;
339 : Set::Scalar m_eta_ref_threshold = NAN;
340 : std::vector<MODEL> models;
341 : IC::IC<Set::Scalar>* ic_eta = nullptr;
342 : IC::IC<Set::Scalar>* ic_psi = nullptr;
343 : IC::IC<Set::Scalar>* ic_trac_normal = nullptr;
344 : BC::BC<Set::Scalar>* bc_psi = nullptr;
345 : BC::BC<Set::Scalar>* bc_trac_normal = nullptr;
346 : bool psi_reset_on_regrid = false;
347 : bool eta_reset_on_regrid = false;
348 :
349 : bool model_neumann_boundary = false;
350 :
351 : Set::Field<Set::Scalar> trac_normal_mf;
352 :
353 : using Base::Mechanics<MODEL>::m_type;
354 : using Base::Mechanics<MODEL>::finest_level;
355 : using Base::Mechanics<MODEL>::geom;
356 : using Base::Mechanics<MODEL>::model_mf;
357 : using Base::Mechanics<MODEL>::psi_mf;
358 : using Base::Mechanics<MODEL>::psi_on;
359 : using Base::Mechanics<MODEL>::rhs_mf;
360 : };
361 :
362 :
363 :
364 :
365 :
366 :
367 :
368 :
369 :
370 :
371 : } // namespace Integrator
372 : #endif
|
