Skip to content

Commit

Permalink
Hari/dielectrics (#7)
Browse files Browse the repository at this point in the history
* adding dielectric surface charging and
dielectric layer boundary conditions

---------

Co-authored-by: Hariswaran Sitaraman <[email protected]>
Co-authored-by: Hariswaran Sitaraman <[email protected]>
  • Loading branch information
3 people authored Nov 1, 2024
1 parent aa5b679 commit f0cf6eb
Show file tree
Hide file tree
Showing 39 changed files with 2,483 additions and 225 deletions.
1 change: 0 additions & 1 deletion Source/ElecEnergySource.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,6 @@ void Vidyut::compute_elecenergy_source(int lev,
efield[1].array(mfi),
efield[2].array(mfi))};

// update residual
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {

//Joule heating
Expand Down
17 changes: 11 additions & 6 deletions Source/Evolve.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@ void Vidyut::Evolve()
Vector<MultiFab> phi_tmp(finest_level+1);

// edge centered efield
Vector< Array<MultiFab,AMREX_SPACEDIM> > efield_ec(finest_level+1);
Vector< Array<MultiFab,AMREX_SPACEDIM> > efield_fc(finest_level+1);

//copy new to old and update time
for(int lev=0;lev<=finest_level;lev++)
Expand Down Expand Up @@ -134,8 +134,8 @@ void Vidyut::Evolve()
flux[lev][idim].define(ba, dmap[lev], 1, 0);
flux[lev][idim].setVal(0.0);

efield_ec[lev][idim].define(ba, dmap[lev], 1, 0);
efield_ec[lev][idim].setVal(0.0);
efield_fc[lev][idim].define(ba, dmap[lev], 1, 0);
efield_fc[lev][idim].setVal(0.0);

gradne_fc[lev][idim].define(ba, dmap[lev], 1, 0);
gradne_fc[lev][idim].setVal(0.0);
Expand Down Expand Up @@ -168,13 +168,13 @@ void Vidyut::Evolve()
flux[lev][idim].setVal(0.0);
gradne_fc[lev][idim].setVal(0.0);
grad_fc[lev][idim].setVal(0.0);
efield_ec[lev][idim].setVal(0.0);
efield_fc[lev][idim].setVal(0.0);
}
expl_src[lev].setVal(0.0);
rxn_src[lev].setVal(0.0);
}

solve_potential(cur_time, Sborder, pot_bc_lo, pot_bc_hi, efield_ec);
solve_potential(cur_time, Sborder, pot_bc_lo, pot_bc_hi, efield_fc);

if(cs_technique)
{
Expand Down Expand Up @@ -237,7 +237,7 @@ void Vidyut::Evolve()
{
compute_elecenergy_source(lev, Sborder[lev],
rxn_src[lev],
efield_ec[lev],
efield_fc[lev],
gradne_fc[lev],
expl_src[lev], cur_time, dt_common,floor_jh);
}
Expand Down Expand Up @@ -302,6 +302,11 @@ void Vidyut::Evolve()
}
}
}

if(track_surf_charge)
{
update_surf_charge(Sborder,cur_time,dt_common);
}
}

if(niter<num_timestep_correctors-1)
Expand Down
1 change: 1 addition & 0 deletions Source/HelperFuncs.H
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@

using namespace amrex;


