Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Ice nucleation #149

Closed
wants to merge 6 commits into from
Closed

Ice nucleation #149

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
ef06864
Updated from main
Mar 19, 2024
7da32c2
Updated from main
Mar 19, 2024
d2fc679
Mobility scale array
Mar 21, 2024
d2721f8
Mobility scale array
Mar 21, 2024
4db4924
Mobility scale array deleted
Mar 24, 2024
fe16e61
Mobility scale array deleted
Mar 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
163 changes: 153 additions & 10 deletions exec/multispec/AdvanceTimestepBousq.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@

#include <AMReX_ParallelDescriptor.H>
#include <AMReX_MultiFabUtil.H>
#include <AMReX_Print.H>


// argv contains the name of the inputs file entered at the command line
Expand All @@ -33,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,
Expand Down Expand Up @@ -122,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<AMREX_SPACEDIM; ++d) {
div_reversible_stress[d].define(convert(ba,nodal_flag_dir[d]), dmap, 1, 0);
}
Expand Down Expand Up @@ -199,7 +204,8 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
if (any_grav) {
Abort("AdvanceTimestepBousq.cpp gravity not implemented");
}


if (use_ice_nucleation ==0) {
if (variance_coef_mass != 0.) {
sMassFlux.fillMassStochastic();
}
Expand All @@ -217,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) {
Expand All @@ -240,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<AMREX_SPACEDIM; ++d) {
MultiFab::Saxpy(gmres_rhs_v[d],1.,div_reversible_stress[d],0,0,1,0);
}

}

if (use_charged_fluid) {

Expand Down Expand Up @@ -358,6 +376,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
eta_ed[d].mult(2.,0,1,0);
}

if (use_ice_nucleation ==0){
//////////////////////////////////////////////
// Step 2: compute reactions and mass fluxes at t^n
//////////////////////////////////////////////
Expand All @@ -368,7 +387,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,

end if
*/

if (advection_type == 1 || advection_type == 2) {

// bds advection
Expand Down Expand Up @@ -415,7 +434,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
Abort("Invalid advection_type");

}

//////////////////////////////////////////////
/// Step 3: density prediction to t^{n+1/2}
//////////////////////////////////////////////
Expand All @@ -440,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);

Expand All @@ -450,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);
Expand All @@ -461,14 +491,19 @@ 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);
}
else {
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}
//////////////////////////////////////////////
Expand Down Expand Up @@ -525,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);

}

Expand Down Expand Up @@ -568,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}
//////////////////////////////////////////////
Expand Down Expand Up @@ -606,6 +641,108 @@ 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);

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<Real>& phiOld = phi_old.array(mfi);
const Array4<Real>& phiPrd = phi_prd.array(mfi);

const amrex::Array4<amrex::Real>& u = umac[0].array(mfi);
const amrex::Array4<amrex::Real>& v = umac[1].array(mfi);
const amrex::Array4<amrex::Real>& u_old = umac_old[0].array(mfi);
const amrex::Array4<amrex::Real>& v_old = umac_old[1].array(mfi);
#if (AMREX_SPACEDIM == 3)
const amrex::Array4<amrex::Real>& w = umac[2].array(mfi);
const amrex::Array4<amrex::Real>& 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<Real>& phiNew = phi_new.array(mfi);
const Array4<Real>& phiOld = phi_old.array(mfi);
const Array4<Real>& phiPrd = phi_prd.array(mfi);

const amrex::Array4<amrex::Real>& u = umac[0].array(mfi);
const amrex::Array4<amrex::Real>& v = umac[1].array(mfi);
const amrex::Array4<amrex::Real>& u_old = umac_old[0].array(mfi);
const amrex::Array4<amrex::Real>& v_old = umac_old[1].array(mfi);
#if (AMREX_SPACEDIM == 3)
const amrex::Array4<amrex::Real>& w = umac[2].array(mfi);
const amrex::Array4<amrex::Real>& 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() << "Mobility = " << Mobility <<"\n";

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}
Expand All @@ -617,14 +754,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
//////////////////////////////////////////////
Expand Down Expand Up @@ -676,7 +817,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<AMREX_SPACEDIM; ++d) {
MultiFab::Saxpy(gmres_rhs_v[d],1.,div_reversible_stress[d],0,0,1,0);
Expand Down Expand Up @@ -783,6 +924,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& 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
Expand All @@ -800,5 +942,6 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
charge_old,grad_Epot_old,Epot,permittivity,0);
}
}
}

}
14 changes: 14 additions & 0 deletions exec/multispec/Checkpoint.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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"));

Expand All @@ -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,
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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"));

Expand Down
2 changes: 1 addition & 1 deletion exec/multispec/GNUmakefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ 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
Expand Down
Loading
Loading