Alamo
Newton.H
Go to the documentation of this file.
1#ifndef SOLVER_NONLOCAL_NEWTON
2#define SOLVER_NONLOCAL_NEWTON
3
4#include <cmath>
5#include "Set/Set.H"
6#include "Operator/Elastic.H"
7#include "Solver/Nonlocal/Linear.H"
8#include "IO/ParmParse.H"
9#include "Numeric/Stencil.H"
10
11namespace Solver
12{
13namespace Nonlocal
14{
15namespace Detail
16{
17template <typename T>
18inline T
19FaceModel(const amrex::Array4<const T>& model,
20 const int i, const int j, const int k, const int face)
21{
22 const int ip = i + (face == 0);
23 const int jp = j + (face == 1);
24 const int kp = k + (face == 2);
25 // Use the same arithmetic face model for the residual and tangent. Tensor
26 // models do not have a component-independent harmonic average.
27 return 0.5 * (model(i, j, k) + model(ip, jp, kp));
28}
29
30template <typename T, typename U>
31inline Set::Matrix
32FaceKinematicVariable(const T& model,
33 const U& u,
34 const int i, const int j, const int k, const int face,
35 const Set::Scalar dx[AMREX_SPACEDIM])
36{
37 Set::Matrix gradu = Numeric::FaceGradient(u, i, j, k, face, dx);
38 if (model.kinvar == Model::Solid::KinematicVariable::gradu)
39 return gradu;
40 if (model.kinvar == Model::Solid::KinematicVariable::epsilon)
41 return 0.5 * (gradu + gradu.transpose());
42 if (model.kinvar == Model::Solid::KinematicVariable::F)
43 return gradu + Set::Matrix::Identity();
44 return Set::Matrix::Zero();
45}
46
47template <typename T, typename U>
48inline Set::Matrix
49FaceStress(const amrex::Array4<const T>& model,
50 const U& u,
51 const int i, const int j, const int k, const int face,
52 const Set::Scalar dx[AMREX_SPACEDIM])
53{
54 const T face_model = FaceModel(model, i, j, k, face);
55 return face_model.DW(
56 FaceKinematicVariable(face_model, u, i, j, k, face, dx));
57}
58
59template <typename T, typename U>
60inline Set::Matrix4<AMREX_SPACEDIM, T::sym>
61FaceTangent(const amrex::Array4<const T>& model,
62 const U& u,
63 const int i, const int j, const int k, const int face,
64 const Set::Scalar dx[AMREX_SPACEDIM])
65{
66 const T face_model = FaceModel(model, i, j, k, face);
67 return face_model.DDW(
68 FaceKinematicVariable(face_model, u, i, j, k, face, dx));
69}
70
71template <typename T, typename U>
72inline Set::Vector
73FaceStressDivergence(const amrex::Array4<const T>& model,
74 const U& u,
75 const int i, const int j, const int k,
76 const Set::Scalar dx[AMREX_SPACEDIM])
77{
78 Set::Vector ret = Set::Vector::Zero();
79 for (int face = 0; face < AMREX_SPACEDIM; ++face)
80 {
81 const int im = i - (face == 0);
82 const int jm = j - (face == 1);
83 const int km = k - (face == 2);
84 const Set::Matrix hi = FaceStress(model, u, i, j, k, face, dx);
85 const Set::Matrix lo = FaceStress(model, u, im, jm, km, face, dx);
86 ret += (hi.col(face) - lo.col(face)) / dx[face];
87 }
88 return ret;
89}
90} // namespace Detail
91
92template <typename T>
93class Newton: public Linear
94{
95public:
96 Newton() {};
97
98 Newton(Operator::Elastic<T::sym>& a_op)
99 {
100 this->Define(a_op);
101 };
102
103 ~Newton()
104 {
105 if (m_defined) Clear();
106 }
107
108 void Define(Operator::Elastic<T::sym>& a_op)
109 {
110 Linear::Define(a_op);
111 m_elastic = dynamic_cast<Operator::Elastic<T::sym> *>(linop);
112 //m_bc = &m_elastic->GetBC();
113 }
114 void Clear()
115 {
116 Linear::Clear();
117 m_elastic = nullptr;
118 //m_bc = nullptr;
119 }
120
121
122 void setNRIters(int a_nriters) { m_nriters = a_nriters; }
123
124 void setPsi(Set::Field<Set::Scalar>& a_psi)
125 {
126 m_psi = &a_psi;
127 if (m_conservative_face_flux && !m_reported_psi_face_flux_conflict)
128 {
129 Util::Warning(INFO, "solver.conservative_face_flux is disabled because a psi field is active");
130 m_reported_psi_face_flux_conflict = true;
131 }
132 }
133
134 [[nodiscard]] bool usesConservativeFaceFlux() const
135 {
136 return m_conservative_face_flux && m_psi == nullptr;
137 }
138
139private:
140 void prepareForSolve(const Set::Field<Set::Scalar>& a_u_mf,
141 const Set::Field<Set::Scalar>& a_b_mf,
142 Set::Field<Set::Scalar>& a_rhs_mf,
143 Set::Field<Set::Matrix>& a_dw_mf,
144 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>>& a_ddw_mf,
145 Set::Field<T>& a_model_mf)
146 {
147 for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
148 {
149 amrex::Box domain(linop->Geom(lev).growPeriodicDomain(1));
150 domain.convert(amrex::IntVect::TheNodeVector());
151 const Set::Scalar* dx = linop->Geom(lev).CellSize();
152 Set::Vector DX(linop->Geom(lev).CellSize());
153 for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
154 {
155 amrex::Box bx = mfi.grownnodaltilebox();
156 bx = bx & domain;
157
158 amrex::Array4<const T> const& model = a_model_mf[lev]->array(mfi);
159 amrex::Array4<const Set::Scalar> const& u = a_u_mf[lev]->array(mfi);
160 amrex::Array4<Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
161 amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, T::sym>> const& ddw = a_ddw_mf[lev]->array(mfi);
162
163 // Set model internal dw and ddw.
164 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
165 auto sten = Numeric::GetStencil(i, j, k, bx);
166
167 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
168
169 if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
170 {
171 dw(i, j, k) = model(i, j, k).DW(gradu);
172 ddw(i, j, k) = model(i, j, k).DDW(gradu);
173 }
174 else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
175 {
176 Set::Matrix eps = 0.5 * (gradu + gradu.transpose());
177 dw(i, j, k) = model(i, j, k).DW(eps);
178 ddw(i, j, k) = model(i, j, k).DDW(eps);
179 }
180 else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
181 {
182 Set::Matrix F = gradu + Set::Matrix::Identity();
183 dw(i, j, k) = model(i, j, k).DW(F);
184 ddw(i, j, k) = model(i, j, k).DDW(F);
185 }
186 });
187 }
188
189 Util::RealFillBoundary(*a_dw_mf[lev], m_elastic->Geom(lev));
190 Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
191 }
192
193 m_elastic->SetModel(a_ddw_mf);
194
195 for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
196 {
197 amrex::Box domain(linop->Geom(lev).growPeriodicDomain(1));
198 domain.convert(amrex::IntVect::TheNodeVector());
199 const Set::Scalar* dx = linop->Geom(lev).CellSize();
200 const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
201 for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
202 {
203 amrex::Box bx = mfi.nodaltilebox();
204 bx = bx & domain;
205 amrex::Array4<const Set::Scalar> const& u = a_u_mf[lev]->array(mfi);
206 amrex::Array4<const Set::Scalar> const& b = a_b_mf[lev]->array(mfi);
207 amrex::Array4<const Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
208 amrex::Array4<Set::Scalar> const& rhs = a_rhs_mf[lev]->array(mfi);
209 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
210 {
211 auto sten = Numeric::GetStencil(i, j, k, bx);
212 // Do this if on the domain boundary
213 if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
214 {
215 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
216 Set::Vector U(AMREX_D_DECL(u(i, j, k, 0), u(i, j, k, 1), u(i, j, k, 2)));
217 Set::Vector ret = m_elastic->GetBC()(U, gradu, dw(i, j, k), i, j, k, bx);
218 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k, d) - ret(d);
219 }
220 else
221 {
222 Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx);
223 for (int p = 0; p < AMREX_SPACEDIM; p++) rhs(i, j, k, p) = b(i, j, k, p) - divdw(p);
224 }
225 });
226 }
227 Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
228 Util::RealFillBoundary(*a_rhs_mf[lev], m_elastic->Geom(lev));
229 }
230 }
231
232 void prepareForSolve(const Set::Field<Set::Vector>& a_u_mf,
233 const Set::Field<Set::Vector>& a_b_mf,
234 Set::Field<Set::Scalar>& a_rhs_mf,
235 Set::Field<Set::Matrix>& a_dw_mf,
236 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>>& a_ddw_mf,
237 Set::Field<T>& a_model_mf)
238 {
239 const bool conservative_solve = usesConservativeFaceFlux();
240 for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
241 {
242 amrex::Box domain(linop->Geom(lev).growPeriodicDomain(2));
243 domain.convert(amrex::IntVect::TheNodeVector());
244 amrex::Box row_domain(linop->Geom(lev).growPeriodicDomain(1));
245 row_domain.convert(amrex::IntVect::TheNodeVector());
246 const Set::Scalar* dx = linop->Geom(lev).CellSize();
247 Set::Vector DX(linop->Geom(lev).CellSize());
248 const amrex::Dim3 hi = amrex::ubound(domain);
249
250 for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
251 {
252 amrex::Box bx = mfi.grownnodaltilebox();
253 bx = bx & domain;
254
255 amrex::Array4<const T> const& model = a_model_mf[lev]->array(mfi);
256 amrex::Array4<const Set::Vector> const& u = a_u_mf[lev]->array(mfi);
257 amrex::Array4<Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
258 amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, T::sym>> const& ddw = a_ddw_mf[lev]->array(mfi);
259
260 if (conservative_solve)
261 {
262 amrex::Box cbx = mfi.grownnodaltilebox(-1, 1) & row_domain;
263 amrex::LoopConcurrentOnCpu(cbx, [=] (int i, int j, int k)
264 {
265 auto sten = Numeric::GetStencil(i, j, k, cbx);
266 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
267 Set::Matrix kinvar;
268 if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
269 kinvar = gradu;
270 else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
271 kinvar = 0.5 * (gradu + gradu.transpose());
272 else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
273 kinvar = gradu + Set::Matrix::Identity();
274 else
275 Util::Abort(INFO, "Invalid model kinematic variable");
276
277 dw(i, j, k) = model(i, j, k).DW(kinvar);
278 ddw(i, j, k, 0) = model(i, j, k).DDW(kinvar);
279 for (int face = 0; face < AMREX_SPACEDIM; ++face)
280 {
281 const int index[3] = {i, j, k};
282 const int upper[3] = {hi.x, hi.y, hi.z};
283 if (index[face] < upper[face])
284 ddw(i, j, k, face + 1) = Detail::FaceTangent(
285 model, u, i, j, k, face, dx);
286 else
287 ddw(i, j, k, face + 1) = ddw(i, j, k, 0);
288 }
289 });
290 continue;
291 }
292
293 // Set model internal dw and ddw.
294 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
295 {
296 auto sten = Numeric::GetStencil(i, j, k, bx);
297
298 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
299 Set::Matrix kinvar;
300 if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
301 kinvar = gradu; // gradu
302 else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
303 kinvar = 0.5 * (gradu + gradu.transpose()); // epsilon
304 else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
305 kinvar = gradu + Set::Matrix::Identity(); // F
306
307 dw(i, j, k) = model(i, j, k).DW(kinvar);
308 ddw(i, j, k) = model(i, j, k).DDW(kinvar);
309
310 });
311 }
312
313 a_dw_mf[lev]->setMultiGhost(true);
314 a_dw_mf[lev]->FillBoundaryAndSync(m_elastic->Geom(lev).periodicity());
315 a_ddw_mf[lev]->setMultiGhost(true);
316 a_ddw_mf[lev]->FillBoundaryAndSync(m_elastic->Geom(lev).periodicity());
317 }
318
319 m_elastic->SetModel(a_ddw_mf);
320 if (m_psi)
321 for (int i = 0; i <= a_b_mf.finest_level; i++)
322 m_elastic->SetPsi(i, *(*m_psi)[i]);
323
324
325 for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
326 {
327 amrex::Box domain(linop->Geom(lev).growPeriodicDomain(2));
328 domain.convert(amrex::IntVect::TheNodeVector());
329 amrex::Box row_domain(linop->Geom(lev).growPeriodicDomain(1));
330 row_domain.convert(amrex::IntVect::TheNodeVector());
331 const Set::Scalar* dx = linop->Geom(lev).CellSize();
332 const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
333 for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
334 {
335 amrex::Box bx = mfi.grownnodaltilebox() & domain;
336 amrex::Array4<const Set::Vector> const& u = a_u_mf[lev]->array(mfi);
337 amrex::Array4<const Set::Vector> const& b = a_b_mf[lev]->array(mfi);
338 amrex::Array4<const Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
339 amrex::Array4<Set::Scalar> const& rhs = a_rhs_mf[lev]->array(mfi);
340 if (conservative_solve)
341 {
342 amrex::Box cbx = mfi.grownnodaltilebox(-1, 1) & row_domain;
343 amrex::Array4<const T> const& model = a_model_mf[lev]->array(mfi);
344 amrex::LoopConcurrentOnCpu(cbx, [=] (int i, int j, int k)
345 {
346 auto sten = Numeric::GetStencil(i, j, k, cbx);
347 const int index[3] = {i, j, k};
348 const int lower[3] = {lo.x, lo.y, lo.z};
349 const int upper[3] = {hi.x, hi.y, hi.z};
350 bool on_boundary = false;
351 for (int dir = 0; dir < AMREX_SPACEDIM; ++dir)
352 on_boundary = on_boundary ||
353 index[dir] == lower[dir] ||
354 index[dir] == upper[dir];
355 if (on_boundary)
356 {
357 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
358 Set::Vector ret = m_elastic->GetBC()(
359 u(i, j, k), gradu, dw(i, j, k), i, j, k, cbx);
360 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - ret(d);
361 }
362 else
363 {
364 Set::Vector divdw = Detail::FaceStressDivergence(
365 model, u, i, j, k, dx);
366 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - divdw(d);
367 }
368 });
369 continue;
370 }
371
372 // This is for if psi is not being used and has not been set
373 if (m_psi == nullptr)
374 {
375 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
376 {
377 auto sten = Numeric::GetStencil(i, j, k, bx);
378 // Do this if on the domain boundary
379 if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
380 {
381 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
382 Set::Vector ret = m_elastic->GetBC()(u(i, j, k), gradu, dw(i, j, k), i, j, k, bx);
383 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - ret(d);
384 }
385 else
386 {
387 Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx, sten);
388 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - divdw(d);
389 }
390 });
391 }
392 // This is EXACTLY THE SAME AS THE ABOVE, excpet that we are now accounting
393 // for psi terms. The reason we do this is because there's no good way (yet) to
394 // initialize patches if psi = nullptr
395 else
396 {
397 const amrex::Dim3 boxlo = amrex::lbound(bx), boxhi = amrex::ubound(bx);
398 amrex::Array4<const Set::Scalar> const& psi = (*m_psi)[lev]->array(mfi);
399 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
400 {
401 auto sten = Numeric::GetStencil(i, j, k, bx);
402 // Do this if on the domain boundary
403 if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
404 {
405 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
406 Set::Scalar psiavg = 1.0;
407 if (AMREX_D_TERM(i > boxlo.x && i<boxhi.x, && j>boxlo.y && j<boxhi.y, && k>boxlo.z && k < boxhi.z))
408 {
409 psiavg = Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0);
410 }
411 Set::Vector ret = m_elastic->GetBC()(u(i, j, k), gradu, dw(i, j, k) * psiavg, i, j, k, bx);
412 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - ret(d);
413 }
414 else
415 {
416 Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx, sten);
417 if (AMREX_D_TERM(i > boxlo.x && i<boxhi.x, && j>boxlo.y && j<boxhi.y, && k>boxlo.z && k < boxhi.z))
418 {
419 divdw *= Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0);
420 Set::Vector gradpsi = Numeric::CellGradientOnNode(psi, i, j, k, 0, dx);
421 divdw += dw(i, j, k) * gradpsi;
422 }
423 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - divdw(d);
424 }
425 });
426 }
427 }
428
429 a_rhs_mf[lev]->setMultiGhost(true);
430 a_rhs_mf[lev]->FillBoundaryAndSync(m_elastic->Geom(lev).periodicity());
431 a_ddw_mf[lev]->setMultiGhost(true);
432 a_ddw_mf[lev]->FillBoundaryAndSync(m_elastic->Geom(lev).periodicity());
433 }
434 }
435
436
437public:
438 Set::Scalar solve(const Set::Field<Set::Vector>& a_u_mf,
439 const Set::Field<Set::Vector>& a_b_mf,
440 Set::Field<T>& a_model_mf,
441 Real a_tol_rel, Real a_tol_abs, const char* checkpoint_file = nullptr)
442 {
443 if (verbose > 0)
444 Util::Message(INFO, "Newton elastic operator: conservative_face_flux=",
445 usesConservativeFaceFlux() ? 1 : 0);
446
447 Set::Field<Set::Scalar> dsol_mf, rhs_mf, u0_mf, utrial_mf;
448 Set::Field<Set::Matrix> dw_mf;
449 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
450 const int ddw_ncomp = usesConservativeFaceFlux() ? AMREX_SPACEDIM + 1 : 1;
451
452 dsol_mf.resize(a_u_mf.finest_level + 1); dsol_mf.finest_level = a_u_mf.finest_level;
453 dw_mf.resize(a_u_mf.finest_level + 1); dw_mf.finest_level = a_u_mf.finest_level;
454 ddw_mf.resize(a_u_mf.finest_level + 1); ddw_mf.finest_level = a_u_mf.finest_level;
455 rhs_mf.resize(a_u_mf.finest_level + 1); rhs_mf.finest_level = a_u_mf.finest_level;
456 u0_mf.resize(a_u_mf.finest_level + 1); u0_mf.finest_level = a_u_mf.finest_level;
457 utrial_mf.resize(a_u_mf.finest_level + 1); utrial_mf.finest_level = a_u_mf.finest_level;
458 for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
459 {
460 dsol_mf.Define(lev, a_u_mf[lev]->boxArray(),
461 a_u_mf[lev]->DistributionMap(),
462 a_u_mf.NComp(),
463 a_u_mf[lev]->nGrow());
464 dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
465 a_b_mf[lev]->DistributionMap(),
466 1,
467 a_b_mf[lev]->nGrow());
468 ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
469 a_b_mf[lev]->DistributionMap(),
470 ddw_ncomp,
471 a_b_mf[lev]->nGrow());
472 rhs_mf.Define(lev, a_b_mf[lev]->boxArray(),
473 a_b_mf[lev]->DistributionMap(),
474 a_b_mf.NComp(),
475 a_b_mf[lev]->nGrow());
476 u0_mf.Define(lev, a_u_mf[lev]->boxArray(),
477 a_u_mf[lev]->DistributionMap(),
478 a_u_mf.NComp(),
479 a_u_mf[lev]->nGrow());
480 utrial_mf.Define(lev, a_u_mf[lev]->boxArray(),
481 a_u_mf[lev]->DistributionMap(),
482 a_u_mf.NComp(),
483 a_u_mf[lev]->nGrow());
484
485 dsol_mf[lev]->setVal(0.0);
486 dw_mf[lev]->setVal(Set::Matrix::Zero());
487 ddw_mf[lev]->setVal(Set::Matrix4<AMREX_SPACEDIM, T::sym>::Zero());
488
489 a_b_mf.Copy(lev, *rhs_mf[lev], 0, 2);
490 //amrex::MultiFab::Copy(*rhs_mf[lev], *a_b_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
491 }
492
493 for (int nriter = 0; nriter < m_nriters; nriter++)
494 {
495 if (verbose > 0 && nriter < m_nriters) Util::Message(INFO, "Newton Iteration ", nriter + 1, " of ", m_nriters);
496
497 prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf, ddw_mf, a_model_mf);
498
499 if (nriter == m_nriters) break;
500
501 if (m_resync_coeffs && nriter > 0)
502 m_elastic->SyncCoefficients();
503
504 const Set::Scalar resnorm0 =
505 m_line_search ? FieldNorm0(rhs_mf) : -1.0;
506
507 for (int lev = 0; lev < dsol_mf.size(); ++lev)
508 dsol_mf[lev]->setVal(0.0);
509 Solver::Nonlocal::Linear::solve(dsol_mf, rhs_mf, a_tol_rel, a_tol_abs, checkpoint_file);
510
511 Set::Scalar cornorm = 0;
512 for (int lev = 0; lev < dsol_mf.size(); ++lev)
513 {
514 for (int comp = 0; comp < AMREX_SPACEDIM; comp++)
515 {
516 Set::Scalar tmpcornorm = dsol_mf[lev]->norm0(comp, 0);
517 if (tmpcornorm > cornorm) cornorm = tmpcornorm;
518 }
519 }
520
521 Set::Scalar alpha = 1.0;
522 Set::Scalar resnorm = resnorm0;
523 if (m_line_search)
524 {
525 for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
526 a_u_mf.Copy(lev, *u0_mf[lev], 0, 2);
527
528 auto restore_baseline = [&]()
529 {
530 for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
531 a_u_mf.CopyFrom(lev, *u0_mf[lev], 0, 2);
532 };
533 constexpr int max_backtrack = 8;
534 bool accepted = false;
535 for (int bt = 0; bt <= max_backtrack; bt++)
536 {
537 for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
538 {
539 amrex::MultiFab::Copy(*utrial_mf[lev], *u0_mf[lev],
540 0, 0, AMREX_SPACEDIM, 2);
541 amrex::MultiFab::Saxpy(*utrial_mf[lev], alpha,
542 *dsol_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
543 a_u_mf.CopyFrom(lev, *utrial_mf[lev], 0, 2);
544 }
545
546 try
547 {
548 prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf,
549 ddw_mf, a_model_mf);
550 }
551 catch (...)
552 {
553 restore_baseline();
554 prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf,
555 ddw_mf, a_model_mf);
556 throw;
557 }
558
559 resnorm = FieldNorm0(rhs_mf);
560
561 // Permit a small residual increase for roundoff-level
562 // changes in the composite AMR norm.
563 if (std::isfinite(resnorm0) && std::isfinite(resnorm) &&
564 (resnorm0 == 0.0 ? resnorm == 0.0 :
565 resnorm / resnorm0 <= 1.0001))
566 {
567 accepted = true;
568 break;
569 }
570
571 alpha *= 0.5;
572 }
573 if (!accepted)
574 {
575 restore_baseline();
576 prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf,
577 ddw_mf, a_model_mf);
578 Util::Abort(INFO, "Newton line search failed after ",
579 max_backtrack, " backtracks: initial residual=", resnorm0,
580 ", final trial residual=", resnorm);
581 }
582 }
583 else
584 {
585 for (int lev = 0; lev < dsol_mf.size(); ++lev)
586 a_u_mf.AddFrom(lev, *dsol_mf[lev], 0, 2);
587 //amrex::MultiFab::Add(*a_u_mf[lev], *dsol_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
588 }
589
590 const Set::Scalar full_update = cornorm;
591 const Set::Scalar accepted_update = alpha * full_update;
592 if (!std::isfinite(full_update) || !std::isfinite(accepted_update))
593 Util::Abort(INFO, "Newton update norm is non-finite");
594 if (verbose > 0)
595 {
596 Util::Message(INFO, "NR iteration ", nriter + 1,
597 ", alpha = ", alpha,
598 ", full max norm(ddisp) = ", full_update,
599 ", accepted max norm(ddisp) = ", accepted_update);
600 if (m_line_search)
601 Util::Message(INFO, "NR residual: initial = ", resnorm0,
602 ", accepted = ", resnorm);
603 }
604
605 // Test the undamped Newton correction so backtracking cannot
606 // manufacture convergence merely by choosing a very small alpha.
607 if (m_nrtolerance > 0.0 && full_update < m_nrtolerance)
608 return full_update;
609 if (m_nrtolerance > 0.0 && nriter + 1 == m_nriters)
610 {
611 if (m_nr_abort_on_fail)
612 Util::Abort(INFO, "Newton solver failed to converge in ",
613 m_nriters, " iterations: final full update=", full_update,
614 ", final accepted update=", accepted_update,
615 ", tolerance=", m_nrtolerance);
616 Util::Warning(INFO, "Newton solver failed to converge in ",
617 m_nriters, " iterations: final full update=", full_update,
618 ", final accepted update=", accepted_update,
619 ", tolerance=", m_nrtolerance);
620 return full_update;
621 }
622
623 }
624
625 return 0.0;
626 }
627
628 Set::Scalar solve(const Set::Field<Set::Scalar>& a_u_mf,
629 const Set::Field<Set::Scalar>& a_b_mf,
630 Set::Field<T>& a_model_mf,
631 Real a_tol_rel, Real a_tol_abs, const char* checkpoint_file = nullptr)
632 {
633 Set::Field<Set::Scalar> dsol_mf, rhs_mf;
634 Set::Field<Set::Matrix> dw_mf;
635 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
636
637 dsol_mf.resize(a_u_mf.finest_level + 1); dsol_mf.finest_level = a_u_mf.finest_level;
638 dw_mf.resize(a_u_mf.finest_level + 1); dw_mf.finest_level = a_u_mf.finest_level;
639 ddw_mf.resize(a_u_mf.finest_level + 1); ddw_mf.finest_level = a_u_mf.finest_level;
640 rhs_mf.resize(a_u_mf.finest_level + 1); rhs_mf.finest_level = a_u_mf.finest_level;
641 for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
642 {
643 dsol_mf.Define(lev, a_u_mf[lev]->boxArray(),
644 a_u_mf[lev]->DistributionMap(),
645 a_u_mf[lev]->nComp(),
646 a_u_mf[lev]->nGrow());
647 dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
648 a_b_mf[lev]->DistributionMap(),
649 1,
650 a_b_mf[lev]->nGrow());
651 ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
652 a_b_mf[lev]->DistributionMap(),
653 1,
654 a_b_mf[lev]->nGrow());
655 rhs_mf.Define(lev, a_b_mf[lev]->boxArray(),
656 a_b_mf[lev]->DistributionMap(),
657 a_b_mf[lev]->nComp(),
658 a_b_mf[lev]->nGrow());
659
660 dsol_mf[lev]->setVal(0.0);
661 dw_mf[lev]->setVal(Set::Matrix::Zero());
662 ddw_mf[lev]->setVal(Set::Matrix4<AMREX_SPACEDIM, T::sym>::Zero());
663
664 amrex::MultiFab::Copy(*rhs_mf[lev], *a_b_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
665 }
666
667 for (int nriter = 0; nriter < m_nriters; nriter++)
668 {
669 if (verbose > 0 && nriter < m_nriters) Util::Message(INFO, "Newton Iteration ", nriter + 1, " of ", m_nriters);
670
671 prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf, ddw_mf, a_model_mf);
672
673
674 if (nriter == m_nriters) break;
675
676 Solver::Nonlocal::Linear::solve(dsol_mf, rhs_mf, a_tol_rel, a_tol_abs, checkpoint_file);
677
678 Set::Scalar cornorm = 0, solnorm = 0;
679 for (int lev = 0; lev < dsol_mf.size(); ++lev)
680 {
681 for (int comp = 0; comp < AMREX_SPACEDIM; comp++)
682 {
683 Set::Scalar tmpcornorm = dsol_mf[lev]->norm0(comp, 0);
684 if (tmpcornorm > cornorm) cornorm = tmpcornorm;
685
686 Set::Scalar tmpsolnorm = a_u_mf[lev]->norm0(comp, 0);
687 if (tmpsolnorm > solnorm) solnorm = tmpsolnorm;
688 }
689
690 }
691 Set::Scalar relnorm;
692 if (solnorm == 0) relnorm = cornorm;
693 else relnorm = cornorm / solnorm;
694 if (verbose > 0) Util::Message(INFO, "NR iteration ", nriter + 1, ", relative norm(ddisp) = ", relnorm);
695
696 for (int lev = 0; lev < dsol_mf.size(); ++lev)
697 amrex::MultiFab::Add(*a_u_mf[lev], *dsol_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
698
699 if (relnorm < m_nrtolerance)
700 return relnorm;
701
702 }
703
704 return 0.0;
705 }
706 Set::Scalar solve(const Set::Field<Set::Scalar>& a_u_mf,
707 const Set::Field<Set::Scalar>& a_b_mf,
708 Set::Field<T>& a_model_mf)
709 {
710 return solve(a_u_mf, a_b_mf, a_model_mf, tol_rel, tol_abs);
711 }
712
713 void compResidual(Set::Field<Set::Scalar>& a_res_mf,
714 Set::Field<Set::Scalar>& a_u_mf,
715 Set::Field<Set::Scalar>& a_b_mf,
716 Set::Field<T>& a_model_mf)
717 {
718 Set::Field<Set::Matrix> dw_mf;
719 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
720 dw_mf.resize(a_u_mf.size());
721 ddw_mf.resize(a_u_mf.size());
722 for (int lev = 0; lev < a_u_mf.size(); lev++)
723 {
724 dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
725 a_b_mf[lev]->DistributionMap(),
726 1, a_b_mf[lev]->nGrow());
727 ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
728 a_b_mf[lev]->DistributionMap(),
729 1, a_b_mf[lev]->nGrow());
730 dw_mf[lev]->setVal(Set::Matrix::Zero());
731 }
732
733 prepareForSolve(a_u_mf, a_b_mf, a_res_mf, dw_mf, ddw_mf, a_model_mf);
734 }
735
736 void compResidual(Set::Field<Set::Vector>& a_res_mf,
737 Set::Field<Set::Vector>& a_u_mf,
738 Set::Field<Set::Vector>& a_b_mf,
739 Set::Field<T>& a_model_mf)
740 {
741 Set::Field<Set::Matrix> dw_mf;
742 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
743 Set::Field<Set::Scalar> res_mf(a_res_mf.size());
744 dw_mf.resize(a_u_mf.size());
745 ddw_mf.resize(a_u_mf.size());
746 res_mf.resize(a_u_mf.size());
747 const int ddw_ncomp = usesConservativeFaceFlux() ? AMREX_SPACEDIM + 1 : 1;
748 for (int lev = 0; lev < a_u_mf.size(); lev++)
749 {
750 dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
751 a_b_mf[lev]->DistributionMap(),
752 1, a_b_mf[lev]->nGrow());
753 ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
754 a_b_mf[lev]->DistributionMap(),
755 ddw_ncomp, a_b_mf[lev]->nGrow());
756 res_mf.Define(lev, a_b_mf[lev]->boxArray(),
757 a_b_mf[lev]->DistributionMap(),
758 AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
759 dw_mf[lev]->setVal(Set::Matrix::Zero());
760 }
761
762 prepareForSolve(a_u_mf, a_b_mf, res_mf, dw_mf, ddw_mf, a_model_mf);
763
764 if (usesConservativeFaceFlux())
765 for (int lev = a_u_mf.finest_level - 1; lev >= 0; --lev)
766 m_elastic->Reflux(lev, *res_mf[lev], *res_mf[lev], *res_mf[lev],
767 *res_mf[lev + 1], *res_mf[lev + 1], *res_mf[lev + 1]);
768
769 for (int lev = 0; lev < a_res_mf.size(); ++lev)
770 {
771 Util::RealFillBoundary(*res_mf[lev], m_elastic->Geom(lev));
772 a_res_mf.CopyFrom(lev, *res_mf[lev], 0, 2);
773 }
774 }
775
776 void compLinearSolverResidual(Set::Field<Set::Vector>& a_res_mf,
777 Set::Field<Set::Vector>& a_u_mf,
778 Set::Field<Set::Vector>& a_b_mf)
779 {
780 Util::Assert(INFO, TEST(a_res_mf.finest_level == a_u_mf.finest_level));
781 Util::Assert(INFO, TEST(a_u_mf.finest_level == a_b_mf.finest_level));
782 Set::Field<Set::Scalar> res_mf(a_res_mf.finest_level + 1), sol_mf(a_u_mf.finest_level + 1), rhs_mf(a_b_mf.finest_level + 1);
783 for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
784 {
785 res_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
786 sol_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
787 rhs_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
788
789 a_u_mf.Copy(lev, *sol_mf[lev], 0, 2);
790 a_b_mf.Copy(lev, *rhs_mf[lev], 0, 2);
791 }
792
793 mlmg->compResidual(amrex::GetVecOfPtrs(res_mf), amrex::GetVecOfPtrs(sol_mf), amrex::GetVecOfConstPtrs(rhs_mf));
794
795 for (int lev = 0; lev <= a_res_mf.finest_level; ++lev)
796 {
797 a_res_mf.CopyFrom(lev, *res_mf[lev], 0, 2);
798 }
799 }
800
801
802 void W(Set::Field<Set::Scalar>& a_w_mf,
803 Set::Field<Set::Scalar>& a_u_mf,
804 Set::Field<T>& a_model_mf)
805 {
806 for (int lev = 0; lev < a_u_mf.size(); lev++)
807 {
808 BL_PROFILE("Solver::Nonlocal::Newton::DW()");
809
810 const amrex::Real* DX = linop->Geom(lev).CellSize();
811 amrex::Box domain(linop->Geom(lev).Domain());
812 domain.convert(amrex::IntVect::TheNodeVector());
813
814 for (MFIter mfi(*a_u_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
815 {
816 const Box& bx = mfi.tilebox();
817 amrex::Array4<T> const& C = a_model_mf[lev]->array(mfi);
818 amrex::Array4<amrex::Real> const& w = a_w_mf[lev]->array(mfi);
819 amrex::Array4<const amrex::Real> const& u = a_u_mf[lev]->array(mfi);
820 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
821 {
822 Set::Matrix gradu;
823
824 auto sten = Numeric::GetStencil(i, j, k, bx);
825
826 // Fill gradu
827 for (int p = 0; p < AMREX_SPACEDIM; p++)
828 {
829 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
830 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
831 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
832 }
833
834 if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
835 w(i, j, k) = C(i, j, k).W(gradu);
836 else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
837 w(i, j, k) = C(i, j, k).W(0.5 * (gradu + gradu.transpose()));
838 else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
839 w(i, j, k) = C(i, j, k).W(gradu + Set::Matrix::Identity());
840 });
841 }
842 }
843 }
844
845 void DW(Set::Field<Set::Scalar>& a_dw_mf,
846 Set::Field<Set::Scalar>& a_u_mf,
847 Set::Field<T>& a_model_mf)
848 {
849 for (int lev = 0; lev < a_u_mf.size(); lev++)
850 {
851 BL_PROFILE("Solver::Nonlocal::Newton::DW()");
852
853 const amrex::Real* DX = linop->Geom(lev).CellSize();
854 amrex::Box domain(linop->Geom(lev).Domain());
855 domain.convert(amrex::IntVect::TheNodeVector());
856
857 for (MFIter mfi(*a_u_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
858 {
859 const Box& bx = mfi.tilebox();
860 amrex::Array4<T> const& C = a_model_mf[lev]->array(mfi);
861 amrex::Array4<amrex::Real> const& dw = a_dw_mf[lev]->array(mfi);
862 amrex::Array4<const amrex::Real> const& u = a_u_mf[lev]->array(mfi);
863 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
864 {
865 Set::Matrix gradu;
866
867 auto sten = Numeric::GetStencil(i, j, k, bx);
868
869 // Fill gradu
870 for (int p = 0; p < AMREX_SPACEDIM; p++)
871 {
872 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
873 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
874 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
875 }
876
877 Set::Matrix sig = Set::Matrix::Zero();
878
879 if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
880 sig = C(i, j, k).DW(gradu);
881 else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
882 sig = C(i, j, k).DW(0.5 * (gradu + gradu.transpose()));
883 else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
884 sig = C(i, j, k).DW(gradu + Set::Matrix::Identity());
885
886 // = C(i,j,k)(gradu,m_homogeneous);
887
888 AMREX_D_PICK(dw(i, j, k, 0) = sig(0, 0);
889 ,
890 dw(i, j, k, 0) = sig(0, 0); dw(i, j, k, 1) = sig(0, 1);
891 dw(i, j, k, 2) = sig(1, 0); dw(i, j, k, 3) = sig(1, 1);
892 ,
893 dw(i, j, k, 0) = sig(0, 0); dw(i, j, k, 1) = sig(0, 1); dw(i, j, k, 2) = sig(0, 2);
894 dw(i, j, k, 3) = sig(1, 0); dw(i, j, k, 4) = sig(1, 1); dw(i, j, k, 5) = sig(1, 2);
895 dw(i, j, k, 6) = sig(2, 0); dw(i, j, k, 7) = sig(2, 1); dw(i, j, k, 8) = sig(2, 2););
896
897 });
898 }
899 }
900 }
901
902
903public:
904 int m_nriters = 1;
905 Set::Scalar m_nrtolerance = 0.0;
906 bool m_line_search = false;
907 bool m_resync_coeffs = false;
908 bool m_nr_abort_on_fail = false;
909 Operator::Elastic<T::sym>* m_elastic;
910 //BC::Operator::Elastic::Elastic *m_bc;
911
912 Set::Field<Set::Scalar>* m_psi = nullptr;
913 bool m_conservative_face_flux = false;
914
915private:
916 bool m_reported_psi_face_flux_conflict = false;
917
918 Set::Scalar FieldNorm0(const Set::Field<Set::Scalar>& a_mf) const
919 {
920 amrex::Vector<std::unique_ptr<amrex::MultiFab>> composite(a_mf.size());
921 for (int lev = 0; lev < a_mf.size(); ++lev)
922 {
923 composite[lev] = std::make_unique<amrex::MultiFab>(
924 a_mf[lev]->boxArray(), a_mf[lev]->DistributionMap(),
925 a_mf[lev]->nComp(), a_mf[lev]->nGrowVect());
926 amrex::MultiFab::Copy(*composite[lev], *a_mf[lev], 0, 0,
927 a_mf[lev]->nComp(), a_mf[lev]->nGrowVect());
928 }
929
930 for (int lev = static_cast<int>(composite.size()) - 2; lev >= 0; --lev)
931 m_elastic->Reflux(lev, *composite[lev], *composite[lev],
932 *composite[lev], *composite[lev + 1], *composite[lev + 1],
933 *composite[lev + 1]);
934
935 Set::Scalar norm = 0.0;
936 for (int lev = 0; lev < composite.size(); ++lev)
937 {
938 const Set::Scalar tmp =
939 m_elastic->normInf(lev, *composite[lev], false);
940 if (tmp > norm) norm = tmp;
941 }
942 return norm;
943 }
944
945public:
946 // These paramters control a standard Newton-Raphson solve.
947 //
948 // **Note**:
949 // This class inherits all of the linear solve paramters
950 // from its parent class, :ref:`Solver::Nonlocal::Linear`
951 static void Parse(Newton<T>& value, amrex::ParmParse& pp)
952 {
953 Linear::Parse(value, pp);
954
955 // Number of newton-raphson iterations.
956 pp_query("nriters", value.m_nriters);
957
958 // Tolerance for the undamped absolute correction norm in vector
959 // mechanics solves, and the relative correction norm in scalar solves.
960 pp_query("nrtolerance", value.m_nrtolerance);
961
962 // Backtrack nonlinear updates until the residual is within a small
963 // relative allowance of its initial value.
964 pp_query("line_search", value.m_line_search);
965
966 // Rebuild coarse MG coefficients and the smoother diagonal after each
967 // nonlinear relinearization. MLMG otherwise prepares them only once.
968 pp_query("resync_coeffs", value.m_resync_coeffs);
969
970 // Abort instead of warning if Newton iterations do not converge.
971 pp_query("nr_abort_on_fail", value.m_nr_abort_on_fail);
972
973 // Discretize the elastic residual with coefficient-aware face fluxes.
974 pp_query("conservative_face_flux", value.m_conservative_face_flux);
975 }
976
977};
978} // namespace Nonlocal
979} // namespace Solver
980
981
982#endif
pp_query
#define pp_query(...)
Definition ParmParse.H:129
TEST
#define TEST(x)
Definition Util.H:25
INFO
#define INFO
Definition Util.H:24
Operator::Elastic::SetModel
void SetModel(Set::Matrix4< AMREX_SPACEDIM, SYM > &a_model)
Definition Elastic.cpp:67
Operator::Elastic::GetBC
::BC::Operator::Elastic::Elastic & GetBC()
Definition Elastic.H:68
Operator::Elastic::SetPsi
void SetPsi(int amrlev, const amrex::MultiFab &a_psi)
Definition Elastic.cpp:137
Operator::Operator< Grid::Node >::Reflux
void Reflux(int crse_amrlev, MultiFab &res, const MultiFab &crse_sol, const MultiFab &crse_rhs, MultiFab &fine_res, MultiFab &fine_sol, const MultiFab &fine_rhs)
Definition Operator.H:51
Operator::Operator< Grid::Node >::Geom
const Geometry & Geom(int amr_lev, int mglev=0) const noexcept
Definition Operator.H:49
Set::Field::Copy
void Copy(int, amrex::MultiFab &, int, int) const
Definition Set.H:68
Set::Field::AddFrom
void AddFrom(int, amrex::MultiFab &, int, int) const
Definition Set.H:71
Set::Field::finest_level
int finest_level
Definition Set.H:67
Set::Field::Define
void Define(int a_levs, const amrex::Vector< amrex::BoxArray > &a_grids, const amrex::Vector< amrex::DistributionMapping > &a_dmap, int a_ncomp, int a_nghost)
Definition Set.H:52
Set::Field::NComp
int NComp() const
Definition Set.H:72
Set::Field::CopyFrom
void CopyFrom(int, amrex::MultiFab &, int, int) const
Definition Set.H:69
Solver::Nonlocal::Linear
Multigrid Linear solver for multicomponent, multi-level operators.
Definition Linear.H:19
Solver::Nonlocal::Linear::solve
Set::Scalar solve(amrex::Vector< std::unique_ptr< amrex::MultiFab > > &a_sol, amrex::Vector< std::unique_ptr< amrex::MultiFab > > &a_rhs, Real a_tol_rel, Real a_tol_abs, const char *checkpoint_file=nullptr)
Definition Linear.H:131
Solver::Nonlocal::Linear::tol_rel
Set::Scalar tol_rel
Definition Linear.H:246
Solver::Nonlocal::Linear::Parse
static void Parse(Linear &value, amrex::ParmParse &pp)
Definition Linear.H:296
Solver::Nonlocal::Linear::Define
void Define(Operator::Operator< Grid::Node > &a_lp)
Definition Linear.H:36
Solver::Nonlocal::Linear::mlmg
amrex::MLMG * mlmg
Definition Linear.H:253
Solver::Nonlocal::Linear::tol_abs
Set::Scalar tol_abs
Definition Linear.H:247
Solver::Nonlocal::Linear::linop
Operator::Operator< Grid::Node > * linop
Definition Linear.H:252
Solver::Nonlocal::Newton::setPsi
void setPsi(Set::Field< Set::Scalar > &a_psi)
Definition Newton.H:124
Solver::Nonlocal::Newton::solve
Set::Scalar solve(const Set::Field< Set::Scalar > &a_u_mf, const Set::Field< Set::Scalar > &a_b_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:706
Solver::Nonlocal::Newton::m_reported_psi_face_flux_conflict
bool m_reported_psi_face_flux_conflict
Definition Newton.H:916
Solver::Nonlocal::Newton::Newton
Newton(Operator::Elastic< T::sym > &a_op)
Definition Newton.H:98
Solver::Nonlocal::Newton::solve
Set::Scalar solve(const Set::Field< Set::Scalar > &a_u_mf, const Set::Field< Set::Scalar > &a_b_mf, Set::Field< T > &a_model_mf, Real a_tol_rel, Real a_tol_abs, const char *checkpoint_file=nullptr)
Definition Newton.H:628
Solver::Nonlocal::Newton::m_nrtolerance
Set::Scalar m_nrtolerance
Definition Newton.H:905
Solver::Nonlocal::Newton::compResidual
void compResidual(Set::Field< Set::Vector > &a_res_mf, Set::Field< Set::Vector > &a_u_mf, Set::Field< Set::Vector > &a_b_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:736
Solver::Nonlocal::Newton::FieldNorm0
Set::Scalar FieldNorm0(const Set::Field< Set::Scalar > &a_mf) const
Definition Newton.H:918
Solver::Nonlocal::Newton::Define
void Define(Operator::Elastic< T::sym > &a_op)
Definition Newton.H:108
Solver::Nonlocal::Newton::usesConservativeFaceFlux
bool usesConservativeFaceFlux() const
Definition Newton.H:134
Solver::Nonlocal::Newton::m_psi
Set::Field< Set::Scalar > * m_psi
Definition Newton.H:912
Solver::Nonlocal::Newton::W
void W(Set::Field< Set::Scalar > &a_w_mf, Set::Field< Set::Scalar > &a_u_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:802
Solver::Nonlocal::Newton::prepareForSolve
void prepareForSolve(const Set::Field< Set::Scalar > &a_u_mf, const Set::Field< Set::Scalar > &a_b_mf, Set::Field< Set::Scalar > &a_rhs_mf, Set::Field< Set::Matrix > &a_dw_mf, Set::Field< Set::Matrix4< AMREX_SPACEDIM, T::sym > > &a_ddw_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:140
Solver::Nonlocal::Newton::Parse
static void Parse(Newton< T > &value, amrex::ParmParse &pp)
Definition Newton.H:951
Solver::Nonlocal::Newton::m_elastic
Operator::Elastic< T::sym > * m_elastic
Definition Newton.H:909
Solver::Nonlocal::Newton::DW
void DW(Set::Field< Set::Scalar > &a_dw_mf, Set::Field< Set::Scalar > &a_u_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:845
Solver::Nonlocal::Newton::prepareForSolve
void prepareForSolve(const Set::Field< Set::Vector > &a_u_mf, const Set::Field< Set::Vector > &a_b_mf, Set::Field< Set::Scalar > &a_rhs_mf, Set::Field< Set::Matrix > &a_dw_mf, Set::Field< Set::Matrix4< AMREX_SPACEDIM, T::sym > > &a_ddw_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:232
Solver::Nonlocal::Newton::compLinearSolverResidual
void compLinearSolverResidual(Set::Field< Set::Vector > &a_res_mf, Set::Field< Set::Vector > &a_u_mf, Set::Field< Set::Vector > &a_b_mf)
Definition Newton.H:776
Solver::Nonlocal::Newton::setNRIters
void setNRIters(int a_nriters)
Definition Newton.H:122
Solver::Nonlocal::Newton::compResidual
void compResidual(Set::Field< Set::Scalar > &a_res_mf, Set::Field< Set::Scalar > &a_u_mf, Set::Field< Set::Scalar > &a_b_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:713
Solver::Nonlocal::Newton::solve
Set::Scalar solve(const Set::Field< Set::Vector > &a_u_mf, const Set::Field< Set::Vector > &a_b_mf, Set::Field< T > &a_model_mf, Real a_tol_rel, Real a_tol_abs, const char *checkpoint_file=nullptr)
Definition Newton.H:438
Numeric::CellGradientOnNode
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Set::Vector CellGradientOnNode(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition Stencil.H:705
Numeric::FaceGradient
Set::Matrix FaceGradient(const amrex::Array4< const Set::Vector > &f, const int i, const int j, const int k, const int face, const Set::Scalar dx[AMREX_SPACEDIM])
Definition Stencil.H:796
Numeric::Divergence
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Set::Vector Divergence(const amrex::Array4< const Set::Matrix > &dw, const int &i, const int &j, const int &k, const Set::Scalar DX[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType())
Definition Stencil.H:596
Numeric::GetStencil
static AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE std::array< StencilType, AMREX_SPACEDIM > GetStencil(const int i, const int j, const int k, const amrex::Box domain)
Definition Stencil.H:51
Numeric::Gradient
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Set::Vector Gradient(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType())
Definition Stencil.H:687
Set::Scalar
amrex::Real Scalar
Definition Base.H:19
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:21
Set::Matrix
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, AMREX_SPACEDIM > Matrix
Definition Base.H:24
Solver::Nonlocal::Detail::FaceStress
Set::Matrix FaceStress(const amrex::Array4< const T > &model, const U &u, const int i, const int j, const int k, const int face, const Set::Scalar dx[AMREX_SPACEDIM])
Definition Newton.H:49
Solver::Nonlocal::Detail::FaceModel
T FaceModel(const amrex::Array4< const T > &model, const int i, const int j, const int k, const int face)
Definition Newton.H:19
Solver::Nonlocal::Detail::FaceTangent
Set::Matrix4< AMREX_SPACEDIM, T::sym > FaceTangent(const amrex::Array4< const T > &model, const U &u, const int i, const int j, const int k, const int face, const Set::Scalar dx[AMREX_SPACEDIM])
Definition Newton.H:61
Solver::Nonlocal::Detail::FaceStressDivergence
Set::Vector FaceStressDivergence(const amrex::Array4< const T > &model, const U &u, const int i, const int j, const int k, const Set::Scalar dx[AMREX_SPACEDIM])
Definition Newton.H:73
Solver::Nonlocal::Detail::FaceKinematicVariable
Set::Matrix FaceKinematicVariable(const T &model, const U &u, const int i, const int j, const int k, const int face, const Set::Scalar dx[AMREX_SPACEDIM])
Definition Newton.H:32
Solver
A bunch of solvers.
Definition CG.H:6
Util::Assert
AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE void Assert(const char *file, const char *func, int line, const char *smt, bool pass, Args const &... args)
Definition Util.H:60
Util::Warning
void Warning(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:213
Util::Message
void Message(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:140
Util::Abort
AMREX_GPU_HOST_DEVICE void Abort(const char *msg)
Definition Util.cpp:406
Util::RealFillBoundary
AMREX_FORCE_INLINE void RealFillBoundary(amrex::FabArray< amrex::BaseFab< T > > &a_mf, const amrex::Geometry &, const int nghost=2)
Definition Util.H:350
Numeric::Interpolate::CellToNodeAverage
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE T CellToNodeAverage(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType())
Definition Stencil.H:1445
Generated by doxygen 1.9.8