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