Alamo
Newton.H
Go to the documentation of this file.
1#ifndef SOLVER_NONLOCAL_NEWTON
2#define SOLVER_NONLOCAL_NEWTON
3
4#include "Set/Set.H"
5#include "Operator/Elastic.H"
6#include "Solver/Nonlocal/Linear.H"
7#include "IO/ParmParse.H"
8#include "Numeric/Stencil.H"
9
10namespace Solver
11{
12namespace Nonlocal
13{
14template <typename T>
15class Newton: public Linear
16{
17public:
18 Newton() {};
19
20 Newton(Operator::Elastic<T::sym>& a_op)
21 {
22 this->Define(a_op);
23 };
24
25 ~Newton()
26 {
27 if (m_defined) Clear();
28 }
29
30 void Define(Operator::Elastic<T::sym>& a_op)
31 {
32 Linear::Define(a_op);
33 m_elastic = dynamic_cast<Operator::Elastic<T::sym> *>(linop);
34 //m_bc = &m_elastic->GetBC();
35 }
36 void Clear()
37 {
38 Linear::Clear();
39 m_elastic = nullptr;
40 //m_bc = nullptr;
41 }
42
43
44 void setNRIters(int a_nriters) { m_nriters = a_nriters; }
45
46 void setPsi(Set::Field<Set::Scalar>& a_psi)
47 {
48 m_psi = &a_psi;
49 }
50
51private:
52 void prepareForSolve(const Set::Field<Set::Scalar>& a_u_mf,
53 const Set::Field<Set::Scalar>& a_b_mf,
54 Set::Field<Set::Scalar>& a_rhs_mf,
55 Set::Field<Set::Matrix>& a_dw_mf,
56 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>>& a_ddw_mf,
57 Set::Field<T>& a_model_mf)
58 {
59 for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
60 {
61 amrex::Box domain(linop->Geom(lev).growPeriodicDomain(1));
62 domain.convert(amrex::IntVect::TheNodeVector());
63 const Set::Scalar* dx = linop->Geom(lev).CellSize();
64 Set::Vector DX(linop->Geom(lev).CellSize());
65 for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
66 {
67 amrex::Box bx = mfi.grownnodaltilebox();
68 bx = bx & domain;
69
70 amrex::Array4<const T> const& model = a_model_mf[lev]->array(mfi);
71 amrex::Array4<const Set::Scalar> const& u = a_u_mf[lev]->array(mfi);
72 amrex::Array4<Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
73 amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, T::sym>> const& ddw = a_ddw_mf[lev]->array(mfi);
74
75 // Set model internal dw and ddw.
76 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
77 auto sten = Numeric::GetStencil(i, j, k, bx);
78
79 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
80
81 if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
82 {
83 dw(i, j, k) = model(i, j, k).DW(gradu);
84 ddw(i, j, k) = model(i, j, k).DDW(gradu);
85 }
86 else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
87 {
88 Set::Matrix eps = 0.5 * (gradu + gradu.transpose());
89 dw(i, j, k) = model(i, j, k).DW(eps);
90 ddw(i, j, k) = model(i, j, k).DDW(eps);
91 }
92 else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
93 {
94 Set::Matrix F = gradu + Set::Matrix::Identity();
95 dw(i, j, k) = model(i, j, k).DW(F);
96 ddw(i, j, k) = model(i, j, k).DDW(F);
97 }
98 });
99 }
100
101 Util::RealFillBoundary(*a_dw_mf[lev], m_elastic->Geom(lev));
102 Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
103 }
104
105 m_elastic->SetModel(a_ddw_mf);
106
107 for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
108 {
109 amrex::Box domain(linop->Geom(lev).growPeriodicDomain(1));
110 domain.convert(amrex::IntVect::TheNodeVector());
111 const Set::Scalar* dx = linop->Geom(lev).CellSize();
112 const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
113 for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
114 {
115 amrex::Box bx = mfi.nodaltilebox();
116 bx = bx & domain;
117 amrex::Array4<const Set::Scalar> const& u = a_u_mf[lev]->array(mfi);
118 amrex::Array4<const Set::Scalar> const& b = a_b_mf[lev]->array(mfi);
119 amrex::Array4<const Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
120 amrex::Array4<Set::Scalar> const& rhs = a_rhs_mf[lev]->array(mfi);
121 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
122 {
123 auto sten = Numeric::GetStencil(i, j, k, bx);
124 // Do this if on the domain boundary
125 if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
126 {
127 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
128 Set::Vector U(AMREX_D_DECL(u(i, j, k, 0), u(i, j, k, 1), u(i, j, k, 2)));
129 Set::Vector ret = m_elastic->GetBC()(U, gradu, dw(i, j, k), i, j, k, bx);
130 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k, d) - ret(d);
131 }
132 else
133 {
134 Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx);
135 for (int p = 0; p < AMREX_SPACEDIM; p++) rhs(i, j, k, p) = b(i, j, k, p) - divdw(p);
136 }
137 });
138 }
139 Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
140 Util::RealFillBoundary(*a_rhs_mf[lev], m_elastic->Geom(lev));
141 }
142 }
143
144 void prepareForSolve(const Set::Field<Set::Vector>& a_u_mf,
145 const Set::Field<Set::Vector>& a_b_mf,
146 Set::Field<Set::Scalar>& a_rhs_mf,
147 Set::Field<Set::Matrix>& a_dw_mf,
148 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>>& a_ddw_mf,
149 Set::Field<T>& a_model_mf)
150 {
151 for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
152 {
153 amrex::Box domain(linop->Geom(lev).growPeriodicDomain(2));
154 domain.convert(amrex::IntVect::TheNodeVector());
155 const Set::Scalar* dx = linop->Geom(lev).CellSize();
156 Set::Vector DX(linop->Geom(lev).CellSize());
157 for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
158 {
159 amrex::Box bx = mfi.grownnodaltilebox();
160 bx = bx & domain;
161
162 amrex::Array4<const T> const& model = a_model_mf[lev]->array(mfi);
163 amrex::Array4<const Set::Vector> const& u = a_u_mf[lev]->array(mfi);
164 amrex::Array4<Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
165 amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, T::sym>> const& ddw = a_ddw_mf[lev]->array(mfi);
166
167 // Set model internal dw and ddw.
168 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
169 {
170 auto sten = Numeric::GetStencil(i, j, k, bx);
171
172 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
173 Set::Matrix kinvar;
174 if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
175 kinvar = gradu; // gradu
176 else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
177 kinvar = 0.5 * (gradu + gradu.transpose()); // epsilon
178 else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
179 kinvar = gradu + Set::Matrix::Identity(); // F
180
181 dw(i, j, k) = model(i, j, k).DW(kinvar);
182 ddw(i, j, k) = model(i, j, k).DDW(kinvar);
183
184 });
185 }
186
187 a_dw_mf[lev]->setMultiGhost(true);
188 a_dw_mf[lev]->FillBoundaryAndSync(m_elastic->Geom(lev).periodicity());
189 a_ddw_mf[lev]->setMultiGhost(true);
190 a_ddw_mf[lev]->FillBoundaryAndSync(m_elastic->Geom(lev).periodicity());
191 }
192
193 m_elastic->SetModel(a_ddw_mf);
194 if (m_psi)
195 for (int i = 0; i <= a_b_mf.finest_level; i++)
196 m_elastic->SetPsi(i, *(*m_psi)[i]);
197
198
199 for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
200 {
201 amrex::Box domain(linop->Geom(lev).growPeriodicDomain(2));
202 domain.convert(amrex::IntVect::TheNodeVector());
203 const Set::Scalar* dx = linop->Geom(lev).CellSize();
204 const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
205 for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
206 {
207 amrex::Box bx = mfi.grownnodaltilebox() & domain;
208 amrex::Array4<const Set::Vector> const& u = a_u_mf[lev]->array(mfi);
209 amrex::Array4<const Set::Vector> const& b = a_b_mf[lev]->array(mfi);
210 amrex::Array4<const Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
211 amrex::Array4<Set::Scalar> const& rhs = a_rhs_mf[lev]->array(mfi);
212
213 // This is for if psi is not being used and has not been set
214 if (m_psi == nullptr)
215 {
216 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
217 {
218 auto sten = Numeric::GetStencil(i, j, k, bx);
219 // Do this if on the domain boundary
220 if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
221 {
222 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
223 Set::Vector ret = m_elastic->GetBC()(u(i, j, k), gradu, dw(i, j, k), i, j, k, bx);
224 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - ret(d);
225 }
226 else
227 {
228 Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx, sten);
229 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - divdw(d);
230 }
231 });
232 }
233 // This is EXACTLY THE SAME AS THE ABOVE, excpet that we are now accounting
234 // for psi terms. The reason we do this is because there's no good way (yet) to
235 // initialize patches if psi = nullptr
236 else
237 {
238 const amrex::Dim3 boxlo = amrex::lbound(bx), boxhi = amrex::ubound(bx);
239 amrex::Array4<const Set::Scalar> const& psi = (*m_psi)[lev]->array(mfi);
240 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
241 {
242 auto sten = Numeric::GetStencil(i, j, k, bx);
243 // Do this if on the domain boundary
244 if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
245 {
246 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
247 Set::Scalar psiavg = 1.0;
248 if (AMREX_D_TERM(i > boxlo.x && i<boxhi.x, && j>boxlo.y && j<boxhi.y, && k>boxlo.z && k < boxhi.z))
249 {
250 psiavg = Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0);
251 }
252 Set::Vector ret = m_elastic->GetBC()(u(i, j, k), gradu, dw(i, j, k) * psiavg, i, j, k, bx);
253 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - ret(d);
254 }
255 else
256 {
257 Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx, sten);
258 if (AMREX_D_TERM(i > boxlo.x && i<boxhi.x, && j>boxlo.y && j<boxhi.y, && k>boxlo.z && k < boxhi.z))
259 {
260 divdw *= Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0);
261 Set::Vector gradpsi = Numeric::CellGradientOnNode(psi, i, j, k, 0, dx);
262 divdw += dw(i, j, k) * gradpsi;
263 }
264 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - divdw(d);
265 }
266 });
267 }
268 }
269
270 a_rhs_mf[lev]->setMultiGhost(true);
271 a_rhs_mf[lev]->FillBoundaryAndSync(m_elastic->Geom(lev).periodicity());
272 a_ddw_mf[lev]->setMultiGhost(true);
273 a_ddw_mf[lev]->FillBoundaryAndSync(m_elastic->Geom(lev).periodicity());
274 }
275 }
276
277
278public:
279 Set::Scalar solve(const Set::Field<Set::Vector>& a_u_mf,
280 const Set::Field<Set::Vector>& a_b_mf,
281 Set::Field<T>& a_model_mf,
282 Real a_tol_rel, Real a_tol_abs, const char* checkpoint_file = nullptr)
283 {
284 Set::Field<Set::Scalar> dsol_mf, rhs_mf;
285 Set::Field<Set::Matrix> dw_mf;
286 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
287
288 dsol_mf.resize(a_u_mf.finest_level + 1); dsol_mf.finest_level = a_u_mf.finest_level;
289 dw_mf.resize(a_u_mf.finest_level + 1); dw_mf.finest_level = a_u_mf.finest_level;
290 ddw_mf.resize(a_u_mf.finest_level + 1); ddw_mf.finest_level = a_u_mf.finest_level;
291 rhs_mf.resize(a_u_mf.finest_level + 1); rhs_mf.finest_level = a_u_mf.finest_level;
292 for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
293 {
294 dsol_mf.Define(lev, a_u_mf[lev]->boxArray(),
295 a_u_mf[lev]->DistributionMap(),
296 a_u_mf.NComp(),
297 a_u_mf[lev]->nGrow());
298 dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
299 a_b_mf[lev]->DistributionMap(),
300 1,
301 a_b_mf[lev]->nGrow());
302 ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
303 a_b_mf[lev]->DistributionMap(),
304 1,
305 a_b_mf[lev]->nGrow());
306 rhs_mf.Define(lev, a_b_mf[lev]->boxArray(),
307 a_b_mf[lev]->DistributionMap(),
308 a_b_mf.NComp(),
309 a_b_mf[lev]->nGrow());
310
311 dsol_mf[lev]->setVal(0.0);
312 dw_mf[lev]->setVal(Set::Matrix::Zero());
313 ddw_mf[lev]->setVal(Set::Matrix4<AMREX_SPACEDIM, T::sym>::Zero());
314
315 a_b_mf.Copy(lev, *rhs_mf[lev], 0, 2);
316 //amrex::MultiFab::Copy(*rhs_mf[lev], *a_b_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
317 }
318
319 for (int nriter = 0; nriter < m_nriters; nriter++)
320 {
321 if (verbose > 0 && nriter < m_nriters) Util::Message(INFO, "Newton Iteration ", nriter + 1, " of ", m_nriters);
322
323 prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf, ddw_mf, a_model_mf);
324
325
326 if (nriter == m_nriters) break;
327
328 Solver::Nonlocal::Linear::solve(dsol_mf, rhs_mf, a_tol_rel, a_tol_abs, checkpoint_file);
329
330 Set::Scalar cornorm = 0, solnorm = 0;
331 for (int lev = 0; lev < dsol_mf.size(); ++lev)
332 {
333 for (int comp = 0; comp < AMREX_SPACEDIM; comp++)
334 {
335 Set::Scalar tmpcornorm = dsol_mf[lev]->norm0(comp, 0);
336 if (tmpcornorm > cornorm) cornorm = tmpcornorm;
337
338 //Set::Scalar tmpsolnorm = a_u_mf[lev]->norm0(comp,0);
339 //if (tmpsolnorm > solnorm) solnorm = tmpsolnorm;
340 }
341
342 }
343 Set::Scalar relnorm;
344 if (solnorm == 0) relnorm = cornorm;
345 else relnorm = cornorm / solnorm;
346 if (verbose > 0) Util::Message(INFO, "NR iteration ", nriter + 1, ", relative norm(ddisp) = ", relnorm);
347
348 for (int lev = 0; lev < dsol_mf.size(); ++lev)
349 a_u_mf.AddFrom(lev, *dsol_mf[lev], 0, 2);
350 //amrex::MultiFab::Add(*a_u_mf[lev], *dsol_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
351
352 if (relnorm < m_nrtolerance)
353 return relnorm;
354
355 }
356
357 return 0.0;
358 }
359
360 Set::Scalar solve(const Set::Field<Set::Scalar>& a_u_mf,
361 const Set::Field<Set::Scalar>& a_b_mf,
362 Set::Field<T>& a_model_mf,
363 Real a_tol_rel, Real a_tol_abs, const char* checkpoint_file = nullptr)
364 {
365 Set::Field<Set::Scalar> dsol_mf, rhs_mf;
366 Set::Field<Set::Matrix> dw_mf;
367 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
368
369 dsol_mf.resize(a_u_mf.finest_level + 1); dsol_mf.finest_level = a_u_mf.finest_level;
370 dw_mf.resize(a_u_mf.finest_level + 1); dw_mf.finest_level = a_u_mf.finest_level;
371 ddw_mf.resize(a_u_mf.finest_level + 1); ddw_mf.finest_level = a_u_mf.finest_level;
372 rhs_mf.resize(a_u_mf.finest_level + 1); rhs_mf.finest_level = a_u_mf.finest_level;
373 for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
374 {
375 dsol_mf.Define(lev, a_u_mf[lev]->boxArray(),
376 a_u_mf[lev]->DistributionMap(),
377 a_u_mf[lev]->nComp(),
378 a_u_mf[lev]->nGrow());
379 dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
380 a_b_mf[lev]->DistributionMap(),
381 1,
382 a_b_mf[lev]->nGrow());
383 ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
384 a_b_mf[lev]->DistributionMap(),
385 1,
386 a_b_mf[lev]->nGrow());
387 rhs_mf.Define(lev, a_b_mf[lev]->boxArray(),
388 a_b_mf[lev]->DistributionMap(),
389 a_b_mf[lev]->nComp(),
390 a_b_mf[lev]->nGrow());
391
392 dsol_mf[lev]->setVal(0.0);
393 dw_mf[lev]->setVal(Set::Matrix::Zero());
394 ddw_mf[lev]->setVal(Set::Matrix4<AMREX_SPACEDIM, T::sym>::Zero());
395
396 amrex::MultiFab::Copy(*rhs_mf[lev], *a_b_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
397 }
398
399 for (int nriter = 0; nriter < m_nriters; nriter++)
400 {
401 if (verbose > 0 && nriter < m_nriters) Util::Message(INFO, "Newton Iteration ", nriter + 1, " of ", m_nriters);
402
403 prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf, ddw_mf, a_model_mf);
404
405
406 if (nriter == m_nriters) break;
407
408 Solver::Nonlocal::Linear::solve(dsol_mf, rhs_mf, a_tol_rel, a_tol_abs, checkpoint_file);
409
410 Set::Scalar cornorm = 0, solnorm = 0;
411 for (int lev = 0; lev < dsol_mf.size(); ++lev)
412 {
413 for (int comp = 0; comp < AMREX_SPACEDIM; comp++)
414 {
415 Set::Scalar tmpcornorm = dsol_mf[lev]->norm0(comp, 0);
416 if (tmpcornorm > cornorm) cornorm = tmpcornorm;
417
418 Set::Scalar tmpsolnorm = a_u_mf[lev]->norm0(comp, 0);
419 if (tmpsolnorm > solnorm) solnorm = tmpsolnorm;
420 }
421
422 }
423 Set::Scalar relnorm;
424 if (solnorm == 0) relnorm = cornorm;
425 else relnorm = cornorm / solnorm;
426 if (verbose > 0) Util::Message(INFO, "NR iteration ", nriter + 1, ", relative norm(ddisp) = ", relnorm);
427
428 for (int lev = 0; lev < dsol_mf.size(); ++lev)
429 amrex::MultiFab::Add(*a_u_mf[lev], *dsol_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
430
431 if (relnorm < m_nrtolerance)
432 return relnorm;
433
434 }
435
436 return 0.0;
437 }
438 Set::Scalar solve(const Set::Field<Set::Scalar>& a_u_mf,
439 const Set::Field<Set::Scalar>& a_b_mf,
440 Set::Field<T>& a_model_mf)
441 {
442 return solve(a_u_mf, a_b_mf, a_model_mf, tol_rel, tol_abs);
443 }
444
445 void compResidual(Set::Field<Set::Scalar>& a_res_mf,
446 Set::Field<Set::Scalar>& a_u_mf,
447 Set::Field<Set::Scalar>& a_b_mf,
448 Set::Field<T>& a_model_mf)
449 {
450 Set::Field<Set::Matrix> dw_mf;
451 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
452 dw_mf.resize(a_u_mf.size());
453 ddw_mf.resize(a_u_mf.size());
454 for (int lev = 0; lev < a_u_mf.size(); lev++)
455 {
456 dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
457 a_b_mf[lev]->DistributionMap(),
458 1, a_b_mf[lev]->nGrow());
459 ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
460 a_b_mf[lev]->DistributionMap(),
461 1, a_b_mf[lev]->nGrow());
462 dw_mf[lev]->setVal(Set::Matrix::Zero());
463 }
464
465 prepareForSolve(a_u_mf, a_b_mf, a_res_mf, dw_mf, ddw_mf, a_model_mf);
466 }
467
468 void compResidual(Set::Field<Set::Vector>& a_res_mf,
469 Set::Field<Set::Vector>& a_u_mf,
470 Set::Field<Set::Vector>& a_b_mf,
471 Set::Field<T>& a_model_mf)
472 {
473 Set::Field<Set::Matrix> dw_mf;
474 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
475 Set::Field<Set::Scalar> res_mf(a_res_mf.size());
476 dw_mf.resize(a_u_mf.size());
477 ddw_mf.resize(a_u_mf.size());
478 res_mf.resize(a_u_mf.size());
479 for (int lev = 0; lev < a_u_mf.size(); lev++)
480 {
481 dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
482 a_b_mf[lev]->DistributionMap(),
483 1, a_b_mf[lev]->nGrow());
484 ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
485 a_b_mf[lev]->DistributionMap(),
486 1, a_b_mf[lev]->nGrow());
487 res_mf.Define(lev, a_b_mf[lev]->boxArray(),
488 a_b_mf[lev]->DistributionMap(),
489 AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
490 dw_mf[lev]->setVal(Set::Matrix::Zero());
491 }
492
493 prepareForSolve(a_u_mf, a_b_mf, res_mf, dw_mf, ddw_mf, a_model_mf);
494
495 for (int lev = 0; lev < a_res_mf.size(); ++lev)
496 {
497 Util::RealFillBoundary(*res_mf[lev], m_elastic->Geom(lev));
498 a_res_mf.CopyFrom(lev, *res_mf[lev], 0, 2);
499 }
500 }
501
502 void compLinearSolverResidual(Set::Field<Set::Vector>& a_res_mf,
503 Set::Field<Set::Vector>& a_u_mf,
504 Set::Field<Set::Vector>& a_b_mf)
505 {
506 Util::Assert(INFO, TEST(a_res_mf.finest_level == a_u_mf.finest_level));
507 Util::Assert(INFO, TEST(a_u_mf.finest_level == a_b_mf.finest_level));
508 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);
509 for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
510 {
511 res_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
512 sol_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
513 rhs_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
514
515 a_u_mf.Copy(lev, *sol_mf[lev], 0, 2);
516 a_b_mf.Copy(lev, *rhs_mf[lev], 0, 2);
517 }
518
519 mlmg->compResidual(amrex::GetVecOfPtrs(res_mf), amrex::GetVecOfPtrs(sol_mf), amrex::GetVecOfConstPtrs(rhs_mf));
520
521 for (int lev = 0; lev <= a_res_mf.finest_level; ++lev)
522 {
523 a_res_mf.CopyFrom(lev, *res_mf[lev], 0, 2);
524 }
525 }
526
527
528 void W(Set::Field<Set::Scalar>& a_w_mf,
529 Set::Field<Set::Scalar>& a_u_mf,
530 Set::Field<T>& a_model_mf)
531 {
532 for (int lev = 0; lev < a_u_mf.size(); lev++)
533 {
534 BL_PROFILE("Solver::Nonlocal::Newton::DW()");
535
536 const amrex::Real* DX = linop->Geom(lev).CellSize();
537 amrex::Box domain(linop->Geom(lev).Domain());
538 domain.convert(amrex::IntVect::TheNodeVector());
539
540 for (MFIter mfi(*a_u_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
541 {
542 const Box& bx = mfi.tilebox();
543 amrex::Array4<T> const& C = a_model_mf[lev]->array(mfi);
544 amrex::Array4<amrex::Real> const& w = a_w_mf[lev]->array(mfi);
545 amrex::Array4<const amrex::Real> const& u = a_u_mf[lev]->array(mfi);
546 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
547 {
548 Set::Matrix gradu;
549
550 auto sten = Numeric::GetStencil(i, j, k, bx);
551
552 // Fill gradu
553 for (int p = 0; p < AMREX_SPACEDIM; p++)
554 {
555 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
556 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
557 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
558 }
559
560 if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
561 w(i, j, k) = C(i, j, k).W(gradu);
562 else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
563 w(i, j, k) = C(i, j, k).W(0.5 * (gradu + gradu.transpose()));
564 else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
565 w(i, j, k) = C(i, j, k).W(gradu + Set::Matrix::Identity());
566 });
567 }
568 }
569 }
570
571 void DW(Set::Field<Set::Scalar>& a_dw_mf,
572 Set::Field<Set::Scalar>& a_u_mf,
573 Set::Field<T>& a_model_mf)
574 {
575 for (int lev = 0; lev < a_u_mf.size(); lev++)
576 {
577 BL_PROFILE("Solver::Nonlocal::Newton::DW()");
578
579 const amrex::Real* DX = linop->Geom(lev).CellSize();
580 amrex::Box domain(linop->Geom(lev).Domain());
581 domain.convert(amrex::IntVect::TheNodeVector());
582
583 for (MFIter mfi(*a_u_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
584 {
585 const Box& bx = mfi.tilebox();
586 amrex::Array4<T> const& C = a_model_mf[lev]->array(mfi);
587 amrex::Array4<amrex::Real> const& dw = a_dw_mf[lev]->array(mfi);
588 amrex::Array4<const amrex::Real> const& u = a_u_mf[lev]->array(mfi);
589 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
590 {
591 Set::Matrix gradu;
592
593 auto sten = Numeric::GetStencil(i, j, k, bx);
594
595 // Fill gradu
596 for (int p = 0; p < AMREX_SPACEDIM; p++)
597 {
598 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
599 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
600 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
601 }
602
603 Set::Matrix sig = Set::Matrix::Zero();
604
605 if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
606 sig = C(i, j, k).DW(gradu);
607 else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
608 sig = C(i, j, k).DW(0.5 * (gradu + gradu.transpose()));
609 else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
610 sig = C(i, j, k).DW(gradu + Set::Matrix::Identity());
611
612 // = C(i,j,k)(gradu,m_homogeneous);
613
614 AMREX_D_PICK(dw(i, j, k, 0) = sig(0, 0);
615 ,
616 dw(i, j, k, 0) = sig(0, 0); dw(i, j, k, 1) = sig(0, 1);
617 dw(i, j, k, 2) = sig(1, 0); dw(i, j, k, 3) = sig(1, 1);
618 ,
619 dw(i, j, k, 0) = sig(0, 0); dw(i, j, k, 1) = sig(0, 1); dw(i, j, k, 2) = sig(0, 2);
620 dw(i, j, k, 3) = sig(1, 0); dw(i, j, k, 4) = sig(1, 1); dw(i, j, k, 5) = sig(1, 2);
621 dw(i, j, k, 6) = sig(2, 0); dw(i, j, k, 7) = sig(2, 1); dw(i, j, k, 8) = sig(2, 2););
622
623 });
624 }
625 }
626 }
627
628
629public:
630 int m_nriters = 1;
631 Set::Scalar m_nrtolerance = 0.0;
632 Operator::Elastic<T::sym>* m_elastic;
633 //BC::Operator::Elastic::Elastic *m_bc;
634
635 Set::Field<Set::Scalar>* m_psi = nullptr;
636
637public:
638 // These paramters control a standard Newton-Raphson solve.
639 //
640 // **Note**:
641 // This class inherits all of the linear solve paramters
642 // from its parent class, :ref:`Solver::Nonlocal::Linear`
643 static void Parse(Newton<T>& value, amrex::ParmParse& pp)
644 {
645 Linear::Parse(value, pp);
646
647 // Number of newton-raphson iterations.
648 pp_query("nriters", value.m_nriters);
649
650 // Tolerance to use for newton-raphson convergence
651 pp_query("nrtolerance", value.m_nrtolerance);
652 }
653
654};
655} // namespace Nonlocal
656} // namespace Solver
657
658
659#endif
pp_query
#define pp_query(...)
Definition ParmParse.H:108
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:59
Operator::Elastic::GetBC
::BC::Operator::Elastic::Elastic & GetBC()
Definition Elastic.H:66
Operator::Elastic::SetPsi
void SetPsi(int amrlev, const amrex::MultiFab &a_psi)
Definition Elastic.cpp:120
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:46
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:438
Solver::Nonlocal::Newton::Newton
Newton(Operator::Elastic< T::sym > &a_op)
Definition Newton.H:20
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:360
Solver::Nonlocal::Newton::m_nrtolerance
Set::Scalar m_nrtolerance
Definition Newton.H:631
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:468
Solver::Nonlocal::Newton::Define
void Define(Operator::Elastic< T::sym > &a_op)
Definition Newton.H:30
Solver::Nonlocal::Newton::m_psi
Set::Field< Set::Scalar > * m_psi
Definition Newton.H:635
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:528
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:52
Solver::Nonlocal::Newton::Parse
static void Parse(Newton< T > &value, amrex::ParmParse &pp)
Definition Newton.H:643
Solver::Nonlocal::Newton::m_elastic
Operator::Elastic< T::sym > * m_elastic
Definition Newton.H:632
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:571
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:144
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:502
Solver::Nonlocal::Newton::setNRIters
void setNRIters(int a_nriters)
Definition Newton.H:44
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:445
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:279
Numeric::GetStencil
static 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:45
Numeric::CellGradientOnNode
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:699
Numeric::Gradient
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:681
Numeric::Divergence
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:590
Set::Scalar
amrex::Real Scalar
Definition Base.H:18
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:19
Set::Matrix
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, AMREX_SPACEDIM > Matrix
Definition Base.H:22
Solver
A bunch of solvers.
Definition CG.H:6
Util::Assert
AMREX_FORCE_INLINE void Assert(std::string file, std::string func, int line, std::string smt, bool pass, Args const &... args)
Definition Util.H:58
Util::Message
void Message(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:129
Util::RealFillBoundary
AMREX_FORCE_INLINE void RealFillBoundary(amrex::FabArray< amrex::BaseFab< T > > &a_mf, const amrex::Geometry &, const int nghost=2)
Definition Util.H:339
Numeric::Interpolate::CellToNodeAverage
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:1390
Generated by doxygen 1.9.8