Skip to content

Commit

Permalink
Hari/features (#14)
Browse files Browse the repository at this point in the history 
1. adding conjugate heat transfer
2. splitting files for split chem evolve for readability
3. gas phase flux computation now restricted only to gas-gas and gas-mostly-gas
faces, increases performance a bit
  • Loading branch information
@hsitaram
hsitaram authored Aug 1, 2024
1 parent 9b0865e commit c0ceab4
Show file tree
Hide file tree
Showing 32 changed files with 1,195 additions and 343 deletions.
2 changes: 1 addition & 1 deletion Source/Make.package
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
CEXE_sources += main.cpp mflo.cpp mflo_evolve.cpp mflo_plot.cpp
CEXE_sources += mflo_gridwork.cpp mflo_timestep.cpp
CEXE_sources += mflo_gridwork.cpp mflo_timestep.cpp mflo_splitchem.cpp

CEXE_headers += mflo.H kernels_3d.H mflo_tagging.H mflo_bcfill.H globalDefines.H mflo_constants.H
190 changes: 107 additions & 83 deletions Source/Src_3d/compute_flux_3d.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,6 +188,13 @@ void compute_flux(int i, int j, int k,int sweepdir,
facevel[XDIR] = half*(phi(left,VELX_INDX) + phi(right,VELX_INDX));
facevel[YDIR] = half*(phi(left,VELY_INDX) + phi(right,VELY_INDX));
facevel[ZDIR] = half*(phi(left,VELZ_INDX) + phi(right,VELZ_INDX));

Real facevfrac;
facevfrac=half*(phi(left,VFRAC_INDX)+phi(right,VFRAC_INDX));

//note:
//if facevfrac >=0.5, this face is a gas cell-fractional cell interface
//if facevfrac < 0.5, this face is solid cell-fractional cell interface

for(int c=0;c<NCVARS;c++)
{
Expand Down Expand Up @@ -241,116 +248,133 @@ void compute_flux(int i, int j, int k,int sweepdir,
}
}

if(hyperbolics_order==1)
{
ausmp_up_flux(ul,ur,fhalf,normal);
}
else if(hyperbolics_order==2)
{
get_higherorder_states(ul2,ur2,ul,ur,ulm1,urp1);
ausmp_up_flux(ul2,ur2,fhalf,normal);
}
else if(hyperbolics_order==4)
{
central_with_diss_flux(ul,ur,ulm1,urp1,fhalf,normal,dissfactor);
}
else if(hyperbolics_order==5)
{
weno_flux(ulm2,ulm1,ul,ur,urp1,urp2,fhalf,normal);
}
else
//do hyperbolic fluxes only at
//gas-cell--gas-cell interface or
//gas-cell--fractional-cell interface
if(facevfrac >= 0.5)
{
amrex::Abort("unknown hyperbolics scheme");
}
if(hyperbolics_order==1)
{
ausmp_up_flux(ul,ur,fhalf,normal);
}
else if(hyperbolics_order==2)
{
get_higherorder_states(ul2,ur2,ul,ur,ulm1,urp1);
ausmp_up_flux(ul2,ur2,fhalf,normal);
}
else if(hyperbolics_order==4)
{
central_with_diss_flux(ul,ur,ulm1,urp1,fhalf,normal,dissfactor);
}
else if(hyperbolics_order==5)
{
weno_flux(ulm2,ulm1,ul,ur,urp1,urp2,fhalf,normal);
}
else
{
amrex::Abort("unknown hyperbolics scheme");
}

for(int c=0;c<NCVARS;c++)
{
fx(i,j,k,c+RHO_INDX) = fhalf[c+RHO_IND];
}
for(int sp=0;sp<NUM_SPECIES;sp++)
{
fx(i,j,k,sp+FLO_NVARS) = mflo_species::advect_flags[sp]*fhalf[NCVARS+sp];
}
for(int c=0;c<NCVARS;c++)
{
fx(i,j,k,c+RHO_INDX) = fhalf[c+RHO_IND];
}

gradu[0] = (phi(right,VELX_INDX) - phi(left,VELX_INDX))/dx[sweepdir];
//this is still correct in this if loop,
//as we dont have any advection in solid
for(int sp=0;sp<NUM_SPECIES;sp++)
{
fx(i,j,k,sp+FLO_NVARS) = mflo_species::advect_flags[sp]*fhalf[NCVARS+sp];
}

gradu[1] = ( half*(phi(top_left,VELX_INDX) + phi(top_right,VELX_INDX))
- half*(phi(bot_left,VELX_INDX) + phi(bot_right,VELX_INDX)) )/(two*dx[trans1dir]);
//viscous fluxes again inside this if loop
//only matters in gas-cell-gas-cell and
//gas-cell-fractional-cell interfaces
gradu[0] = (phi(right,VELX_INDX) - phi(left,VELX_INDX))/dx[sweepdir];

gradu[2] = ( half*(phi(frt_left,VELX_INDX) + phi(frt_right,VELX_INDX))
- half*(phi(bck_left,VELX_INDX) + phi(bck_right,VELX_INDX)) )/(two*dx[trans2dir]);
gradu[1] = ( half*(phi(top_left,VELX_INDX) + phi(top_right,VELX_INDX))
- half*(phi(bot_left,VELX_INDX) + phi(bot_right,VELX_INDX)) )/(two*dx[trans1dir]);

gradv[0] = (phi(right,VELY_INDX) - phi(left,VELY_INDX))/dx[sweepdir];
gradu[2] = ( half*(phi(frt_left,VELX_INDX) + phi(frt_right,VELX_INDX))
- half*(phi(bck_left,VELX_INDX) + phi(bck_right,VELX_INDX)) )/(two*dx[trans2dir]);

gradv[1] = ( half*(phi(top_left,VELY_INDX) + phi(top_right,VELY_INDX))
- half*(phi(bot_left,VELY_INDX) + phi(bot_right,VELY_INDX)) )/(two*dx[trans1dir]);
gradv[0] = (phi(right,VELY_INDX) - phi(left,VELY_INDX))/dx[sweepdir];

gradv[2] = ( half*(phi(frt_left,VELY_INDX) + phi(frt_right,VELY_INDX))
- half*(phi(bck_left,VELY_INDX) + phi(bck_right,VELY_INDX)) )/(two*dx[trans2dir]);
gradv[1] = ( half*(phi(top_left,VELY_INDX) + phi(top_right,VELY_INDX))
- half*(phi(bot_left,VELY_INDX) + phi(bot_right,VELY_INDX)) )/(two*dx[trans1dir]);

gradw[0] = (phi(right,VELZ_INDX) - phi(left,VELZ_INDX))/dx[sweepdir];
gradv[2] = ( half*(phi(frt_left,VELY_INDX) + phi(frt_right,VELY_INDX))
- half*(phi(bck_left,VELY_INDX) + phi(bck_right,VELY_INDX)) )/(two*dx[trans2dir]);

gradw[1] = ( half*(phi(top_left,VELZ_INDX) + phi(top_right,VELZ_INDX))
- half*(phi(bot_left,VELZ_INDX) + phi(bot_right,VELZ_INDX)) )/(two*dx[trans1dir]);
gradw[0] = (phi(right,VELZ_INDX) - phi(left,VELZ_INDX))/dx[sweepdir];

gradw[2] = ( half*(phi(frt_left,VELZ_INDX) + phi(frt_right,VELZ_INDX))
- half*(phi(bck_left,VELZ_INDX) + phi(bck_right,VELZ_INDX)) )/(two*dx[trans2dir]);
gradw[1] = ( half*(phi(top_left,VELZ_INDX) + phi(top_right,VELZ_INDX))
- half*(phi(bot_left,VELZ_INDX) + phi(bot_right,VELZ_INDX)) )/(two*dx[trans1dir]);

dudx=gradu[GET_XDIR(sweepdir)];
dudy=gradu[GET_YDIR(sweepdir)];
dudz=gradu[GET_ZDIR(sweepdir)];
gradw[2] = ( half*(phi(frt_left,VELZ_INDX) + phi(frt_right,VELZ_INDX))
- half*(phi(bck_left,VELZ_INDX) + phi(bck_right,VELZ_INDX)) )/(two*dx[trans2dir]);

dvdx=gradv[GET_XDIR(sweepdir)];
dvdy=gradv[GET_YDIR(sweepdir)];
dvdz=gradv[GET_ZDIR(sweepdir)];
dudx=gradu[GET_XDIR(sweepdir)];
dudy=gradu[GET_YDIR(sweepdir)];
dudz=gradu[GET_ZDIR(sweepdir)];

dwdx=gradw[GET_XDIR(sweepdir)];
dwdy=gradw[GET_YDIR(sweepdir)];
dwdz=gradw[GET_ZDIR(sweepdir)];
dvdx=gradv[GET_XDIR(sweepdir)];
dvdy=gradv[GET_YDIR(sweepdir)];
dvdz=gradv[GET_ZDIR(sweepdir)];

dTdn = (phi(right,TEMP_INDX) - phi(left,TEMP_INDX))/dx[sweepdir];
dwdx=gradw[GET_XDIR(sweepdir)];
dwdy=gradw[GET_YDIR(sweepdir)];
dwdz=gradw[GET_ZDIR(sweepdir)];

delu = dudx + dvdy + dwdz;

tau[XDIR][XDIR]= visc*(two*dudx - two3rd*delu);
tau[XDIR][YDIR]= visc*(dudy + dvdx);
tau[XDIR][ZDIR]= visc*(dudz + dwdx);
delu = dudx + dvdy + dwdz;

tau[YDIR][XDIR]= visc*(dvdx + dudy);
tau[YDIR][YDIR]= visc*(two*dvdy - two3rd*delu);
tau[YDIR][ZDIR]= visc*(dvdz + dwdy);
tau[XDIR][XDIR]= visc*(two*dudx - two3rd*delu);
tau[XDIR][YDIR]= visc*(dudy + dvdx);
tau[XDIR][ZDIR]= visc*(dudz + dwdx);

tau[ZDIR][XDIR]= visc*(dwdx + dudz);
tau[ZDIR][YDIR]= visc*(dwdy + dvdz);
tau[ZDIR][ZDIR]= visc*(two*dwdz - two3rd*delu);
tau[YDIR][XDIR]= visc*(dvdx + dudy);
tau[YDIR][YDIR]= visc*(two*dvdy - two3rd*delu);
tau[YDIR][ZDIR]= visc*(dvdz + dwdy);

fx(i,j,k,RHOU_INDX) -= ( tau[XDIR][XDIR]*normal[XDIR]
+ tau[YDIR][XDIR]*normal[YDIR]
+ tau[ZDIR][XDIR]*normal[ZDIR]);
tau[ZDIR][XDIR]= visc*(dwdx + dudz);
tau[ZDIR][YDIR]= visc*(dwdy + dvdz);
tau[ZDIR][ZDIR]= visc*(two*dwdz - two3rd*delu);

fx(i,j,k,RHOV_INDX) -= ( tau[XDIR][YDIR]*normal[XDIR]
+ tau[YDIR][YDIR]*normal[YDIR]
+ tau[ZDIR][YDIR]*normal[ZDIR]);
fx(i,j,k,RHOU_INDX) -= ( tau[XDIR][XDIR]*normal[XDIR]
+ tau[YDIR][XDIR]*normal[YDIR]
+ tau[ZDIR][XDIR]*normal[ZDIR]);

fx(i,j,k,RHOW_INDX) -= ( tau[XDIR][ZDIR]*normal[XDIR]
+ tau[YDIR][ZDIR]*normal[YDIR]
+ tau[ZDIR][ZDIR]*normal[ZDIR]);
fx(i,j,k,RHOV_INDX) -= ( tau[XDIR][YDIR]*normal[XDIR]
+ tau[YDIR][YDIR]*normal[YDIR]
+ tau[ZDIR][YDIR]*normal[ZDIR]);

fx(i,j,k,RHOE_INDX) -= thcond*dTdn;
fx(i,j,k,RHOW_INDX) -= ( tau[XDIR][ZDIR]*normal[XDIR]
+ tau[YDIR][ZDIR]*normal[YDIR]
+ tau[ZDIR][ZDIR]*normal[ZDIR]);

//this term is the work done on the control volume by viscous forces - tau_{ji} n_j v_i
fx(i,j,k,RHOE_INDX) -= (tau[XDIR][XDIR]*normal[XDIR] +
tau[YDIR][XDIR]*normal[YDIR] +
tau[ZDIR][XDIR]*normal[ZDIR] )*facevel[XDIR];

fx(i,j,k,RHOE_INDX) -= (tau[XDIR][YDIR]*normal[XDIR] +
tau[YDIR][YDIR]*normal[YDIR] +
tau[ZDIR][YDIR]*normal[ZDIR] )*facevel[YDIR];

//this term is the work done on the control volume by viscous forces - tau_{ji} n_j v_i
fx(i,j,k,RHOE_INDX) -= (tau[XDIR][XDIR]*normal[XDIR] +
tau[YDIR][XDIR]*normal[YDIR] +
tau[ZDIR][XDIR]*normal[ZDIR] )*facevel[XDIR];
fx(i,j,k,RHOE_INDX) -= (tau[XDIR][ZDIR]*normal[XDIR] +
tau[YDIR][ZDIR]*normal[YDIR] +
tau[ZDIR][ZDIR]*normal[ZDIR] )*facevel[ZDIR];

fx(i,j,k,RHOE_INDX) -= (tau[XDIR][YDIR]*normal[XDIR] +
tau[YDIR][YDIR]*normal[YDIR] +
tau[ZDIR][YDIR]*normal[ZDIR] )*facevel[YDIR];
}

//keeping thermal conduction outside
//to include conjugate heat transfer
dTdn = (phi(right,TEMP_INDX) - phi(left,TEMP_INDX))/dx[sweepdir];
fx(i,j,k,RHOE_INDX) -= thcond*dTdn;

fx(i,j,k,RHOE_INDX) -= (tau[XDIR][ZDIR]*normal[XDIR] +
tau[YDIR][ZDIR]*normal[YDIR] +
tau[ZDIR][ZDIR]*normal[ZDIR] )*facevel[ZDIR];
//amrex::Print()<<"dTdn:"<<dTdn<<"\t"<<thcond<<"\t"<<fx(i,j,k,RHOE_INDX)<<"\t"
//<<phi(iv,VFRAC_INDX)<<"\t"<<facevel[0]<<"\t"<<facevel[1]<<"\t"<<facevel[2]<<"\n";
}
else
{
Expand Down
41 changes: 34 additions & 7 deletions Source/Src_3d/div_3d.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,8 @@ void update_residual(
Array4<Real> const& flxy,
Array4<Real> const& flxz),
const GpuArray<Real, AMREX_SPACEDIM>& dx,
int spec_in_solid)
int spec_in_solid,
int conj_ht)
{
// remember, we are solve dudt + del.F = 0
dsdt(i, j, k, n) = (flxx(i, j, k, n) - flxx(i + 1, j, k, n)) / dx[0] +
Expand All @@ -33,15 +34,41 @@ void update_residual(

int nvars_to_nullify=(spec_in_solid==1)?FLO_NVARS:FLO_NVARS+NUM_SPECIES;


if(n < nvars_to_nullify)
//check for rhoe because we may have conjugate heat transfer
if(n!=RHOE_INDX)
{
if(n < nvars_to_nullify)
{
if(phi(i,j,k,VFRAC_INDX) < one)
{
dsdt(i, j, k, n)=zeroval;
}
}
}
//if n is rhoe
else
{
if(phi(i,j,k,VFRAC_INDX) < one)
{
dsdt(i, j, k, n)=zeroval;
//if no conjugate heat transfer
//null the dsdt
if(!conj_ht)
{
if(phi(i,j,k,VFRAC_INDX) < one)
{
dsdt(i, j, k, n)=zeroval;
}
}
else
{
//when conjugate heat transfer is there
//null the dsdt in fractional cells
//as these are recalculated from cut cells
if(phi(i,j,k,VFRAC_INDX) > zeroval && phi(i,j,k,VFRAC_INDX) < one)
{
dsdt(i, j, k, n)=zeroval;
}
}
}

if(dsdt(i,j,k,n)!=dsdt(i,j,k,n))
{
amrex::Abort("NaN detected in dsdt\n");
Expand Down
Loading

0 comments on commit c0ceab4

Please sign in to comment.