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

FH Attempt 3 #153

Closed
wants to merge 13 commits into from
Closed

FH Attempt 3 #153

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
481f5bd
tweaks to initial conditions
jbbel Oct 19, 2022
7bebd68
Merge branch 'flory_huggins' of github.com:AMReX-FHD/FHDeX into flory…
jbbel Nov 13, 2023
379d721
Merge branch 'main' into flory_huggins
jbbel Nov 13, 2023
dd332e2
added support for contact angle boundary conditions for flory huggins…
jbbel Dec 1, 2023
5272e05
fixed bug in contact angle boundary conditions
jbbel Jan 11, 2024
19dcda3
added disjoining pressure and fixed contact angle BC
jbbel Jan 26, 2024
34230db
fix disjoining pressure bug
jbbel Jan 26, 2024
905bc01
adjust range in contact anglre bc to use ce to 1-ce
jbbel Feb 2, 2024
2bd05e1
tweaked input
jbbel Feb 15, 2024
97b1b21
set up hemisphere initial data
jbbel Mar 28, 2024
95903b5
compilation fixes
ajnonaka Apr 2, 2024
fa397c7
avoiding merge conflicts
ajnonaka Apr 2, 2024
9eb8dc1
Merge branch 'main' into flory_huggins
ajnonaka Apr 3, 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
60 changes: 60 additions & 0 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 Down Expand Up @@ -41,6 +42,8 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,

BL_PROFILE_VAR("AdvanceTimestepBousq()",AdvanceTimestepBousq);

// int use_disjoin_pres = 0;

BoxArray ba = rho_old.boxArray();
DistributionMapping dmap = rho_old.DistributionMap();

Expand Down Expand Up @@ -88,6 +91,8 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
// only used when use_charged_fluid=1
std::array< MultiFab, AMREX_SPACEDIM > Lorentz_force;

std::array< MultiFab, AMREX_SPACEDIM > disjoining_pressure;

// only used when use_multiphase=1
std::array< MultiFab, AMREX_SPACEDIM > div_reversible_stress;

Expand Down Expand Up @@ -127,6 +132,11 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
div_reversible_stress[d].define(convert(ba,nodal_flag_dir[d]), dmap, 1, 0);
}
}
if (use_disjoin_pres == 1) {
for (int d=0; d<AMREX_SPACEDIM; ++d) {
disjoining_pressure[d].define(convert(ba,nodal_flag_dir[d]), dmap, 1, 0);
}
}

// make a copy of umac at t^n
for (int d=0; d<AMREX_SPACEDIM; ++d) {
Expand Down Expand Up @@ -253,6 +263,28 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,

}

if (use_disjoin_pres == 1) {

//amrex::Print() << " computing disjoining pressure " << std::endl;
// compute old disjoining pressure
ComputeDisjoiningPressure(disjoining_pressure,rhotot_old,rho_old,geom);

// add dijoining pressure gradient to gmres_rhs_v
for (int d=0; d<AMREX_SPACEDIM; ++d) {
//Print() <<" BEFORE RHS PRED " << d << " " << std::endl;
//print_state(gmres_rhs_v[d],IntVect(255,25));
//Print() <<" PRES GRAD BEFORE PRED " << std::endl;
//print_state(disjoining_pressure[d],IntVect(255,25));

MultiFab::Saxpy(gmres_rhs_v[d],1.0,disjoining_pressure[d],0,0,1,0);

//Print() <<" RHS AFTER PRED " << std::endl;
//print_state(gmres_rhs_v[d],IntVect(255,25));
//Print() <<" " << std::endl;
}

}

// compute grad pi^{n-1/2}
ComputeGrad(pi,gradpi,0,0,1,PRES_BC_COMP,geom);

Expand Down Expand Up @@ -536,6 +568,15 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,

}

if (use_disjoin_pres == 1) {

//amrex::Print() << " dj for corrector " << std::endl;

// compute old disjoining pressure
ComputeDisjoiningPressure(disjoining_pressure,rhotot_new,rho_new,geom);

}

// compute chemical rates m_i*R^{n+1/2}_i
/*
if (nreactions > 0) then
Expand Down Expand Up @@ -692,6 +733,25 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,

}

if (use_disjoin_pres == 1) {

//amrex::Print() << " add dj in corrector " << std::endl;
// add dijoining pressure gradient to gmres_rhs_v
for (int d=0; d<AMREX_SPACEDIM; ++d) {
//Print() <<" BEFORE RHS " << d << " " << std::endl;
//print_state(gmres_rhs_v[d],IntVect(255,25));
//Print() <<" PRES GRAD BEFORE " << std::endl;
//print_state(disjoining_pressure[d],IntVect(255,25));

MultiFab::Saxpy(gmres_rhs_v[d],1.0,disjoining_pressure[d],0,0,1,0);

//Print() <<" RHS AFTER " << std::endl;
//print_state(gmres_rhs_v[d],IntVect(255,25));
//Print() <<" " << std::endl;
}

}

// compute grad pi^{n+1/2,*}
ComputeGrad(pi,gradpi,0,0,1,PRES_BC_COMP,geom);

Expand Down
76 changes: 74 additions & 2 deletions exec/multispec/Init.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,13 +45,23 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac,
*/
Real rad = L[0] / 4.;
rad = radius_cyl;

