From df2f03348695433323c03c6e66eca434c4005f0e Mon Sep 17 00:00:00 2001 From: Sreehari Perumanath Date: Wed, 13 Mar 2024 08:43:01 +0000 Subject: [PATCH] Laplace Working --- exec/multispec/AdvanceTimestepBousq.cpp | 167 ++++- exec/multispec/AdvanceTimestepBousq.cpp~ | 805 ++++++++++++++++++++++ exec/multispec/Checkpoint.cpp | 14 + exec/multispec/GNUmakefile | 4 +- exec/multispec/Init.cpp | 17 +- exec/multispec/WritePlotFile.cpp | 26 +- exec/multispec/main_driver.cpp | 28 +- exec/multispec/multispec_test_functions.H | 15 +- src_multispec/InitialProjection.cpp | 5 +- src_multispec/multispec_functions.cpp | 6 + src_multispec/multispec_namespace.H | 2 + 11 files changed, 1058 insertions(+), 31 deletions(-) create mode 100644 exec/multispec/AdvanceTimestepBousq.cpp~ diff --git a/exec/multispec/AdvanceTimestepBousq.cpp b/exec/multispec/AdvanceTimestepBousq.cpp index f3f381b67..69e91ee19 100644 --- a/exec/multispec/AdvanceTimestepBousq.cpp +++ b/exec/multispec/AdvanceTimestepBousq.cpp @@ -1,4 +1,4 @@ -//SCRTP + #include "hydro_functions.H" #include "common_functions.H" #include "gmres_functions.H" @@ -34,6 +34,10 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, MultiFab& permittivity, StochMassFlux& sMassFlux, StochMomFlux& sMomFlux, + MultiFab& phi_old, + MultiFab& phi_new, + MultiFab& phitot_old, + MultiFab& phitot_new, const Real& dt, const Real& time, const int& istep, @@ -123,7 +127,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, } } - if (use_multiphase || use_flory_huggins) { + if (use_multiphase || use_flory_huggins || use_ice_nucleation) { for (int d=0; d& umac, if (any_grav) { Abort("AdvanceTimestepBousq.cpp gravity not implemented"); } - + +if (use_ice_nucleation ==0) { if (variance_coef_mass != 0.) { sMassFlux.fillMassStochastic(); } @@ -218,7 +223,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, ComputeMassFluxdiv(rho_old,rhotot_old,Temp,diff_mass_fluxdiv,stoch_mass_fluxdiv, diff_mass_flux,stoch_mass_flux,sMassFlux,0.5*dt,time,geom,weights_mass, charge_old,grad_Epot_old,Epot,permittivity); - + // here is a reasonable place to call something to compute in reversible stress term // in this case want to get divergence so it looks like a add to rhs for stokes solver if (use_multiphase) { @@ -241,6 +246,18 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, } } + } else { + + ComputeRhotot(phi_old,phitot_old); + // compute reversible stress tensor ---added term + ComputeDivFHReversibleStress(div_reversible_stress,phitot_old,phi_old,geom); + + // add divergence of reversible stress to gmres_rhs_v + for (int d=0; d& umac, eta_ed[d].mult(2.,0,1,0); } +if (use_ice_nucleation ==0){ ////////////////////////////////////////////// // Step 2: compute reactions and mass fluxes at t^n ////////////////////////////////////////////// @@ -369,7 +387,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, end if */ - + if (advection_type == 1 || advection_type == 2) { // bds advection @@ -416,7 +434,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, Abort("Invalid advection_type"); } - + ////////////////////////////////////////////// /// Step 3: density prediction to t^{n+1/2} ////////////////////////////////////////////// @@ -441,7 +459,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, end if */ } - + // compute rhotot^{n+1/2} from rho^{n+1/2} in VALID REGION ComputeRhotot(rho_new,rhotot_new); @@ -451,6 +469,17 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, // average rho_i^{n+1/2} to faces AverageCCToFace(rho_new,rho_fc,0,nspecies,SPEC_BC_COMP,geom); +} else { + // compute rhotot^{n+1/2} from rho^{n+1/2} in VALID REGION + ComputeRhotot(rho_old,rhotot_new); + + // fill rho and rhotot ghost cells at t^{n+1/2} + FillRhoRhototGhost(rho_old,rhotot_new,geom); + + // average rho_i^{n+1/2} to faces + //AverageCCToFace(rho_old,rho_fc,0,nspecies,SPEC_BC_COMP,geom); +} + if (use_charged_fluid) { // compute total charge at t^{n+1/2} DotWithZ(rho_new,charge_new); @@ -462,7 +491,11 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, } // compute (eta,kappa)^{n+1/2} +if (use_ice_nucleation ==0){ ComputeEta(rho_new, rhotot_new, eta); +} else { + ComputeEta(phi_old, phitot_old, eta); +} if (AMREX_SPACEDIM == 2) { AverageCCToNode(eta,eta_ed[0],0,1,SPEC_BC_COMP,geom); } @@ -470,6 +503,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, AverageCCToEdge(eta,eta_ed,0,1,SPEC_BC_COMP,geom); } +if (use_ice_nucleation ==0){ ////////////////////////////////////////////// // Step 4: compute mass fluxes and reactions at t^{n+1/2} ////////////////////////////////////////////// @@ -526,7 +560,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, } else if (use_flory_huggins ==1){ // compute reversible stress tensor ---added term - ComputeDivFHReversibleStress(div_reversible_stress,rhotot_new,rho_new,geom); + ComputeDivFHReversibleStress(div_reversible_stress,rhotot_new,rho_new,geom); } @@ -569,7 +603,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, MkAdvSFluxdiv(umac_old,rho_fc,adv_mass_fluxdiv,geom,0,nspecies,false); MkAdvSFluxdiv(umac ,rho_fc,adv_mass_fluxdiv,geom,0,nspecies,true); } - + ////////////////////////////////////////////// // Step 5: density integration to t^{n+1} ////////////////////////////////////////////// @@ -607,6 +641,111 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, // average rho^{n+1} to faces AverageCCToFace(rhotot_new,rhotot_fc_new,0,1,RHO_BC_COMP,geom); + +} else { + AverageCCToFace(rhotot_old,rhotot_fc_new,0,1,RHO_BC_COMP,geom); + + amrex::Real EnScale = rhobar[0]*k_B*T_init[0]*fh_chi(1,0)/monomer_mass; // 2.79e09 + amrex::Real GradEnCoef = rhobar[0]*k_B*T_init[0]*fh_kappa(1,0)/monomer_mass; // 2.73e-6 + amrex::Real PotWellDepr = -2.687e-05*T_init[0]*T_init[0]+0.009943*T_init[0]-0.7128; // For polynomial potential + amrex::Real Mobility = 4.997e6*exp(-4438.0/T_init[0]); // AC Mobility from MATLAB simulations + + MultiFab phi_prd(ba,dmap,1,2); + + FillRhoRhototGhost(phi_old,phitot_old,geom); + + // predict phi at t+1/2 + for (MFIter mfi(phi_prd); mfi.isValid(); ++mfi ) + { + Box bx = mfi.tilebox(); + const Array4& phiOld = phi_old.array(mfi); + const Array4& phiPrd = phi_prd.array(mfi); + + const amrex::Array4& u = umac[0].array(mfi); + const amrex::Array4& v = umac[1].array(mfi); + const amrex::Array4& u_old = umac_old[0].array(mfi); + const amrex::Array4& v_old = umac_old[1].array(mfi); +#if (AMREX_SPACEDIM == 3) + const amrex::Array4& w = umac[2].array(mfi); + const amrex::Array4& w_old = umac_old[2].array(mfi); +#endif + int n=0; + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + + phiPrd(i,j,k) = phiOld(i,j,k,n) + + dt/2 * Mobility * 2 * GradEnCoef * + ( (phiOld(i+1,j,k,n) - 2.*phiOld(i,j,k,n) + phiOld(i-1,j,k,n)) / (dx[0]*dx[0]) + +(phiOld(i,j+1,k,n) - 2.*phiOld(i,j,k,n) + phiOld(i,j-1,k,n)) / (dx[1]*dx[1]) + ) + - dt/2 * ((u(i+1,j,k)+u_old(i+1,j,k))/2*(phiOld(i+1,j,k,n)+phiOld(i,j,k,n))/2-(u(i,j,k)+u_old(i,j,k))/2*(phiOld(i,j,k,n)+phiOld(i-1,j,k,n))/2)/dx[0] + - dt/2 * ((v(i,j+1,k)+v_old(i,j+1,k))/2*(phiOld(i,j+1,k,n)+phiOld(i,j,k,n))/2-(v(i,j,k)+v_old(i,j,k))/2*(phiOld(i,j,k,n)+phiOld(i,j-1,k,n))/2)/dx[1] +#if (AMREX_SPACEDIM == 3) + + dt/2 * Mobility * 2 * GradEnCoef * (phiOld(i,j,k+1,n) - 2.*phiOld(i,j,k,n) + phiOld(i,j,k-1,n)) / (dx[2]*dx[2]) + - dt/2 * ((w(i,j,k+1)+w_old(i,j,k+1))/2*(phiOld(i,j,k+1,n)+phiOld(i,j,k,n))/2-(w(i,j,k)+w_old(i,j,k))/2*(phiOld(i,j,k,n)+phiOld(i,j,k-1,n))/2)/dx[2] + + variance_coef_mass*sqrt(2.0*k_B*T_init[0]*Mobility*dt/(dx[0]*dx[1]*dx[2]))*amrex::RandomNormal(0,1)/sqrt(2) +#elif (AMREX_SPACEDIM == 2) + + variance_coef_mass*sqrt(2.0*k_B*T_init[0]*Mobility*dt/(dx[0]*dx[1]*cell_depth))*amrex::RandomNormal(0,1)/sqrt(2) +#endif + - dt/2 * Mobility * EnScale*(2*phiOld(i,j,k,n)*(phiOld(i,j,k,n)-1)*(phiOld(i,j,k,n)-1)+2*phiOld(i,j,k,n)*phiOld(i,j,k,n)*(phiOld(i,j,k,n)-1) - PotWellDepr) + ; + }); + } + + phi_prd.FillBoundary(geom.periodicity()); + // advance phi to t+1 + for (MFIter mfi(phi_new); mfi.isValid(); ++mfi ) + { + Box bx = mfi.tilebox(); + const Array4& phiNew = phi_new.array(mfi); + const Array4& phiOld = phi_old.array(mfi); + const Array4& phiPrd = phi_prd.array(mfi); + + const amrex::Array4& u = umac[0].array(mfi); + const amrex::Array4& v = umac[1].array(mfi); + const amrex::Array4& u_old = umac_old[0].array(mfi); + const amrex::Array4& v_old = umac_old[1].array(mfi); +#if (AMREX_SPACEDIM == 3) + const amrex::Array4& w = umac[2].array(mfi); + const amrex::Array4& w_old = umac_old[2].array(mfi); +#endif + int n=0; + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + + phiNew(i,j,k,n) = phiOld(i,j,k,n) + + dt * Mobility * 2 * GradEnCoef * + ( (phiPrd(i+1,j,k) - 2.*phiPrd(i,j,k) + phiPrd(i-1,j,k)) / (dx[0]*dx[0]) + +(phiPrd(i,j+1,k) - 2.*phiPrd(i,j,k) + phiPrd(i,j-1,k)) / (dx[1]*dx[1]) + ) + - dt * ((u(i+1,j,k)+u_old(i+1,j,k))/2*(phiPrd(i+1,j,k)+phiPrd(i,j,k))/2-(u(i,j,k)+u_old(i,j,k))/2*(phiPrd(i,j,k)+phiPrd(i-1,j,k))/2)/dx[0] + - dt * ((v(i,j+1,k)+v_old(i,j+1,k))/2*(phiPrd(i,j+1,k)+phiPrd(i,j,k))/2-(v(i,j,k)+v_old(i,j,k))/2*(phiPrd(i,j,k)+phiPrd(i,j-1,k))/2)/dx[1] +#if (AMREX_SPACEDIM == 3) + + dt * Mobility * 2 * GradEnCoef * (phiPrd(i,j,k+1) - 2.*phiPrd(i,j,k) + phiPrd(i,j,k-1)) / (dx[2]*dx[2]) + - dt * ((w(i,j,k+1)+w_old(i,j,k+1))/2*(phiPrd(i,j,k+1)+phiPrd(i,j,k))/2-(w(i,j,k)+w_old(i,j,k))/2*(phiPrd(i,j,k)+phiPrd(i,j,k-1))/2)/dx[2] + + variance_coef_mass*sqrt(2.0*k_B*T_init[0]*Mobility*dt/(dx[0]*dx[1]*dx[2]))*amrex::RandomNormal(0,1) +#elif (AMREX_SPACEDIM == 2) + + variance_coef_mass*sqrt(2.0*k_B*T_init[0]*Mobility*dt/(dx[0]*dx[1]*cell_depth))*amrex::RandomNormal(0,1) +#endif + - dt * Mobility * EnScale*(2*phiPrd(i,j,k)*(phiPrd(i,j,k)-1)*(phiPrd(i,j,k)-1)+2*phiPrd(i,j,k)*phiPrd(i,j,k)*(phiPrd(i,j,k)-1) - PotWellDepr) + ; + + phiNew(i,j,k,n+1) = 1-phiNew(i,j,k,n); + }); + } + + ComputeRhotot(phi_new,phitot_new); + // compute reversible stress tensor ---added term + ComputeDivFHReversibleStress(div_reversible_stress,phitot_new,phi_new,geom); + + MultiFab::Copy(phi_old,phi_new,0,0,nspecies,0); + + amrex::Print() << "Advanced phi (phi1_avg = " << phi_new.sum(0)/(n_cells[0]*n_cells[1] +#if (AMREX_SPACEDIM == 3) + *n_cells[2] +#endif + ) <<")\n"; +} if (use_charged_fluid) { // compute total charge at t^{n+1} @@ -618,14 +757,18 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, } } +if (use_ice_nucleation ==0){ ComputeEta(rho_new, rhotot_new, eta); +} else { + ComputeEta(phi_new, phitot_new, eta); +} if (AMREX_SPACEDIM == 2) { AverageCCToNode(eta,eta_ed[0],0,1,SPEC_BC_COMP,geom); } else { AverageCCToEdge(eta,eta_ed,0,1,SPEC_BC_COMP,geom); } - + ////////////////////////////////////////////// // Step 6: solve for v^{n+1} and pi^{n+1/2} using GMRES ////////////////////////////////////////////// @@ -677,7 +820,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, // call mk_grav_force_bousq(mla,gmres_rhs_v,.true.,rho_fc,the_bc_tower) } - if (use_multiphase || use_flory_huggins) { + if (use_multiphase || use_flory_huggins || use_ice_nucleation) { // add divergence of reversible stress to gmres_rhs_v for (int d=0; d& umac, eta_ed[d].mult(2.,0,1,0); } +if (use_ice_nucleation ==0) { // if writing a plotfile for electro-explicit option, compute Epot using new densities // use existing machinery to compute all mass fluxes - yes this is overkill // but better than writing a separate routine for now @@ -801,5 +945,6 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, charge_old,grad_Epot_old,Epot,permittivity,0); } } + } } diff --git a/exec/multispec/AdvanceTimestepBousq.cpp~ b/exec/multispec/AdvanceTimestepBousq.cpp~ new file mode 100644 index 000000000..1157e08a9 --- /dev/null +++ b/exec/multispec/AdvanceTimestepBousq.cpp~ @@ -0,0 +1,805 @@ + +#include "hydro_functions.H" +#include "common_functions.H" +#include "gmres_functions.H" +#include "multispec_functions.H" + +#include "StochMomFlux.H" + + +#include +#include +#include + + +// argv contains the name of the inputs file entered at the command line +void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, + MultiFab& rho_old, + MultiFab& rho_new, + MultiFab& rhotot_old, + MultiFab& rhotot_new, + MultiFab& pi, + MultiFab& eta, + std::array< MultiFab, NUM_EDGE >& eta_ed, + MultiFab& kappa, MultiFab& Temp, + std::array< MultiFab, NUM_EDGE >& Temp_ed, + MultiFab& diff_mass_fluxdiv, + MultiFab& stoch_mass_fluxdiv, + std::array< MultiFab, AMREX_SPACEDIM >& stoch_mass_flux, + std::array< MultiFab, AMREX_SPACEDIM >& grad_Epot_old, + std::array< MultiFab, AMREX_SPACEDIM >& grad_Epot_new, + MultiFab& charge_old, + MultiFab& charge_new, + MultiFab& Epot, + MultiFab& permittivity, + StochMassFlux& sMassFlux, + StochMomFlux& sMomFlux, + const Real& dt, + const Real& time, + const int& istep, + const Geometry& geom) +{ + + BL_PROFILE_VAR("AdvanceTimestepBousq()",AdvanceTimestepBousq); + + BoxArray ba = rho_old.boxArray(); + DistributionMapping dmap = rho_old.DistributionMap(); + + const Real* dx = geom.CellSize(); + + Vector weights_mom(1); + Vector weights_mass(2); + + if (barodiffusion_type != 0) { + Abort("AdvanceTimestepBousq: barodiffusion not supported yet"); + } + + // For BDS we need to project edge states onto constraints + // int proj_type = (use_charged_fluid && electroneutral) ? 4 : 3; + int proj_type = 3; + + Real theta_alpha = 1./dt; + + Real relxn_param_charge_in; + Real norm_pre_rhs; + + MultiFab adv_mass_fluxdiv(ba,dmap,nspecies,0); + MultiFab gmres_rhs_p (ba,dmap, 1,0); + MultiFab dpi (ba,dmap, 1,1); + + std::array< MultiFab, AMREX_SPACEDIM > umac_old; + std::array< MultiFab, AMREX_SPACEDIM > mtemp; + std::array< MultiFab, AMREX_SPACEDIM > adv_mom_fluxdiv_old; + std::array< MultiFab, AMREX_SPACEDIM > adv_mom_fluxdiv_new; + std::array< MultiFab, AMREX_SPACEDIM > diff_mom_fluxdiv_old; + std::array< MultiFab, AMREX_SPACEDIM > diff_mom_fluxdiv_new; + std::array< MultiFab, AMREX_SPACEDIM > stoch_mom_fluxdiv; + std::array< MultiFab, AMREX_SPACEDIM > gmres_rhs_v; + std::array< MultiFab, AMREX_SPACEDIM > dumac; + std::array< MultiFab, AMREX_SPACEDIM > gradpi; + std::array< MultiFab, AMREX_SPACEDIM > rho_fc; + std::array< MultiFab, AMREX_SPACEDIM > rhotot_fc_old; + std::array< MultiFab, AMREX_SPACEDIM > rhotot_fc_new; + std::array< MultiFab, AMREX_SPACEDIM > diff_mass_flux; + + // only used when variance_coef_mass>0 and midpoint_stoch_mass_flux_type=2 + MultiFab stoch_mass_fluxdiv_old; + std::array< MultiFab, AMREX_SPACEDIM > stoch_mass_flux_old; + + // only used when use_charged_fluid=1 + std::array< MultiFab, AMREX_SPACEDIM > Lorentz_force; + + // only used when use_multiphase=1 + std::array< MultiFab, AMREX_SPACEDIM > div_reversible_stress; + + for (int d=0; d umac_tmp; + + ////////////////////////////////////////////// + // Step 1: solve for v^{n+1,*} and pi^{n+1/2,*} using GMRES + ////////////////////////////////////////////// + + // average rho_i^n and rho^n to faces + AverageCCToFace( rho_old, rho_fc ,0,nspecies,SPEC_BC_COMP,geom); + AverageCCToFace(rhotot_old,rhotot_fc_old,0,1 , RHO_BC_COMP,geom); + + // compute mtemp = (rho*v)^n + ConvertMToUmac(rhotot_fc_old,umac,mtemp,0); + + // set gmres_rhs_v = mtemp / dt + for (int d=0; d 0) then + + end if + */ + + if (advection_type == 1 || advection_type == 2) { + + // bds advection + rho_update.define(ba,dmap,nspecies,0); + bds_force.define(ba,dmap,nspecies,1); + for (int d=0; d 0) { + // call multifab_plus_plus_c(rho_update(n),1,chem_rate(n),1,nspecies,0) + } + */ + + // set to zero to make sure ghost cells behind physical boundaries don't have NaNs + bds_force.setVal(0,0,1,1); + MultiFab::Copy(bds_force,rho_update,0,0,nspecies,0); + bds_force.FillBoundary(geom.periodicity()); + + // bds increments rho_update with the advection term + BDS(rho_update, nspecies, SPEC_BC_COMP, rho_old, umac_tmp, bds_force, geom, 0.5*dt, proj_type); + + } else if (advection_type == 0) { + + // compute adv_mass_fluxdiv = -rho_i^n * v^n and then + // increment adv_mass_fluxdiv by -rho_i^n * v^{n+1,*} + MkAdvSFluxdiv(umac_old,rho_fc,adv_mass_fluxdiv,geom,0,nspecies,false); + MkAdvSFluxdiv(umac ,rho_fc,adv_mass_fluxdiv,geom,0,nspecies,true); + } else { + + Abort("Invalid advection_type"); + + } + + ////////////////////////////////////////////// + /// Step 3: density prediction to t^{n+1/2} + ////////////////////////////////////////////// + + if (advection_type == 1 || advection_type == 2) { + + MultiFab::LinComb(rho_new,1.,rho_old,0,0.5*dt,rho_update,0,0,nspecies,0); + + } else { + + // compute rho_i^{n+1/2} (store in rho_new) + // multiply adv_mass_fluxdiv by (1/4) since it contains -rho_i^n * (v^n + v^{n+1,*}) + MultiFab::Copy(rho_new,rho_old,0,0,nspecies,0); + MultiFab::Saxpy(rho_new,0.25*dt, adv_mass_fluxdiv,0,0,nspecies,0); + MultiFab::Saxpy(rho_new,0.50*dt,diff_mass_fluxdiv,0,0,nspecies,0); + if (variance_coef_mass != 0.) { + MultiFab::Saxpy(rho_new,0.50*dt,stoch_mass_fluxdiv,0,0,nspecies,0); + } + /* + if (nreactions > 0) then + call multifab_saxpy_3_cc(rho_new(n),1,0.5d0*dt,chem_rate(n),1,nspecies) + end if + */ + } + + // compute rhotot^{n+1/2} from rho^{n+1/2} in VALID REGION + ComputeRhotot(rho_new,rhotot_new); + + // fill rho and rhotot ghost cells at t^{n+1/2} + FillRhoRhototGhost(rho_new,rhotot_new,geom); + + // average rho_i^{n+1/2} to faces + AverageCCToFace(rho_new,rho_fc,0,nspecies,SPEC_BC_COMP,geom); + + if (use_charged_fluid) { + // compute total charge at t^{n+1/2} + DotWithZ(rho_new,charge_new); + + // compute permittivity at t^{n+1/2} + if (dielectric_type != 0) { + Abort("AdvanceTimestepBousq dielectric_type != 0"); + } + } + + // compute (eta,kappa)^{n+1/2} + ComputeEta(rho_new, rhotot_new, eta); + if (AMREX_SPACEDIM == 2) { + AverageCCToNode(eta,eta_ed[0],0,1,SPEC_BC_COMP,geom); + } + else { + AverageCCToEdge(eta,eta_ed,0,1,SPEC_BC_COMP,geom); + } + + ////////////////////////////////////////////// + // Step 4: compute mass fluxes and reactions at t^{n+1/2} + ////////////////////////////////////////////// + + // compute mass fluxes and reactions at t^{n+1/2} + // For electroneutral, enable charge correction in the corrector + relxn_param_charge=relxn_param_charge_in; // Default value is 1 + + if (midpoint_stoch_mass_flux_type == 1) { + // strato + + // compute diffusive, stochastic, potential mass fluxes + // with barodiffusion and thermodiffusion + // this computes "-F = rho W chi [Gamma grad x... ]" + weights_mass[0] = 1./std::sqrt(2.); + weights_mass[1] = 1./std::sqrt(2.); + ComputeMassFluxdiv(rho_new,rhotot_new,Temp,diff_mass_fluxdiv,stoch_mass_fluxdiv, + diff_mass_flux,stoch_mass_flux,sMassFlux,dt,time,geom,weights_mass, + charge_new,grad_Epot_new,Epot,permittivity,0); + + } else if (midpoint_stoch_mass_flux_type == 2) { + + // ito + if (variance_coef_mass != 0.) { + // for ito interpretation we need to save stoch_mass_fluxdiv_old here + // then later add it to stoch_mass_fluxdiv and multiply by 1/2 + MultiFab::Copy(stoch_mass_fluxdiv_old,stoch_mass_fluxdiv,0,0,nspecies,0); + for (int d=0; d 0) then + + end if + */ + + if (advection_type == 1 || advection_type == 2) { + + // add the diff/stoch/react terms to rho_update + MultiFab::Copy(rho_update,diff_mass_fluxdiv,0,0,nspecies,0); + if (variance_coef_mass != 0.) { + MultiFab::Add(rho_update,stoch_mass_fluxdiv,0,0,nspecies,0); + } + /* + if (nreactions > 0) { + call multifab_plus_plus_c(rho_update(n),1,chem_rate(n),1,nspecies,0) + } + */ + + bds_force.setVal(0,0,1,1); + MultiFab::Copy(bds_force,rho_update,0,0,nspecies,0); + bds_force.FillBoundary(geom.periodicity()); + + BDS(rho_update, nspecies, SPEC_BC_COMP, rho_old, umac_tmp, bds_force, geom, dt, proj_type); + + } else if (advection_type == 0) { + // compute adv_mass_fluxdiv = -rho_i^{n+1/2} * v^n and + // increment adv_mass_fluxdiv by -rho_i^{n+1/2} * v^{n+1,*} + MkAdvSFluxdiv(umac_old,rho_fc,adv_mass_fluxdiv,geom,0,nspecies,false); + MkAdvSFluxdiv(umac ,rho_fc,adv_mass_fluxdiv,geom,0,nspecies,true); + } + + ////////////////////////////////////////////// + // Step 5: density integration to t^{n+1} + ////////////////////////////////////////////// + + if (advection_type == 1 || advection_type == 2) { + + MultiFab::LinComb(rho_new,1.,rho_old,0,0.5*dt,rho_update,0,0,nspecies,0); + + } else { + + // compute rho_i^{n+1} + // multiply adv_mass_fluxdiv by (1/2) since it contains -rho_i^{n+1/2} * (v^n + v^{n+1,*}) + MultiFab::Copy(rho_new,rho_old,0,0,nspecies,0); + MultiFab::Saxpy(rho_new,0.5*dt, adv_mass_fluxdiv,0,0,nspecies,0); + MultiFab::Saxpy(rho_new, dt,diff_mass_fluxdiv,0,0,nspecies,0); + if (variance_coef_mass != 0.) { + MultiFab::Saxpy(rho_new,dt,stoch_mass_fluxdiv,0,0,nspecies,0); + } + /* + if (nreactions > 0) then + call multifab_saxpy_3_cc(rho_new(n),1,dt,chem_rate(n),1,nspecies) + end if + */ + } + + if ( project_eos_int > 0 && istep%project_eos_int == 0) { + ProjectOntoEOS(rho_new); + } + + // compute rhotot^{n+1} from rho^{n+1} in VALID REGION + ComputeRhotot(rho_new,rhotot_new); + + // fill rho and rhotot ghost cells at t^{n+1} + FillRhoRhototGhost(rho_new,rhotot_new,geom); + + // average rho^{n+1} to faces + AverageCCToFace(rhotot_new,rhotot_fc_new,0,1,RHO_BC_COMP,geom); + + if (use_charged_fluid) { + // compute total charge at t^{n+1} + DotWithZ(rho_new,charge_new); + + // compute permittivity at t^{n+1} + if (dielectric_type != 0) { + Abort("AdvanceTimestepInertial dielectric_type != 0"); + } + } + + ComputeEta(rho_new, rhotot_new, eta); + if (AMREX_SPACEDIM == 2) { + AverageCCToNode(eta,eta_ed[0],0,1,SPEC_BC_COMP,geom); + } + else { + AverageCCToEdge(eta,eta_ed,0,1,SPEC_BC_COMP,geom); + } + + ////////////////////////////////////////////// + // Step 6: solve for v^{n+1} and pi^{n+1/2} using GMRES + ////////////////////////////////////////////// + + // compute mtemp = (rho*v)^n + ConvertMToUmac(rhotot_fc_old,umac_old,mtemp,0); + + // set gmres_rhs_v = mtemp / dt + for (int d=0; d 0 && ( (istep%plot_int == 0) || (istep == max_step)) ) { + // compute the new charge and store it in charge_old (it could get modified to + // subtract off the average in compute_mass_fluxdiv) + DotWithZ(rho_new,charge_old); + + ComputeMassFluxdiv(rho_new,rhotot_new,Temp,diff_mass_fluxdiv,stoch_mass_fluxdiv, + diff_mass_flux,stoch_mass_flux,sMassFlux,dt,time,geom,weights_mass, + charge_old,grad_Epot_old,Epot,permittivity,0); + } + } + +} diff --git a/exec/multispec/Checkpoint.cpp b/exec/multispec/Checkpoint.cpp index 06d20503c..4d46609c1 100644 --- a/exec/multispec/Checkpoint.cpp +++ b/exec/multispec/Checkpoint.cpp @@ -23,6 +23,8 @@ void WriteCheckPoint(int step, const amrex::Real dt, const MultiFab& rho, const MultiFab& rhotot, + const MultiFab& phi, + const MultiFab& phitot, const MultiFab& pi, std::array< MultiFab, AMREX_SPACEDIM >& umac, const MultiFab& Epot, @@ -145,6 +147,10 @@ void WriteCheckPoint(int step, amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "rho")); VisMF::Write(rhotot, amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "rhotot")); + VisMF::Write(phi, + amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "phi")); + VisMF::Write(phitot, + amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "phitot")); VisMF::Write(pi, amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "pi")); @@ -159,6 +165,8 @@ void ReadCheckPoint(int& step, amrex::Real& dt, MultiFab& rho, MultiFab& rhotot, + MultiFab& phi, + MultiFab& phitot, MultiFab& pi, std::array< MultiFab, AMREX_SPACEDIM >& umac, MultiFab& Epot, @@ -237,6 +245,8 @@ void ReadCheckPoint(int& step, rho .define(ba, dmap, nspecies, ng_s); rhotot.define(ba, dmap, 1, ng_s); + phi .define(ba, dmap, nspecies, ng_s); + phitot.define(ba, dmap, 1, ng_s); pi .define(ba, dmap, 1, 1); if (use_charged_fluid) { @@ -327,6 +337,10 @@ void ReadCheckPoint(int& step, amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "rho")); VisMF::Read(rhotot, amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "rhotot")); + VisMF::Read(phi, + amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "phi")); + VisMF::Read(phitot, + amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "phitot")); VisMF::Read(pi, amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "pi")); diff --git a/exec/multispec/GNUmakefile b/exec/multispec/GNUmakefile index 11584cd9d..f94107b34 100644 --- a/exec/multispec/GNUmakefile +++ b/exec/multispec/GNUmakefile @@ -7,13 +7,13 @@ USE_MPI = TRUE USE_OMP = FALSE USE_CUDA = FALSE COMP = gnu -DIM = 3 +DIM = 2 DSMC = FALSE MAX_SPEC = 8 # MAX_ELEM needs to be MAX_SPEC*(MAX_SPEC-1)/2 MAX_ELEM = 28 -TINY_PROFILE = TRUE +TINY_PROFILE = FALSE USE_PARTICLES = FALSE include $(AMREX_HOME)/Tools/GNUMake/Make.defs diff --git a/exec/multispec/Init.cpp b/exec/multispec/Init.cpp index 5a5e02e1b..da4671d4c 100644 --- a/exec/multispec/Init.cpp +++ b/exec/multispec/Init.cpp @@ -7,6 +7,7 @@ using namespace multispec; void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac, MultiFab& rho_in, + MultiFab& phi_in, const Geometry& geom) { @@ -35,7 +36,7 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac, Box bx = mfi.tilebox(); const Array4& c = conc.array(mfi); - + if (prob_type == 1) { /* @@ -78,7 +79,7 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac, }); - } else if (prob_type == 5) { + } else if (prob_type == 5) { /* bubble with radius = 1/4 of domain in x @@ -680,12 +681,12 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac, } else if (prob_type == 4) { /* - bubble with radius = 1/4 of domain in x + slab with thickness = 2*radius_cyl c=c_init_1(:) inside, c=c_init_2(:) outside can be discontinous or smooth depending on smoothing_width */ - Real l1 = L[1] / 4.; - Real l2 = 3.*l1; + Real l1 = L[1]*0.5-radius_cyl; + Real l2 = L[1]*0.5+radius_cyl; amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -796,6 +797,7 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac, const Array4& c = conc.array(mfi); const Array4& rho = rho_in.array(mfi); + const Array4& p = phi_in.array(mfi); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -820,10 +822,11 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac, } rho_total = 1./sumtot; } - - // calculate rho_i + + // calculate rho_i and make phi=conc for (int n=0; n grad_Epot_old; + /***/ + MultiFab phi_old; + MultiFab phitot_old; + if (restart > 0) { ReadCheckPoint(init_step,time,dt, - rho_old,rhotot_old,pi,umac,Epot,grad_Epot_old, + rho_old,rhotot_old,phi_old,phitot_old,pi,umac,Epot,grad_Epot_old, ba,dmap); } else { @@ -171,6 +175,10 @@ void main_driver(const char* argv) if (use_charged_fluid) { Epot.define(ba, dmap, 1, 1); } + + /***/ + phi_old.define (ba, dmap, nspecies, ng_s); + phitot_old.define (ba, dmap, 1 , ng_s); } // moved this to here so can change dt from value in checkpoint @@ -191,8 +199,9 @@ void main_driver(const char* argv) if (restart < 0) { + /***/ // initialize rho and umac in valid region only - InitRhoUmac(umac,rho_old,geom); + InitRhoUmac(umac,rho_old,phi_old,geom); // initialize pi, including ghost cells pi.setVal(0.); @@ -207,6 +216,9 @@ void main_driver(const char* argv) // pressure ghost cells pi.FillBoundary(geom.periodicity()); MultiFabPhysBC(pi,geom,0,1,PRES_BC_COMP); + /***/ + ComputeRhotot(phi_old,phitot_old); + FillRhoRhototGhost(phi_old,phitot_old,geom); //======================================================= // Build multifabs for all the variables @@ -219,6 +231,10 @@ void main_driver(const char* argv) MultiFab eta (ba, dmap, 1 , 1); MultiFab kappa (ba, dmap, 1 , 1); + /***/ + MultiFab phi_new (ba, dmap, nspecies, ng_s); + MultiFab phitot_new (ba, dmap, 1 , ng_s); + ///////////////////////////////////////// // eta and Temp on nodes (2d) or edges (3d) @@ -506,7 +522,7 @@ void main_driver(const char* argv) // write initial plotfile and structure factor if (plot_int > 0) { - WritePlotFile(0,0.,geom,umac,rhotot_old,rho_old,pi,charge_old,Epot); + WritePlotFile(0,0.,geom,umac,rhotot_old,rho_old,pi,charge_old,Epot,phitot_old,phi_old); if (n_steps_skip == 0 && struct_fact_int > 0) { structFact.WritePlotFile(0,0.,geom,"plt_SF"); } @@ -548,12 +564,14 @@ void main_driver(const char* argv) } else if (algorithm_type == 6) { // boussinesq + /***/ AdvanceTimestepBousq(umac,rho_old,rho_new,rhotot_old,rhotot_new, pi,eta,eta_ed,kappa,Temp,Temp_ed, diff_mass_fluxdiv,stoch_mass_fluxdiv,stoch_mass_flux, grad_Epot_old,grad_Epot_new, charge_old,charge_new,Epot,permittivity, sMassFlux,sMomFlux, + phi_old,phi_new,phitot_old,phitot_new, dt,time,istep,geom); } else { @@ -591,7 +609,7 @@ void main_driver(const char* argv) // write plotfile at specific intervals if (plot_int > 0 && istep%plot_int == 0) { - WritePlotFile(istep,time,geom,umac,rhotot_new,rho_new,pi,charge_new,Epot); + WritePlotFile(istep,time,geom,umac,rhotot_new,rho_new,pi,charge_new,Epot,phitot_new,phi_new); if (istep > n_steps_skip && struct_fact_int > 0) { structFact.WritePlotFile(istep,time,geom,"plt_SF"); } @@ -599,7 +617,7 @@ void main_driver(const char* argv) // write checkpoint at specific intervals if (chk_int > 0 && istep%chk_int == 0) { - WriteCheckPoint(istep,time,dt,rho_new,rhotot_new,pi,umac,Epot,grad_Epot_new); + WriteCheckPoint(istep,time,dt,rho_new,rhotot_new,phi_new,phitot_new,pi,umac,Epot,grad_Epot_new); } // set old state to new state diff --git a/exec/multispec/multispec_test_functions.H b/exec/multispec/multispec_test_functions.H index 2c0cb2b89..0dac6e50e 100644 --- a/exec/multispec/multispec_test_functions.H +++ b/exec/multispec/multispec_test_functions.H @@ -60,6 +60,10 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac, MultiFab& permittivity, StochMassFlux& sMassFlux, StochMomFlux& sMomFlux, + MultiFab& phi_old, + MultiFab& phi_new, + MultiFab& phitot_old, + MultiFab& phitot_new, const Real& dt, const Real& time, const int& istep, @@ -73,6 +77,8 @@ void WriteCheckPoint(int step, const amrex::Real dt, const MultiFab& rho, const MultiFab& rhotot, + const MultiFab& phi, + const MultiFab& phitot, const MultiFab& pi, std::array< MultiFab, AMREX_SPACEDIM >& umac, const MultiFab& Epot, @@ -83,6 +89,8 @@ void ReadCheckPoint(int& step, amrex::Real& dt, MultiFab& rho, MultiFab& rhotot, + MultiFab& phi, + MultiFab& phitot, MultiFab& pi, std::array< MultiFab, AMREX_SPACEDIM >& umac, MultiFab& Epot, @@ -98,6 +106,7 @@ void ReadFile(const std::string& filename, Vector& charBuf, void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac, MultiFab& rho_in, + MultiFab& phi_old, const Geometry& geom); /////////////////////////// @@ -108,10 +117,12 @@ void WritePlotFile(const int step, const Geometry geom, std::array< MultiFab, AMREX_SPACEDIM >& umac, const MultiFab& rhotot, - const MultiFab& rho, + const MultiFab& rho, const MultiFab& pres, const MultiFab& charge, - const MultiFab& Epot); + const MultiFab& Epot, + const MultiFab& phitot, + const MultiFab& phi); /////////////////////////// diff --git a/src_multispec/InitialProjection.cpp b/src_multispec/InitialProjection.cpp index 4aa1a869a..7f7567df5 100644 --- a/src_multispec/InitialProjection.cpp +++ b/src_multispec/InitialProjection.cpp @@ -95,11 +95,12 @@ void InitialProjection(std::array< MultiFab, AMREX_SPACEDIM >& umac, if (variance_coef_mass != 0. && algorithm_type != 6) { sMassFlux.fillMassStochastic(); } - + +if (use_ice_nucleation ==0) { ComputeMassFluxdiv(rho,rhotot,Temp,diff_mass_fluxdiv,stoch_mass_fluxdiv, diff_mass_flux,stoch_mass_flux,sMassFlux,dt_eff,time,geom,weights, charge_old,grad_Epot_old,Epot,permittivity); - +} // assumble total fluxes to be used in reservoirs // // diff --git a/src_multispec/multispec_functions.cpp b/src_multispec/multispec_functions.cpp index 398fa4e52..cdaf91b80 100644 --- a/src_multispec/multispec_functions.cpp +++ b/src_multispec/multispec_functions.cpp @@ -18,6 +18,8 @@ AMREX_GPU_MANAGED int multispec::is_ideal_ int multispec::use_lapack; AMREX_GPU_MANAGED int multispec::use_multiphase; AMREX_GPU_MANAGED int multispec::use_flory_huggins; +AMREX_GPU_MANAGED int multispec::use_ice_nucleation; +AMREX_GPU_MANAGED int multispec::use_log_potential; AMREX_GPU_MANAGED amrex::Real multispec::kc_tension; AMREX_GPU_MANAGED amrex::Real multispec::alpha_gex; AMREX_GPU_MANAGED amrex::Array2D multispec::fh_kappa; @@ -83,6 +85,8 @@ void InitializeMultispecNamespace() { use_lapack = 0; // Use LAPACK or iterative method for diffusion matrix (recommend False) use_multiphase = 0; // for RTIL use_flory_huggins = 0; // for flory huggins + use_ice_nucleation = 0; // for ice nucleation simulations + use_log_potential = 0; // use quartic potential by default kc_tension = 0; // for RTIL alpha_gex = 0; // for RTIL n_gex = 1; // for RTIL @@ -192,6 +196,8 @@ void InitializeMultispecNamespace() { pp.query("use_lapack",use_lapack); pp.query("use_multiphase",use_multiphase); pp.query("use_flory_huggins",use_flory_huggins); + pp.query("use_ice_nucleation",use_ice_nucleation); + pp.query("use_log_potential",use_log_potential); pp.query("monomer_mass",monomer_mass); pp.query("kc_tension",kc_tension); pp.query("fh_tension",fh_tension); diff --git a/src_multispec/multispec_namespace.H b/src_multispec/multispec_namespace.H index 6e4819833..806f1b7a3 100644 --- a/src_multispec/multispec_namespace.H +++ b/src_multispec/multispec_namespace.H @@ -21,6 +21,8 @@ namespace multispec { extern AMREX_GPU_MANAGED int use_lapack; extern AMREX_GPU_MANAGED int use_multiphase; extern AMREX_GPU_MANAGED int use_flory_huggins; + extern AMREX_GPU_MANAGED int use_ice_nucleation; + extern AMREX_GPU_MANAGED int use_log_potential; extern AMREX_GPU_MANAGED amrex::Array2D fh_kappa; extern AMREX_GPU_MANAGED amrex::Array2D fh_chi; extern AMREX_GPU_MANAGED amrex::GpuArray fh_monomers;