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