LCOV - code coverage report
Current view: top level - src/Solver/Nonlocal - Newton.H (source / functions) Coverage Total Hit
Test: coverage_merged.info Lines: 98.0 % 151 148
Test Date: 2025-04-03 04:02:21 Functions: 32.9 % 216 71

            Line data    Source code
       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           25 :     Newton() {};
      19              : 
      20              :     Newton(Operator::Elastic<T::sym>& a_op)
      21              :     {
      22              :         this->Define(a_op);
      23              :     };
      24              : 
      25           25 :     ~Newton()
      26              :     {
      27           25 :         if (m_defined) Clear();
      28           25 :     }
      29              : 
      30          426 :     void Define(Operator::Elastic<T::sym>& a_op)
      31              :     {
      32          426 :         Linear::Define(a_op);
      33          426 :         m_elastic = dynamic_cast<Operator::Elastic<T::sym> *>(linop);
      34              :         //m_bc = &m_elastic->GetBC();
      35          426 :     }
      36          426 :     void Clear()
      37              :     {
      38          426 :         Linear::Clear();
      39          426 :         m_elastic = nullptr;
      40              :         //m_bc = nullptr;
      41          426 :     }
      42              : 
      43              : 
      44              :     void setNRIters(int a_nriters) { m_nriters = a_nriters; }
      45              : 
      46           14 :     void setPsi(Set::Field<Set::Scalar>& a_psi)
      47              :     {
      48           14 :         m_psi = &a_psi;
      49           14 :     }
      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         1326 :     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         2698 :         for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
     152              :         {
     153         1372 :             amrex::Box domain(linop->Geom(lev).Domain());
     154         1372 :             domain.convert(amrex::IntVect::TheNodeVector());
     155         1372 :             const Set::Scalar* dx = linop->Geom(lev).CellSize();
     156         1372 :             Set::Vector DX(linop->Geom(lev).CellSize());
     157         2877 :             for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
     158              :             {
     159         1505 :                 amrex::Box bx = mfi.grownnodaltilebox();
     160         1505 :                 bx = bx & domain;
     161              : 
     162         1505 :                 amrex::Array4<const T>           const& model = a_model_mf[lev]->array(mfi);
     163         1505 :                 amrex::Array4<const Set::Vector> const& u = a_u_mf[lev]->array(mfi);
     164         1505 :                 amrex::Array4<Set::Matrix>       const& dw = a_dw_mf[lev]->array(mfi);
     165         1505 :                 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       412439 :                 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
     169              :                 {
     170       209757 :                     auto sten = Numeric::GetStencil(i, j, k, bx);
     171              : 
     172       209757 :                     Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
     173       209757 :                     Set::Matrix kinvar;
     174       419514 :                     if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
     175       138723 :                         kinvar = gradu; // gradu
     176       142068 :                     else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
     177            0 :                         kinvar = 0.5 * (gradu + gradu.transpose()); // epsilon
     178       142068 :                     else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
     179        71034 :                         kinvar = gradu + Set::Matrix::Identity(); // F
     180              : 
     181       629271 :                     dw(i, j, k) = model(i, j, k).DW(kinvar);
     182       629271 :                     ddw(i, j, k) = model(i, j, k).DDW(kinvar);
     183              : 
     184              :                 });
     185              :             }
     186              : 
     187         1372 :             Util::RealFillBoundary(*a_dw_mf[lev], m_elastic->Geom(lev));
     188         1372 :             Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
     189              :         }
     190              : 
     191         1326 :         m_elastic->SetModel(a_ddw_mf);
     192         1326 :         if (m_psi)
     193           58 :             for (int i = 0; i <= a_b_mf.finest_level; i++)
     194           44 :                 m_elastic->SetPsi(i, *(*m_psi)[i]);
     195              : 
     196              : 
     197         2698 :         for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
     198              :         {
     199         1372 :             amrex::Box domain(linop->Geom(lev).Domain());
     200         1372 :             domain.convert(amrex::IntVect::TheNodeVector());
     201         1372 :             const Set::Scalar* dx = linop->Geom(lev).CellSize();
     202         1372 :             const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
     203         2877 :             for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
     204              :             {
     205         1505 :                 amrex::Box bx = mfi.grownnodaltilebox() & domain;
     206         1505 :                 amrex::Array4<const Set::Vector>  const& u = a_u_mf[lev]->array(mfi);
     207         1505 :                 amrex::Array4<const Set::Vector>  const& b = a_b_mf[lev]->array(mfi);
     208         1505 :                 amrex::Array4<const Set::Matrix>  const& dw = a_dw_mf[lev]->array(mfi);
     209         1505 :                 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         1505 :                 if (m_psi == nullptr)
     213              :                 {
     214       122666 :                     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
     215              :                     {
     216       121275 :                         auto sten = Numeric::GetStencil(i, j, k, bx);
     217              :                         // Do this if on the domain boundary
     218       121275 :                         if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
     219              :                         {
     220        46805 :                             Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
     221       140415 :                             Set::Vector ret = m_elastic->GetBC()(u(i, j, k), gradu, dw(i, j, k), i, j, k, bx);
     222       292825 :                             for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - ret(d);
     223        46805 :                         }
     224              :                         else
     225              :                         {
     226        74470 :                             Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx, sten);
     227       388550 :                             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          114 :                     const amrex::Dim3 boxlo = amrex::lbound(bx), boxhi = amrex::ubound(bx);
     237          114 :                     amrex::Array4<const Set::Scalar>  const& psi = (*m_psi)[lev]->array(mfi);
     238        88596 :                     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
     239              :                     {
     240        88482 :                         auto sten = Numeric::GetStencil(i, j, k, bx);
     241              :                         // Do this if on the domain boundary
     242        88482 :                         if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
     243              :                         {
     244         2560 :                             Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
     245         2560 :                             Set::Scalar psiavg = 1.0;
     246         2560 :                             if (AMREX_D_TERM(i > boxlo.x && i<boxhi.x, && j>boxlo.y && j<boxhi.y, && k>boxlo.z && k < boxhi.z))
     247              :                             {
     248            0 :                                 psiavg = Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0);
     249              :                             }
     250         7680 :                             Set::Vector ret = m_elastic->GetBC()(u(i, j, k), gradu, dw(i, j, k) * psiavg, i, j, k, bx);
     251        12800 :                             for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - ret(d);
     252         2560 :                         }
     253              :                         else
     254              :                         {
     255        85922 :                             Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx, sten);
     256        85922 :                             if (AMREX_D_TERM(i > boxlo.x && i<boxhi.x, && j>boxlo.y && j<boxhi.y, && k>boxlo.z && k < boxhi.z))
     257              :                             {
     258       153860 :                                 divdw *= Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0);
     259        76930 :                                 Set::Vector gradpsi = Numeric::CellGradientOnNode(psi, i, j, k, 0, dx);
     260       153860 :                                 divdw += dw(i, j, k) * gradpsi;
     261              :                             }
     262       429610 :                             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         1372 :             Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
     269         1372 :             Util::RealFillBoundary(*a_rhs_mf[lev], m_elastic->Geom(lev));
     270              :         }
     271         1326 :     }
     272              : 
     273              : 
     274              : public:
     275          426 :     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          426 :         Set::Field<Set::Scalar> dsol_mf, rhs_mf;
     281          426 :         Set::Field<Set::Matrix> dw_mf;
     282          426 :         Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
     283              : 
     284          426 :         dsol_mf.resize(a_u_mf.finest_level + 1); dsol_mf.finest_level = a_u_mf.finest_level;
     285          426 :         dw_mf.resize(a_u_mf.finest_level + 1); dw_mf.finest_level = a_u_mf.finest_level;
     286          426 :         ddw_mf.resize(a_u_mf.finest_level + 1); ddw_mf.finest_level = a_u_mf.finest_level;
     287          426 :         rhs_mf.resize(a_u_mf.finest_level + 1); rhs_mf.finest_level = a_u_mf.finest_level;
     288          898 :         for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
     289              :         {
     290          944 :             dsol_mf.Define(lev, a_u_mf[lev]->boxArray(),
     291          472 :                 a_u_mf[lev]->DistributionMap(),
     292              :                 a_u_mf.NComp(),
     293          472 :                 a_u_mf[lev]->nGrow());
     294          944 :             dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
     295          472 :                 a_b_mf[lev]->DistributionMap(),
     296              :                 1,
     297          472 :                 a_b_mf[lev]->nGrow());
     298          944 :             ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
     299          472 :                 a_b_mf[lev]->DistributionMap(),
     300              :                 1,
     301          472 :                 a_b_mf[lev]->nGrow());
     302          944 :             rhs_mf.Define(lev, a_b_mf[lev]->boxArray(),
     303          472 :                 a_b_mf[lev]->DistributionMap(),
     304              :                 a_b_mf.NComp(),
     305          472 :                 a_b_mf[lev]->nGrow());
     306              : 
     307          472 :             dsol_mf[lev]->setVal(0.0);
     308          472 :             dw_mf[lev]->setVal(Set::Matrix::Zero());
     309          472 :             ddw_mf[lev]->setVal(Set::Matrix4<AMREX_SPACEDIM, T::sym>::Zero());
     310              : 
     311          472 :             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         1750 :         for (int nriter = 0; nriter < m_nriters; nriter++)
     316              :         {
     317         6630 :             if (verbose > 0 && nriter < m_nriters) Util::Message(INFO, "Newton Iteration ", nriter + 1, " of ", m_nriters);
     318              : 
     319         1326 :             prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf, ddw_mf, a_model_mf);
     320              : 
     321              : 
     322         1326 :             if (nriter == m_nriters) break;
     323              : 
     324         1326 :             Solver::Nonlocal::Linear::solve(dsol_mf, rhs_mf, a_tol_rel, a_tol_abs, checkpoint_file);
     325              : 
     326         1326 :             Set::Scalar cornorm = 0, solnorm = 0;
     327         2698 :             for (int lev = 0; lev < dsol_mf.size(); ++lev)
     328              :             {
     329         4416 :                 for (int comp = 0; comp < AMREX_SPACEDIM; comp++)
     330              :                 {
     331         3044 :                     Set::Scalar tmpcornorm = dsol_mf[lev]->norm0(comp, 0);
     332         3044 :                     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         1326 :             if (solnorm == 0) relnorm = cornorm;
     341            0 :             else relnorm = cornorm / solnorm;
     342         6630 :             if (verbose > 0) Util::Message(INFO, "NR iteration ", nriter + 1, ", relative norm(ddisp) = ", relnorm);
     343              : 
     344         2698 :             for (int lev = 0; lev < dsol_mf.size(); ++lev)
     345         1372 :                 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         1326 :             if (relnorm < m_nrtolerance)
     349            2 :                 return relnorm;
     350              : 
     351              :         }
     352              : 
     353          424 :         return 0.0;
     354          426 :     }
     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            4 :     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           28 :         Util::Assert(INFO, TEST(a_res_mf.finest_level == a_u_mf.finest_level));
     504           28 :         Util::Assert(INFO, TEST(a_u_mf.finest_level == a_b_mf.finest_level));
     505            4 :         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           18 :         for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
     507              :         {
     508           14 :             res_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
     509           14 :             sol_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
     510           14 :             rhs_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
     511              : 
     512           14 :             a_u_mf.Copy(lev, *sol_mf[lev], 0, 2);
     513           14 :             a_b_mf.Copy(lev, *rhs_mf[lev], 0, 2);
     514              :         }
     515              : 
     516            4 :         mlmg->compResidual(amrex::GetVecOfPtrs(res_mf), amrex::GetVecOfPtrs(sol_mf), amrex::GetVecOfConstPtrs(rhs_mf));
     517              : 
     518           18 :         for (int lev = 0; lev <= a_res_mf.finest_level; ++lev)
     519              :         {
     520           14 :             a_res_mf.CopyFrom(lev, *res_mf[lev], 0, 2);
     521              :         }
     522            4 :     }
     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           19 :     static void Parse(Newton<T>& value, amrex::ParmParse& pp)
     641              :     {
     642           19 :         Linear::Parse(value, pp);
     643              : 
     644              :         // Number of newton-raphson iterations.
     645           19 :         pp_query("nriters", value.m_nriters);
     646              : 
     647              :         // Tolerance to use for newton-raphson convergence
     648           19 :         pp_query("nrtolerance", value.m_nrtolerance); 
     649           19 :     }
     650              : 
     651              : };
     652              : } // namespace Nonlocal
     653              : } // namespace Solver
     654              : 
     655              : 
     656              : #endif
        

Generated by: LCOV version 2.0-1