AMREX_GPU_HOST_DEVICE
AMREX_FORCE_INLINE
amrex::Real get_efield_alongdir(int i,int j,int k,int dir,
Expand Down
2 changes: 1 addition & 1 deletion Source/Make.package
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
CEXE_sources += main.cpp Vidyut.cpp Evolve.cpp Plot.cpp
CEXE_sources += GridWork.cpp TimeStep.cpp PotentialSolve.cpp ScalarSolve.cpp
CEXE_sources += ElecEnergySource.cpp PlasmaChem.cpp
CEXE_sources += ElecEnergySource.cpp PlasmaChem.cpp SurfCharge.cpp
CEXE_headers += Vidyut.H Tagging.H BCfill.H Constants.H UnivConstants.H VarDefines.H
CEXE_headers += compute_explicit_flux.H PlasmaChem.H HelperFuncs.H
137 changes: 137 additions & 0 deletions Source/SurfCharge.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,137 @@
#include <Vidyut.H>
#include <UserSources.H>
#include <Chemistry.H>
#include <UserFunctions.H>
#include <PlasmaChem.H>
#include <ProbParm.H>
#include <AMReX_MLABecLaplacian.H>
#include <HelperFuncs.H>

AMREX_GPU_HOST_DEVICE
AMREX_FORCE_INLINE
amrex::Real charge_flux(int i,int j,int k, int sgn, int dir,
const GpuArray<Real,AMREX_SPACEDIM> dx,
amrex::Real gastemp,
Array4<Real> const& sb_arr)
{
IntVect icell{AMREX_D_DECL(i,j,k)};
IntVect icell_prvs{AMREX_D_DECL(i,j,k)};
amrex::Real outward_normal[AMREX_SPACEDIM]={0.0};
outward_normal[dir]=sgn;
int cell_adjust;

cell_adjust = (sgn==1)?-1:0;
icell[dir]+=cell_adjust;
//i+1 at domain lo and i-1 at
//domain hi
icell_prvs[dir]=icell[dir]-sgn;


amrex::Real etemp=sb_arr(icell,ETEMP_ID);
amrex::Real ebyn=sb_arr(icell,REF_ID);
amrex::Real efield_mag=sb_arr(icell,EFX_ID)*sb_arr(icell,EFX_ID);
#if AMREX_SPACEDIM > 1
efield_mag += sb_arr(icell,EFY_ID)*sb_arr(icell,EFY_ID);
#if AMREX_SPACEDIM == 3
efield_mag += sb_arr(icell,EFZ_ID)*sb_arr(icell,EFZ_ID);
#endif
#endif
efield_mag=std::sqrt(efield_mag);
amrex::Real q_times_flux=0.0;
amrex::Real efield_n=sb_arr(icell,EFX_ID+dir)*outward_normal[dir];
amrex::Real numdens=0.0;
for(int sp=0; sp<NUM_SPECIES; sp++) numdens += sb_arr(icell,sp);
for(int sp=0;sp<NUM_SPECIES;sp++)
{
amrex::Real chrg=plasmachem::get_charge(sp);

if(amrex::Math::abs(chrg) > 0)
{
amrex::Real specdens=sb_arr(icell,sp);
amrex::Real gradn_n=(sb_arr(icell,sp)-sb_arr(icell_prvs,sp))/dx[dir];
amrex::Real mu = specMob(sp, etemp, numdens, efield_mag, gastemp);
amrex::Real dcoeff = specDiff(sp, etemp, numdens, efield_mag, gastemp);

amrex::Real flx=mu*specdens*efield_n-dcoeff*gradn_n;

//add only when flux is directed towards the surface
if(flx > 0.0)
{
q_times_flux += chrg*flx;
}
}
}

q_times_flux*=ECHARGE;
return(q_times_flux);
}

void Vidyut::update_surf_charge(Vector<MultiFab>& Sborder,
Real current_time, Real dt)
{
amrex::Real gastemp=gas_temperature;
amrex::Real tstep=dt;
ProbParm const* localprobparm = d_prob_parm;
for (int ilev = 0; ilev <= finest_level; ilev++)
{
// set boundary conditions
for (MFIter mfi(phi_new[ilev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
const auto dx = geom[ilev].CellSizeArray();
auto prob_lo = geom[ilev].ProbLoArray();
auto prob_hi = geom[ilev].ProbHiArray();
const Box& domain = geom[ilev].Domain();

Array4<Real> sb_arr = Sborder[ilev].array(mfi);
Array4<Real> phi_arr = phi_new[ilev].array(mfi);
Real time = current_time; // for GPU capture

for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
{
//note: bdryLo/bdryHi grabs the face indices from bx that are the boundary
//since they are face indices, the bdry normal index is 0/n+1, n is number of cells
//so the ghost cell index at left side is i-1 while it is i on the right
if (bx.smallEnd(idim) == domain.smallEnd(idim))
{
int sign=-1;
amrex::ParallelFor(amrex::bdryLo(bx, idim),
[=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {

int dielectricflag=user_transport::is_dielectric(i, j, k, idim, sign,
prob_lo, prob_hi, dx,
time,*localprobparm);
if(dielectricflag)
{
IntVect icell{AMREX_D_DECL(i,j,k)};
int cell_adjust = (sign==1)?-1:0;
icell[idim]+=cell_adjust;
amrex::Real q_times_flux=charge_flux(i,j,k,sign,idim,
dx,gastemp,sb_arr);
phi_arr(icell,SRFCH_ID)+= q_times_flux*tstep;
}
});
}
if (bx.bigEnd(idim) == domain.bigEnd(idim))
{
int sign=1;
amrex::ParallelFor(amrex::bdryHi(bx, idim),
[=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {

int dielectricflag=user_transport::is_dielectric(i, j, k, idim, sign,
prob_lo, prob_hi, dx,
time,*localprobparm);
if(dielectricflag)
{
IntVect icell{AMREX_D_DECL(i,j,k)};
int cell_adjust = (sign==1)?-1:0;
icell[idim]+=cell_adjust;
amrex::Real q_times_flux=charge_flux(i,j,k,sign,idim,dx,gastemp,sb_arr);
phi_arr(icell,SRFCH_ID)+= q_times_flux*tstep;
}
});
}
}
}
}
}
3 changes: 2 additions & 1 deletion Source/VarDefines.H
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,8 @@
#define EIH_ID (NUM_SPECIES+7) //electron inelastic heating
#define EEH_ID (NUM_SPECIES+8) //electron elastic heating
#define REF_ID (NUM_SPECIES+9)
#define NUM_PLASMAVARS (REF_ID-NUM_SPECIES+1)
#define SRFCH_ID (NUM_SPECIES+10)
#define NUM_PLASMAVARS (SRFCH_ID-NUM_SPECIES+1)
#define NVAR (NUM_SPECIES+NUM_PLASMAVARS)

#define PERBC 0
Expand Down
4 changes: 4 additions & 0 deletions Source/Vidyut.H
Original file line number Diff line number Diff line change
Expand Up @@ -206,6 +206,7 @@ public:
int hyp_order=1;
int ngrow_for_fillpatch=2;
int weno_scheme=3; //WENO-Z by default
int track_surf_charge=0;

// Include physics
int elecenergy_solve=0;
Expand Down Expand Up @@ -241,6 +242,9 @@ public:
void update_cs_technique_potential();
void potential_gradlimiter(amrex::Vector<MultiFab>& Sborder);

void update_surf_charge(Vector<MultiFab>& Sborder,
Real current_time, Real dt);


// how often each level regrids the higher levels of refinement
// (after a level advances that many time steps)
Expand Down
2 changes: 2 additions & 0 deletions Source/Vidyut.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ Vidyut::Vidyut()
plasma_param_names[7]="Electron_inelasticHeat";
plasma_param_names[8]="Electron_elasticHeat";
plasma_param_names[9]="ReducedEF";
plasma_param_names[10]="SurfaceCharge";

allvarnames.resize(NVAR);
for (int i = 0; i < NUM_SPECIES; i++)
Expand Down Expand Up @@ -325,6 +326,7 @@ void Vidyut::ReadParameters()
pp.queryarr("bg_species_ids",bg_specid_list);

pp.query("weno_scheme",weno_scheme);
pp.query("track_surf_charge",track_surf_charge);

if(hyp_order==1) //first order upwind
{
Expand Down
17 changes: 14 additions & 3 deletions test/Reactor/CO2_H2_Ar_mechanism_reduced/UserFunctions.H
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,17 @@
using namespace amrex;
namespace user_transport
{
AMREX_GPU_HOST_DEVICE
AMREX_FORCE_INLINE
int is_dielectric(int i, int j, int k, int dir, int sign,
const GpuArray<Real,AMREX_SPACEDIM> prob_lo,
const GpuArray<Real,AMREX_SPACEDIM> prob_hi,
const GpuArray<Real,AMREX_SPACEDIM> dx,
Real time)
{
return(0);
}

AMREX_GPU_DEVICE AMREX_INLINE
void potential_bc(int i, int j, int k,
int dir, int sgn,
Expand All @@ -28,7 +39,7 @@ namespace user_transport
IntVect ghost_cell(i,j,k);
amrex::Real outward_normal[AMREX_SPACEDIM]={0.0};
outward_normal[dir]=sgn;

int gcell_adjust,cell_adjust;

//ghost_cell is one behind
Expand Down Expand Up @@ -85,7 +96,7 @@ namespace user_transport
IntVect ghost_cell(i,j,k);
amrex::Real outward_normal[AMREX_SPACEDIM]={0.0};
outward_normal[dir]=sgn;

int gcell_adjust,cell_adjust;

//ghost_cell is one behind
Expand All @@ -95,7 +106,7 @@ namespace user_transport

ghost_cell[dir]+=gcell_adjust;
cell_int[dir]+=cell_adjust;

//default
robin_a(ghost_cell) = 0.0;
robin_b(ghost_cell) = 1.0;
Expand Down
2 changes: 1 addition & 1 deletion test/clean.sh
Original file line number Diff line number Diff line change
@@ -1 +1 @@
rm -rf plt* Backtrace* core.*
rm -rf plt* Backtrace* core.* chk*
22 changes: 17 additions & 5 deletions test/models/Ar_DC_1d/UserFunctions.H
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,18 @@
using namespace amrex;
namespace user_transport
{
AMREX_GPU_HOST_DEVICE
AMREX_FORCE_INLINE
int is_dielectric(int i, int j, int k, int dir, int sign,
const GpuArray<Real,AMREX_SPACEDIM> prob_lo,
const GpuArray<Real,AMREX_SPACEDIM> prob_hi,
const GpuArray<Real,AMREX_SPACEDIM> dx,
Real time,
ProbParm const& prob_parm)
{
return(0);
}

AMREX_GPU_DEVICE AMREX_INLINE
void potential_bc(int i, int j, int k,
int dir, int sgn,
Expand All @@ -28,7 +40,7 @@ namespace user_transport
IntVect ghost_cell{AMREX_D_DECL(i,j,k)};
amrex::Real outward_normal[AMREX_SPACEDIM]={0.0};
outward_normal[dir]=sgn;

int gcell_adjust,cell_adjust;

//ghost_cell is one behind
Expand Down Expand Up @@ -81,7 +93,7 @@ namespace user_transport
IntVect ghost_cell{AMREX_D_DECL(i,j,k)};
amrex::Real outward_normal[AMREX_SPACEDIM]={0.0};
outward_normal[dir]=Real(sgn);

int gcell_adjust,cell_adjust;

//ghost_cell is one behind
Expand All @@ -91,12 +103,12 @@ namespace user_transport

ghost_cell[dir]+=gcell_adjust;
cell_int[dir]+=cell_adjust;

//default
robin_a(ghost_cell) = 0.0;
robin_b(ghost_cell) = 1.0;
robin_f(ghost_cell) = 0.0;

amrex::Real ndens = 0.0;
for(int sp=0; sp<NUM_SPECIES; sp++) ndens += phi(cell_int,sp);

Expand Down Expand Up @@ -188,7 +200,7 @@ namespace user_transport
amrex::Real efieldmag=std::sqrt( std::pow(phi(cellid,EFX_ID),2.0)+
std::pow(phi(cellid,EFY_ID),2.0)+
std::pow(phi(cellid,EFZ_ID),2.0) );

amrex::Real ndens = 0.0;
for(int sp=0; sp<NUM_SPECIES; sp++) ndens += phi(cellid,sp);

Expand Down
Loading

0 comments on commit f0cf6eb

Please sign in to comment.