Skip to content

Commit

Permalink
fixed bug in contact angle boundary conditions
Browse files Browse the repository at this point in the history
  • Loading branch information
@jbbel
jbbel committed Jan 11, 2024
1 parent dd332e2 commit 5272e05
Show file tree
Hide file tree
Showing 6 changed files with 65 additions and 37 deletions.
1 change: 1 addition & 0 deletions exec/multispec/AdvanceTimestepBousq.cpp
View file
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 Down
31 changes: 29 additions & 2 deletions src_common/BCPhysToMath.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,35 @@ void BCPhysToMath(int bccomp, amrex::Vector<int>& bc_lo, amrex::Vector<int>& bc_
}
}
}
else if (bccomp == RHO_BC_COMP ||
bccomp == SPEC_BC_COMP ||
else if (bccomp == RHO_BC_COMP ) { // density
for (int i=0; i<AMREX_SPACEDIM; ++i) {
if (bc_mass_lo[i] == 1) {
// wall -> first-order extrapolation
bc_lo[i] = FOEXTRAP;
}
else if (bc_mass_lo[i] == 2) {
// reservoir -> dirichlet
bc_lo[i] = EXT_DIR;
}
else if (bc_mass_lo[i] == 4) {
// reservoir -> dirichlet
bc_lo[i] = FOEXTRAP;
}
if (bc_mass_hi[i] == 1) {
// wall -> first-order extrapolation
bc_hi[i] = FOEXTRAP;
}
else if (bc_mass_hi[i] == 2) {
// reservoir -> dirichlet
bc_hi[i] = EXT_DIR;
}
else if (bc_mass_hi[i] == 4) {
// reservoir -> dirichlet
bc_hi[i] = FOEXTRAP;
}
}
}
else if (bccomp == SPEC_BC_COMP ||
bccomp == MOLFRAC_BC_COMP) { // density, species, or mole fraction
for (int i=0; i<AMREX_SPACEDIM; ++i) {
if (bc_mass_lo[i] == 1) {
Expand Down
2 changes: 1 addition & 1 deletion src_common/MultiFabPhysBC.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,7 +250,7 @@ void MultiFabPhysBC(MultiFab& phi, const Geometry& geom, int scomp, int ncomp, i
}
});

} esle if (bc_lo[2] == FOEXTRAP) {
} else if (bc_lo[2] == FOEXTRAP) {
amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
if (k < lo) {
Expand Down
40 changes: 20 additions & 20 deletions src_multispec/MultiFabPhysBC_FH.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ void MultiFabPhysBCFH(MultiFab& phi, const Geometry& geom, int scomp, int ncomp,

// bccomp definitions are in BCPhysToMath.cpp

amrex::Print() << " entering special FH bc routine " << std::endl;
// amrex::Print() << " entering special FH bc routine " << std::endl;

if (geom.isAllPeriodic() || phi.nGrow() == 0) {
return;
Expand Down Expand Up @@ -55,7 +55,7 @@ void MultiFabPhysBCFH(MultiFab& phi, const Geometry& geom, int scomp, int ncomp,
if (bx.smallEnd(0) < lo) {
Real x = prob_lo[0];
if (bc_mass_lo[0] == 1) {
amrex::Print() << " entering special FH bc routine lo " << bc_mass_lo[0] << std::endl;
// amrex::Print() << " entering special FH bc routine lo " << bc_mass_lo[0] << std::endl;
amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
if (i < lo) {
Expand All @@ -66,27 +66,27 @@ void MultiFabPhysBCFH(MultiFab& phi, const Geometry& geom, int scomp, int ncomp,
});
}
else if (bc_mass_lo[0] == 4) {
amrex::Print() << " entering special FH bc routine lo " << bc_mass_lo[0] << std::endl;
// amrex::Print() << " entering special FH bc routine lo " << bc_mass_lo[0] << std::endl;
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
if (i < lo) {
data(i,j,k,scomp+0) = data(lo,j,k,scomp+0) + dx[0]*std::cos(contact_angle_lo[0])*data(lo,j,k,scomp+0)*data(lo,j,k,scomp+1)*coeff;
data(i,j,k,scomp+0) = amrex::min(1.,amrex::max(0.,data(i,j,k,scomp+0)));
data(i,j,k,scomp+1) = 1.-data(i,j,k,scomp+0);
amrex::Print() << " Fh data left " << j << " " << data(i,j,k,scomp) <<" " << data(i,j,k,scomp+1) << std::endl;
// amrex::Print() << " Fh data left " << j << " " << data(i,j,k,scomp) <<" " << data(i,j,k,scomp+1) << std::endl;

if(j == 18){
Real foo = std::cos(contact_angle_lo[0])* coeff;
amrex::Print() << "coeff and cos " << coeff << " " << contact_angle_lo[0] << " " << foo << std::endl;
}
// if(j == 18){
// Real foo = std::cos(contact_angle_lo[0])* coeff;
// amrex::Print() << "coeff and cos " << coeff << " " << contact_angle_lo[0] << " " << foo << std::endl;
// }

//data(i,j,k,scomp+1) = data(lo,j,k,scomp+1) + 0.5*dx[0]*data(lo,j,k,scomp+0)*data(lo,j,k,scomp+1);
}
});

}
else{
amrex::Print{} << "Should not be here " << std::endl;
else if ( bc_mass_lo[0] != -1) {
amrex::Print{} << "Should not be here lo x " << bc_mass_lo[0] << std::endl;
Abort("Wrong parameter in FH boundary condition");
}
}
Expand All @@ -104,22 +104,22 @@ void MultiFabPhysBCFH(MultiFab& phi, const Geometry& geom, int scomp, int ncomp,
});
}
else if (bc_mass_hi[0] == 4) {
amrex::Print() << " entering special FH bc routine hi " << bc_mass_hi[0] << std::endl;
// amrex::Print() << " entering special FH bc routine hi " << bc_mass_hi[0] << std::endl;
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
if (i > hi) {
data(i,j,k,scomp+0) = data(hi,j,k,scomp+0) + dx[0]*std::cos(contact_angle_hi[0])*data(hi,j,k,scomp+0)*data(hi,j,k,scomp+1)*coeff;
data(i,j,k,scomp+0) = amrex::min(1.,amrex::max(0.,data(i,j,k,scomp+0)));
data(i,j,k,scomp+1) = 1.-data(i,j,k,scomp+0);
amrex::Print() << " Fh data right " << j << " " << data(i,j,k,scomp) << " " << data(i,j,k,scomp+1) << std::endl;
// amrex::Print() << " Fh data right " << j << " " << data(i,j,k,scomp) << " " << data(i,j,k,scomp+1) << std::endl;

//data(i,j,k,scomp+1) = data(lo,j,k,scomp+1) + 0.5*dx[0]*data(lo,j,k,scomp+0)*data(lo,j,k,scomp+1);
}
});

}
else{
amrex::Print{} << "Should not be here " << std::endl;
else if ( bc_mass_hi[0] != -1) {
amrex::Print{} << "Should not be here hi x " << bc_mass_hi[0] << std::endl;
Abort("Wrong parameter in FH boundary condition");
}
}
Expand Down Expand Up @@ -156,8 +156,8 @@ void MultiFabPhysBCFH(MultiFab& phi, const Geometry& geom, int scomp, int ncomp,
});

}
else{
amrex::Print{} << "Should not be here " << std::endl;
else if ( bc_mass_lo[1] != -1) {
amrex::Print{} << "Should not be here lo y " << bc_mass_lo[1] << std::endl;
Abort("Wrong parameter in FH boundary condition");
}
}
Expand Down Expand Up @@ -187,8 +187,8 @@ void MultiFabPhysBCFH(MultiFab& phi, const Geometry& geom, int scomp, int ncomp,
});

}
else{
amrex::Print{} << "Should not be here " << std::endl;
else if ( bc_mass_hi[1] != -1) {
amrex::Print{} << "Should not be here hi y " << bc_mass_hi[1] << std::endl;
Abort("Wrong parameter in FH boundary condition");
}
}
Expand Down Expand Up @@ -226,7 +226,7 @@ void MultiFabPhysBCFH(MultiFab& phi, const Geometry& geom, int scomp, int ncomp,
});

}
else{
else if ( bc_mass_lo[2] != -1) {
amrex::Print{} << "Should not be here " << std::endl;
Abort("Wrong parameter in FH boundary condition");
}
Expand Down Expand Up @@ -257,7 +257,7 @@ void MultiFabPhysBCFH(MultiFab& phi, const Geometry& geom, int scomp, int ncomp,
});

}
else{
else if ( bc_mass_hi[2] != -1) {
amrex::Print{} << "Should not be here " << std::endl;
Abort("Wrong parameter in FH boundary condition");
}
Expand Down
2 changes: 0 additions & 2 deletions src_multispec/RhoUtil.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -196,7 +196,6 @@ void FillRhototGhost(MultiFab& rhotot_in, const MultiFab& conc_in, const Geometr
rhoinv += conc(lo-1,j,k,n)/rhobar_gpu[n];
}
rhotot(i,j,k) = 1./rhoinv;
amrex::Print() << " rhotot lo" << j << " " << rhotot(i,j,k) << std::endl;
}
});
}
Expand All @@ -212,7 +211,6 @@ void FillRhototGhost(MultiFab& rhotot_in, const MultiFab& conc_in, const Geometr
rhoinv += conc(hi+1,j,k,n)/rhobar_gpu[n];
}
rhotot(i,j,k) = 1./rhoinv;
amrex::Print() << " rhotot hi" << j << " " << rhotot(i,j,k) << std::endl;
}
});
}
Expand Down
26 changes: 14 additions & 12 deletions src_multispec/StochMassFlux.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,11 +70,11 @@ void StochMassFlux::StochMassFluxBC() {
BL_PROFILE_VAR("StochMassFluxBC()",StochMassFluxBC);

// lo-x domain boundary
if (bc_mass_lo[0] == 1 || bc_mass_lo[0] == 2) {
if (bc_mass_lo[0] == 1 || bc_mass_lo[0] == 2 || bc_mass_lo[0] == 4) {

// 1 = wall : multiply fluxes on wall by 0
// 2 = reservoir : multiply fluxes on wall by sqrt(2)
Real factor = (bc_mass_lo[0] == 1) ? 0. : sqrt(2.);
Real factor = (bc_mass_lo[0] == 1 || bc_mass_lo[0] == 4) ? 0. : sqrt(2.);

// domain grown nodally based on stoch_W_fc_weighted[0] nodality (x)
const Box& dom_x = amrex::convert(geom.Domain(), stoch_W_fc_weighted[0].ixType());
Expand All @@ -100,11 +100,11 @@ void StochMassFlux::StochMassFluxBC() {
}

// hi-x domain boundary
if (bc_mass_hi[0] == 1 || bc_mass_hi[0] == 2) {
if (bc_mass_hi[0] == 1 || bc_mass_hi[0] == 2 || bc_mass_hi[0] == 4) {

// 1 = wall : multiply fluxes on wall by 0
// 2 = reservoir : multiply fluxes on wall by sqrt(2)
Real factor = (bc_mass_hi[0] == 1) ? 0. : sqrt(2.);
Real factor = (bc_mass_hi[0] == 1 || bc_mass_hi[0] == 4) ? 0. : sqrt(2.);

// domain grown nodally based on stoch_W_fc_weighted[0] nodality (x)
const Box& dom_x = amrex::convert(geom.Domain(), stoch_W_fc_weighted[0].ixType());
Expand All @@ -131,11 +131,11 @@ void StochMassFlux::StochMassFluxBC() {
}

// lo-y domain boundary
if (bc_mass_lo[1] == 1 || bc_mass_lo[1] == 2) {
if (bc_mass_lo[1] == 1 || bc_mass_lo[1] == 2 || bc_mass_lo[1] == 4) {

// 1 = wall : multiply fluxes on wall by 0
// 2 = reservoir : multiply fluxes on wall by sqrt(2)
Real factor = (bc_mass_lo[1] == 1) ? 0. : sqrt(2.);
Real factor = (bc_mass_lo[1] == 1 || bc_mass_lo[1] == 4) ? 0. : sqrt(2.);

// domain grown nodally based on stoch_W_fc_weighted[1] nodality (y)
const Box& dom_y = amrex::convert(geom.Domain(), stoch_W_fc_weighted[1].ixType());
Expand All @@ -161,11 +161,11 @@ void StochMassFlux::StochMassFluxBC() {
}

// hi-y domain boundary
if (bc_mass_hi[1] == 1 || bc_mass_hi[1] == 2) {
if (bc_mass_hi[1] == 1 || bc_mass_hi[1] == 2 || bc_mass_hi[1] == 4) {

// 1 = wall : multiply fluxes on wall by 0
// 2 = reservoir : multiply fluxes on wall by sqrt(2)
Real factor = (bc_mass_hi[1] == 1) ? 0. : sqrt(2.);
Real factor = (bc_mass_hi[1] == 1 || bc_mass_hi[1] == 4) ? 0. : sqrt(2.);

// domain grown nodally based on stoch_W_fc_weighted[1] nodality (y)
const Box& dom_y = amrex::convert(geom.Domain(), stoch_W_fc_weighted[1].ixType());
Expand Down Expand Up @@ -193,11 +193,11 @@ void StochMassFlux::StochMassFluxBC() {
#if (AMREX_SPACEDIM == 3)

// lo-z domain boundary
if (bc_mass_lo[2] == 1 || bc_mass_lo[2] == 2) {
if (bc_mass_lo[2] == 1 || bc_mass_lo[2] == 2 || bc_mass_lo[2] == 4) {

// 1 = wall : multiply fluxes on wall by 0
// 2 = reservoir : multiply fluxes on wall by sqrt(2)
Real factor = (bc_mass_lo[2] == 1) ? 0. : sqrt(2.);
Real factor = (bc_mass_lo[2] == 1 || bc_mass_lo[2] == 4) ? 0. : sqrt(2.);

// domain grown nodally based on stoch_W_fc_weighted[2] nodality (z)
const Box& dom_z = amrex::convert(geom.Domain(), stoch_W_fc_weighted[2].ixType());
Expand All @@ -223,11 +223,11 @@ void StochMassFlux::StochMassFluxBC() {
}

// hi-z domain boundary
if (bc_mass_hi[2] == 1 || bc_mass_hi[2] == 2) {
if (bc_mass_hi[2] == 1 || bc_mass_hi[2] == 2 || bc_mass_hi[2] == 4) {

// 1 = wall : multiply fluxes on wall by 0
// 2 = reservoir : multiply fluxes on wall by sqrt(2)
Real factor = (bc_mass_hi[2] == 1) ? 0. : sqrt(2.);
Real factor = (bc_mass_hi[2] == 1 || bc_mass_hi[2] == 4) ? 0. : sqrt(2.);

// domain grown nodally based on stoch_W_fc_weighted[2] nodality (z)
const Box& dom_z = amrex::convert(geom.Domain(), stoch_W_fc_weighted[2].ixType());
Expand Down Expand Up @@ -280,9 +280,11 @@ void StochMassFlux::StochMassFluxDiv(const MultiFab& rho,
MultiFab::Copy(stoch_mass_flux[d], stoch_W_fc_weighted[d], 0, 0, nspecies, 0);
}


const Real* dx = geom.CellSize();
Real dVol = (AMREX_SPACEDIM==2) ? dx[0]*dx[1]*cell_depth : dx[0]*dx[1]*dx[2];
Real variance = sqrt(2.*k_B*variance_coef_mass/(dVol*dt));


// compute variance X sqrtLonsager_fc X stoch_mass_flux X variance
for (int d=0; d<AMREX_SPACEDIM; ++d) {
Expand Down

0 comments on commit 5272e05

Please sign in to comment.