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

Hari/features #14

Merged
merged 6 commits into from
Aug 1, 2024
Merged

Hari/features #14

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
0734211
moving split chem into different file
hsitaram Jul 29, 2024
1d9e288
adding split chem file
hsitaram Jul 29, 2024
666b8bf
added conj heat transfer, hasnt broken anything
hsitaram Jul 30, 2024
4e9d3c3
adding conj heat transfer test
hsitaram Jul 30, 2024
f5f2394
updating all cases with new thermo functions and testing
hsitaram Aug 1, 2024
c49f8a3
removing print statement in div to keep gpu build happy
hsitaram Aug 1, 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
2 changes: 1 addition & 1 deletion Source/Make.package
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading