Line data    Source code

       1              : #include <AMReX_MLPoisson.H>
       2              : 
       3              : #ifdef ALAMO_FFT
       4              : #include <AMReX_FFT.H>
       5              : #endif
       6              : 
       7              : 
       8              : #include "PFC.H"
       9              : #include "BC/Constant.H"
      10              : #include "IO/ParmParse.H"
      11              : #include "IC/Random.H"
      12              : #include "Set/Set.H"
      13              : 
      14              : namespace Integrator
      15              : {
      16            0 : PFC::PFC() : Integrator()
      17              : {
      18            0 : }
      19            0 : PFC::~PFC()
      20              : {
      21            0 :     delete ic;
      22            0 :     delete bc;
      23            0 : }
      24              : 
      25            0 : void PFC::Parse(PFC &value, IO::ParmParse &pp)
      26              : {
      27              :     // frequency term
      28            0 :     pp.query_required("q0",value.q0);
      29              :     // chemical potential width
      30            0 :     pp.query_required("eps",    value.eps);
      31              : 
      32              :     // initial condition for :math:`\eta`
      33            0 :     pp.select_default<IC::Random>("eta.ic", value.ic, value.geom);
      34              :     // boundary condition for :math:`\eta`
      35            0 :     pp.select_default<BC::Constant>("eta.bc", value.bc, 1);
      36              : 
      37            0 :     value.RegisterNewFab(value.eta_mf, value.bc, 1, 1, "eta",true);
      38            0 :     value.RegisterNewFab(value.grad_chempot_mf, value.bc, 1, 1, "grad_chempot",true);
      39            0 : }
      40              : 
      41              : 
      42              : #ifdef ALAMO_FFT
      43              : void
      44              : PFC::Advance (int lev, Set::Scalar /*time*/, Set::Scalar dt)
      45              : {
      46              :     //
      47              :     // FFT Boilerplate
      48              :     //
      49              :     amrex::FFT::R2C my_fft(this->geom[lev].Domain());
      50              :     auto const &[cba, cdm] = my_fft.getSpectralDataLayout();
      51              :     const Set::Scalar* DX = geom[lev].CellSize();
      52              :     amrex::Box const & domain = this->geom[lev].Domain();
      53              :     Set::Scalar
      54              :         AMREX_D_DECL(
      55              :             pi_Lx = 2.0 * Set::Constant::Pi / geom[lev].Domain().length(0) / DX[0],
      56              :             pi_Ly = 2.0 * Set::Constant::Pi / geom[lev].Domain().length(1) / DX[1],
      57              :             pi_Lz = 2.0 * Set::Constant::Pi / geom[lev].Domain().length(2) / DX[2]);
      58              :     Set::Scalar scaling = 1.0 / geom[lev].Domain().d_numPts();
      59              :     
      60              :     //
      61              :     // Compute the gradient of the chemical potential in realspace
      62              :     //
      63              :     for ( amrex::MFIter mfi(*eta_mf[lev],true); mfi.isValid(); ++mfi )
      64              :     {
      65              :         const amrex::Box& bx = mfi.tilebox();
      66              :         amrex::Array4<const amrex::Real> const& eta = eta_mf[lev]->array(mfi);
      67              :         amrex::Array4<amrex::Real> const& grad_chempot    = grad_chempot_mf[lev]->array(mfi);
      68              :         amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
      69              :         {
      70              :             grad_chempot(i, j, k) = eta(i, j, k) * eta(i, j, k) * eta(i, j, k);
      71              :         });
      72              :     }
      73              : 
      74              :     grad_chempot_mf[lev]->FillBoundary();
      75              : 
      76              :     //
      77              :     // FFT of eta
      78              :     // 
      79              :     amrex::FabArray<amrex::BaseFab<amrex::GpuComplex<Set::Scalar> > > eta_hat_mf(cba, cdm, 1, 0);
      80              :     my_fft.forward(*eta_mf[lev], eta_hat_mf);
      81              : 
      82              :     //
      83              :     // FFT of chemical potential gradient
      84              :     //
      85              :     amrex::FabArray<amrex::BaseFab<amrex::GpuComplex<Set::Scalar> > > chempot_hat_mf(cba, cdm, 1, 0);
      86              :     my_fft.forward(*grad_chempot_mf[lev], chempot_hat_mf);
      87              : 
      88              :     //
      89              :     // Perform update in spectral coordinatees
      90              :     //
      91              :     for (amrex::MFIter mfi(eta_hat_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
      92              :     {
      93              :         const amrex::Box &bx = mfi.tilebox();
      94              : 
      95              :         amrex::Array4<amrex::GpuComplex<Set::Scalar>> const & eta_hat   =  eta_hat_mf.array(mfi);
      96              :         amrex::Array4<amrex::GpuComplex<Set::Scalar>> const & N_hat     =  chempot_hat_mf.array(mfi);
      97              : 
      98              :         amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int m, int n, int p) {
      99              : 
     100              :             // Get spectral coordinates
     101              :             AMREX_D_TERM(
     102              :                 Set::Scalar k1 = m * pi_Lx;,
     103              :                 Set::Scalar k2 = (n < domain.length(1)/2 ? n * pi_Ly : (n - domain.length(1)) * pi_Ly);,
     104              :                 Set::Scalar k3 = (p < domain.length(2)/2 ? p * pi_Lz : (p - domain.length(2)) * pi_Lz););
     105              : 
     106              :             Set::Scalar omega2 = AMREX_D_TERM(k1 * k1, + k2 * k2, + k3 * k3);
     107              :             Set::Scalar omega4 = omega2*omega2;
     108              :             Set::Scalar omega6 = omega2*omega2*omega2;
     109              : 
     110              :             eta_hat(m, n, p) = eta_hat(m, n, p) - dt * omega2 * N_hat(m, n, p);
     111              :             eta_hat(m,n,p) /= 1.0 + dt * ((q0*q0*q0*q0 - eps)*omega2  - 2.0* q0*q0 * omega4 + omega6);
     112              :             eta_hat(m,n,p) *= scaling;
     113              :         });
     114              :     }
     115              : 
     116              :     //
     117              :     // Transform solution back to realspace
     118              :     //
     119              :     my_fft.backward(eta_hat_mf, *eta_mf[lev]);
     120              : }
     121              : #else
     122              : void
     123            0 : PFC::Advance (int, Set::Scalar, Set::Scalar)
     124              : {   
     125            0 :     Util::Abort(INFO,"Alamo must be compiled with fft");
     126            0 : }
     127              : #endif 
     128              : 
     129              : void
     130            0 : PFC::Initialize (int lev)
     131              : {
     132            0 :     ic->Initialize(lev,eta_mf);
     133            0 :     grad_chempot_mf[lev]->setVal(0.0);
     134            0 : }
     135              : 
     136              : 
     137              : void
     138            0 : PFC::TagCellsForRefinement (int /*lev*/, amrex::TagBoxArray& /*a_tags*/, Set::Scalar /*time*/, int /*ngrow*/)
     139              : {
     140            0 : }
     141              : 
     142              : 
     143              : }