amrex::Real alpha;

alpha = 3.*3.14159265/4.;

alpha = contact_angle_lo[1];

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real x,y,z;
// AMREX_D_TERM(x = prob_lo[0] + (i+0.5)*dx[0] - center[0];,
// y = prob_lo[1] + (j+0.5)*dx[1] - center[1];,
// z = prob_lo[2] + (k+0.5)*dx[2] - center[2];);
AMREX_D_TERM(x = prob_lo[0] + (i+0.5)*dx[0] - center[0];,
y = prob_lo[1] + (j+0.5)*dx[1] - center[1];,
y = prob_lo[1] + (j+0.5)*dx[1] ;,
z = prob_lo[2] + (k+0.5)*dx[2] - center[2];);
//y = prob_lo[1] + (j+0.5)*dx[1] - rad*std::cos(alpha) ;,

Real r = (AMREX_SPACEDIM == 1) ? std::sqrt(x*x+y*y) : std::sqrt(x*x+y*y+z*z);

Expand Down Expand Up @@ -493,6 +503,68 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac,
}
});

} else if (prob_type == 21) {

/*
cylinder with radius = 1/4 of domain in x
c=c_init_1(:) inside, c=c_init_2(:) outside
can be discontinous or smooth depending on smoothing_width
*/
//Real rad = L[0] / 8.;
int nsub = 100;
Real factor = nsub;
Real dxsub = dx[0]/factor;
Real dysub = dx[1]/factor;
Real x,y,z;
amrex::Print() << "smoothing width " << smoothing_width << " film_thickness " << film_thickness << std::endl;

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
for (int n=0; n<nspecies; ++n) {
c(i,j,k,n) = 0.;
}
Real x,y,z;

for(int i1=0; i1<nsub; ++i1) {
for(int j1=0; j1<nsub; ++j1) {

AMREX_D_TERM(x = prob_lo[0] + i*dx[0] + (i1+0.5)*dxsub ;,
y = prob_lo[1] + j*dx[1] + (j1+0.5)*dysub ;,
z = prob_lo[2] + (k+0.5)*dx[2] - center[2];);

// Real height = film_thickness + 0.5*film_thickness*amrex::Math::cospi(2.*x/L[0]);
Real height = film_thickness + x;
//height = film_thickness ;

// amrex::Print{} << "x,thickness,height " << x << " " << film_thickness << " " << height << std::endl;

if (smoothing_width == 0.) {

// discontinuous interface
if (y < height) {
for (int n=0; n<nspecies; ++n) {
c(i,j,k,n) += c_init_1[n];
}
} else {
for (int n=0; n<nspecies; ++n) {
c(i,j,k,n) += c_init_2[n];
}
}

} else {
// smooth interface
for (int n=0; n<nspecies; ++n) {
c(i,j,k,n) += c_init_1[n] + (c_init_2[n]-c_init_1[n]) *
0.5*(1. + std::tanh((y-height)/(smoothing_width*dx[0])));
}
}
}
}
for (int n=0; n<nspecies; ++n) {
c(i,j,k,n) = c(i,j,k,n)/(factor*factor);
}
});



} else if (prob_type == 3) {
Expand Down Expand Up @@ -654,7 +726,7 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac,
for (int n=0; n<nspecies; ++n) {
c(i,j,k,n) = c_init_1[n] +
(c_init_2[n] - c_init_1[n])*
(1./(1.+std::exp(-smoothing_width*(y-l1))) - 1./(1.+std::exp(-smoothing_width*(y-l2))));
(1./(1.+std::exp(-(y-l1)/(smoothing_width*dx[0]))) - 1./(1.+std::exp(-(y-l2)/(smoothing_width*dx[0]))));
}
}
});
Expand Down
39 changes: 37 additions & 2 deletions src_common/BCPhysToMath.cpp
View file
Open in desktop
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] = amrex::BCType::foextrap;
}
else if (bc_mass_lo[i] == 2) {
// reservoir -> dirichlet
bc_lo[i] = amrex::BCType::ext_dir;
}
else if (bc_mass_lo[i] == 4) {
// reservoir -> dirichlet
bc_lo[i] = amrex::BCType::foextrap;
}
if (bc_mass_hi[i] == 1) {
// wall -> first-order extrapolation
bc_hi[i] = amrex::BCType::foextrap;
}
else if (bc_mass_hi[i] == 2) {
// reservoir -> dirichlet
bc_hi[i] = amrex::BCType::ext_dir;
}
else if (bc_mass_hi[i] == 4) {
// reservoir -> dirichlet
bc_hi[i] = amrex::BCType::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 All @@ -43,6 +70,10 @@ void BCPhysToMath(int bccomp, amrex::Vector<int>& bc_lo, amrex::Vector<int>& bc_
// reservoir -> dirichlet
bc_lo[i] = amrex::BCType::ext_dir;
}
else if (bc_mass_lo[i] == 4) {
// reservoir -> dirichlet
bc_lo[i] = SPEC_CONTACT_BC;
}
if (bc_mass_hi[i] == 1) {
// wall -> first-order extrapolation
bc_hi[i] = amrex::BCType::foextrap;
Expand All @@ -51,6 +82,10 @@ void BCPhysToMath(int bccomp, amrex::Vector<int>& bc_lo, amrex::Vector<int>& bc_
// reservoir -> dirichlet
bc_hi[i] = amrex::BCType::ext_dir;
}
else if (bc_mass_hi[i] == 4) {
// reservoir -> dirichlet
bc_hi[i] = SPEC_CONTACT_BC;
}
}
}
else if (bccomp == TEMP_BC_COMP) { // TEMPERATURE
Expand Down
12 changes: 6 additions & 6 deletions src_common/ComputeDivAndGrad.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ void ComputeGrad(const MultiFab & phi_in, std::array<MultiFab, AMREX_SPACEDIM> &
// boundary conditions
// note: at physical boundaries,
// alter stencil at boundary since ghost value represents value at boundary
if (bc_lo[0] == amrex::BCType::foextrap || bc_lo[0] == amrex::BCType::ext_dir) {
if (bc_lo[0] == amrex::BCType::foextrap || bc_lo[0] == amrex::BCType::ext_dir || bc_lo[0] == SPEC_CONTACT_BC) {
if (bx_x.smallEnd(0) <= dom.smallEnd(0)) {
int lo = dom.smallEnd(0);
amrex::ParallelFor(bx_x, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -106,7 +106,7 @@ void ComputeGrad(const MultiFab & phi_in, std::array<MultiFab, AMREX_SPACEDIM> &
}
}

if (bc_hi[0] == amrex::BCType::foextrap || bc_hi[0] == amrex::BCType::ext_dir) {
if (bc_hi[0] == amrex::BCType::foextrap || bc_hi[0] == amrex::BCType::ext_dir || bc_hi[0] == SPEC_CONTACT_BC) {
if (bx_x.bigEnd(0) >= dom.bigEnd(0)+1) {
int hi = dom.bigEnd(0)+1;
amrex::ParallelFor(bx_x, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -120,7 +120,7 @@ void ComputeGrad(const MultiFab & phi_in, std::array<MultiFab, AMREX_SPACEDIM> &
}
}

if (bc_lo[1] == amrex::BCType::foextrap || bc_lo[1] == amrex::BCType::ext_dir) {
if (bc_lo[1] == amrex::BCType::foextrap || bc_lo[1] == amrex::BCType::ext_dir || bc_lo[1] == SPEC_CONTACT_BC) {
if (bx_y.smallEnd(1) <= dom.smallEnd(1)) {
int lo = dom.smallEnd(1);
amrex::ParallelFor(bx_y, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -134,7 +134,7 @@ void ComputeGrad(const MultiFab & phi_in, std::array<MultiFab, AMREX_SPACEDIM> &
}
}

if (bc_hi[1] == amrex::BCType::foextrap || bc_hi[1] == amrex::BCType::ext_dir) {
if (bc_hi[1] == amrex::BCType::foextrap || bc_hi[1] == amrex::BCType::ext_dir || bc_hi[1] == SPEC_CONTACT_BC) {
if (bx_y.bigEnd(1) >= dom.bigEnd(1)+1) {
int hi = dom.bigEnd(1)+1;
amrex::ParallelFor(bx_y, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -149,7 +149,7 @@ void ComputeGrad(const MultiFab & phi_in, std::array<MultiFab, AMREX_SPACEDIM> &
}

#if (AMREX_SPACEDIM == 3)
if (bc_lo[2] == amrex::BCType::foextrap || bc_lo[2] == amrex::BCType::ext_dir) {
if (bc_lo[2] == amrex::BCType::foextrap || bc_lo[2] == amrex::BCType::ext_dir || bc_lo[2] == SPEC_CONTACT_BC) {
if (bx_z.smallEnd(2) <= dom.smallEnd(2)) {
int lo = dom.smallEnd(2);
amrex::ParallelFor(bx_z, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -163,7 +163,7 @@ void ComputeGrad(const MultiFab & phi_in, std::array<MultiFab, AMREX_SPACEDIM> &
}
}

if (bc_hi[2] == amrex::BCType::foextrap || bc_hi[2] == amrex::BCType::ext_dir) {
if (bc_hi[2] == amrex::BCType::foextrap || bc_hi[2] == amrex::BCType::ext_dir || bc_hi[2] == SPEC_CONTACT_BC) {
if (bx_z.bigEnd(2) >= dom.bigEnd(2)+1) {
int hi = dom.bigEnd(2)+1;
amrex::ParallelFor(bx_z, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand Down
12 changes: 6 additions & 6 deletions src_common/ConvertStag.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ void AverageCCToFace(const MultiFab& cc_in, std::array<MultiFab, AMREX_SPACEDIM>
// note: at physical boundaries,
// the value in the ghost cells represent the value ON the boundary
// so we simply copy the ghost cell value into the value on the domain (and ghost faces too)
if (bc_lo[0] == amrex::BCType::foextrap || bc_lo[0] == amrex::BCType::ext_dir) {
if (bc_lo[0] == amrex::BCType::foextrap || bc_lo[0] == amrex::BCType::ext_dir || bc_lo[0] == SPEC_CONTACT_BC) {
if (bx_x.smallEnd(0) <= dom.smallEnd(0)) {
int lo = dom.smallEnd(0);
amrex::ParallelFor(bx_x, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -93,7 +93,7 @@ void AverageCCToFace(const MultiFab& cc_in, std::array<MultiFab, AMREX_SPACEDIM>
}
}

if (bc_hi[0] == amrex::BCType::foextrap || bc_hi[0] == amrex::BCType::ext_dir) {
if (bc_hi[0] == amrex::BCType::foextrap || bc_hi[0] == amrex::BCType::ext_dir || bc_hi[0] == SPEC_CONTACT_BC) {
if (bx_x.bigEnd(0) >= dom.bigEnd(0)+1) {
int hi = dom.bigEnd(0)+1;
amrex::ParallelFor(bx_x, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -105,7 +105,7 @@ void AverageCCToFace(const MultiFab& cc_in, std::array<MultiFab, AMREX_SPACEDIM>
}
}

if (bc_lo[1] == amrex::BCType::foextrap || bc_lo[1] == amrex::BCType::ext_dir) {
if (bc_lo[1] == amrex::BCType::foextrap || bc_lo[1] == amrex::BCType::ext_dir || bc_lo[1] == SPEC_CONTACT_BC) {
if (bx_y.smallEnd(1) <= dom.smallEnd(1)) {
int lo = dom.smallEnd(1);
amrex::ParallelFor(bx_y, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -117,7 +117,7 @@ void AverageCCToFace(const MultiFab& cc_in, std::array<MultiFab, AMREX_SPACEDIM>
}
}

if (bc_hi[1] == amrex::BCType::foextrap || bc_hi[1] == amrex::BCType::ext_dir) {
if (bc_hi[1] == amrex::BCType::foextrap || bc_hi[1] == amrex::BCType::ext_dir || bc_hi[1] == SPEC_CONTACT_BC) {
if (bx_y.bigEnd(1) >= dom.bigEnd(1)+1) {
int hi = dom.bigEnd(1)+1;
amrex::ParallelFor(bx_y, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -130,7 +130,7 @@ void AverageCCToFace(const MultiFab& cc_in, std::array<MultiFab, AMREX_SPACEDIM>
}

#if (AMREX_SPACEDIM == 3)
if (bc_lo[2] == amrex::BCType::foextrap || bc_lo[2] == amrex::BCType::ext_dir) {
if (bc_lo[2] == amrex::BCType::foextrap || bc_lo[2] == amrex::BCType::ext_dir || bc_lo[2] == SPEC_CONTACT_BC) {
if (bx_z.smallEnd(2) <= dom.smallEnd(2)) {
int lo = dom.smallEnd(2);
amrex::ParallelFor(bx_z, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -142,7 +142,7 @@ void AverageCCToFace(const MultiFab& cc_in, std::array<MultiFab, AMREX_SPACEDIM>
}
}

if (bc_hi[2] == amrex::BCType::foextrap || bc_hi[2] == amrex::BCType::ext_dir) {
if (bc_hi[2] == amrex::BCType::foextrap || bc_hi[2] == amrex::BCType::ext_dir || bc_hi[2] == SPEC_CONTACT_BC) {
if (bx_z.bigEnd(2) >= dom.bigEnd(2)+1) {
int hi = dom.bigEnd(2)+1;
amrex::ParallelFor(bx_z, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand Down
Loading
Loading