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