Line data    Source code

       1              : #include <AMReX_MLPoisson.H>
       2              : 
       3              : #include "CahnHilliard.H"
       4              : #include "BC/Nothing.H"
       5              : #include "BC/Constant.H"
       6              : #include "IO/ParmParse.H"
       7              : #include "Numeric/Stencil.H"
       8              : #include "Set/Set.H"
       9              : 
      10              : namespace Integrator
      11              : {
      12            1 : CahnHilliard::CahnHilliard() : Integrator()
      13              : {
      14            1 : }
      15            0 : CahnHilliard::~CahnHilliard()
      16              : {
      17            0 :     delete ic;
      18            0 :     delete bc;
      19            0 : }
      20              : 
      21            1 : void CahnHilliard::Parse(CahnHilliard &value, IO::ParmParse &pp)
      22              : {
      23              :     // Interface energy
      24            6 :     pp.query_default("gamma",value.gamma, 0.0005);
      25              :     // Mobility
      26            6 :     pp.query_default("L",    value.L,     1.0);
      27              :     // Regridding criterion
      28            5 :     pp.query_default("refinement_threshold",value.refinement_threshold, 1E100);
      29              : 
      30              :     // initial condition for :math:`\eta`
      31            2 :     pp.select_default<IC::Random>("eta.ic", value.ic, value.geom);
      32              :     // boundary condition for :math:`\eta`
      33            2 :     pp.select_default<BC::Constant>("eta.bc", value.bc, 1);
      34              : 
      35            3 :     value.RegisterNewFab(value.etanew_mf, value.bc, 1, 1, "eta",true);
      36            3 :     value.RegisterNewFab(value.etaold_mf, value.bc, 1, 1, "eta_old",false);
      37            3 :     value.RegisterNewFab(value.intermediate, value.bc, 1, 1, "int",false);
      38            1 : }
      39              : 
      40              : 
      41              : void
      42        13881 : CahnHilliard::Advance (int lev, Set::Scalar /*time*/, Set::Scalar dt)
      43              : {
      44        13881 :     std::swap(etaold_mf[lev], etanew_mf[lev]);
      45        13881 :     const amrex::Real* DX = geom[lev].CellSize();
      46        27762 :     for ( amrex::MFIter mfi(*etanew_mf[lev],true); mfi.isValid(); ++mfi )
      47              :     {
      48        13881 :         const amrex::Box& bx = mfi.tilebox();
      49        13881 :         amrex::Array4<const amrex::Real> const& eta = etaold_mf[lev]->array(mfi);
      50        13881 :         amrex::Array4<amrex::Real> const& inter    = intermediate[lev]->array(mfi);
      51        13881 :         amrex::Array4<amrex::Real> const& etanew    = etanew_mf[lev]->array(mfi);
      52              : 
      53        13881 :         amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
      54              :         {
      55    170569728 :             Set::Scalar lap_eta = Numeric::Laplacian(eta,i,j,k,0,DX);
      56              :             
      57              : 
      58    341139456 :             inter(i,j,k) =
      59    511709184 :                 eta(i,j,k)*eta(i,j,k)*eta(i,j,k)
      60    170569728 :                 - eta(i,j,k)
      61    170569728 :                 - gamma*lap_eta;
      62              : 
      63              : 
      64    511709184 :             etanew(i,j,k) = eta(i,j,k) - dt*inter(i,j,k); // Allen Cahn
      65    170569728 :         });
      66              : 
      67        13881 :         amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k){
      68    170569728 :             Set::Scalar lap_inter = Numeric::Laplacian(inter,i,j,k,0,DX);
      69              : 
      70    341139456 :             etanew(i,j,k) = eta(i,j,k) + dt*lap_inter;
      71    170569728 :         });
      72        13881 :     }
      73        13881 : }
      74              : 
      75              : void
      76            2 : CahnHilliard::Initialize (int lev)
      77              : {
      78            2 :     intermediate[lev]->setVal(0.0);
      79            2 :     ic->Initialize(lev,etanew_mf);
      80            2 :     ic->Initialize(lev,etaold_mf);
      81            2 : }
      82              : 
      83              : 
      84              : void
      85         2316 : CahnHilliard::TagCellsForRefinement (int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/)
      86              : {
      87         2316 :         const Set::Scalar* DX = geom[lev].CellSize();
      88         2316 :         Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
      89              : 
      90         4632 :         for (amrex::MFIter mfi(*etanew_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
      91              :         {
      92         2316 :             const amrex::Box& bx = mfi.tilebox();
      93         2316 :             amrex::Array4<char> const&     tags = a_tags.array(mfi);
      94         2316 :             Set::Patch<const Set::Scalar>   eta = (*etanew_mf[lev]).array(mfi);
      95              : 
      96         2316 :             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
      97              :             {
      98      9486336 :                 Set::Vector grad = Numeric::Gradient(eta, i, j, k, 0, DX);
      99      9486336 :                 if (grad.lpNorm<2>() * dr > refinement_threshold)
     100     18953182 :                     tags(i, j, k) = amrex::TagBox::SET;
     101      9486336 :             });
     102         2316 :         }
     103         2316 : }
     104              : 
     105              : 
     106              : }