Commit 365967e0 authored by Udo Ziegler's avatar Udo Ziegler
Browse files

HLLD low Mach pressure correction term added

parent d450be87
......@@ -34,6 +34,7 @@
#define REACTIONS_TOL 1.e-5 /* REACTIONS SOLVER: REL. ERR TOL */
#define REACTIONS_ATOL_X 0. /* REACTIONS SOLVER: ABS. ERR TOL_X */
#define REACTIONS_ATOL_T 0.01 /* REACTIONS SOLVER: ABS. ERR TOL_T */
#define HLLD_PRESSURE_CORR YES /* HLLD: LOW MACH PRESSURE CORRECTION */
/* PHYSICS-RELATED PARAMETERS ----------------------------------------------- */
......@@ -42,8 +43,8 @@
/* NCCM-RELATED PARAMETERS -------------------------------------------------- */
#define CR_IH 1.75e-16 /* H CR IONISATION RATE STD: 2.5e-17,
Indriolo & McCall: 1.75e-16) */
#define CR_IH 1.75e-16 /* H CR IONISATION RATE STD: 2.5e-17
Indriolo & McCall: 1.75e-16 */
#define TDUST 10. /* DUST TEMPERATURE */
#define ADUST 1. /* DUST ABUNDANCE REL. TO MILKY WAY */
#define ISRF_G 0.03 /* ISRF HABING STRENGTH (STD=1.13) */
......
......@@ -72,13 +72,15 @@ void FVcomputeInterfaceState(GRD *g, flag_t coord_dir, double *hz, double **uL,
*
* @author Oliver Gressel, QMUL (adopted from the implementation by S. Fromang
* & P. Hennebelle)
* - rewritten by Udo Ziegler, AIP
* - algorithm rewritten by Udo Ziegler, AIP
* - low Mach number pressure correction term added by Udo Ziegler, AIP
* (taken from Minoshima+ 21, JCP)
*/
{
double sg,b1,rcL,rcR,rcs,sLs,sRs,sM,ptotS,e_den,e_nom,srhoLs,srhoRs,u2S,u3S,
b2S,b3S,rhoL,u1L,u2L,u3L,b2L,b3L,ptotL,rhoR,u1R,u2R,u3R,b2R,b3R,
ptotR,rhoLs,u2Ls,u3Ls,b2Ls,b3Ls,vdotbLs,rhoRs,u2Rs,u3Rs,b2Rs,b3Rs,
vdotbRs,uLs,uRs,pi,smax;
vdotbRs,uLs,uRs,pi,smax,cuL,cuR,phi;
int i,is,ie,i1,i2,i3,i4,i5,i6,iBx,iBy,iBz,ipr,iu;
......@@ -153,7 +155,7 @@ for(i=is; i<=ie; i++)
if(_C.mf)
{
b1=0.5*(uL[i4][i]+uR[i4][i]);
if(fabs(b1)<DEPS*sqrt((rhoL+rhoR)/pi)*smax) b1=0.;
// if(fabs(b1)<DEPS*sqrt((rhoL+rhoR)/pi)*smax) b1=0.;
b2L=uL[i5][i];
b2R=uR[i5][i];
......@@ -179,22 +181,33 @@ for(i=is; i<=ie; i++)
rhoRs=rcR/(sR[i]-sM);
srhoRs=sqrt(rhoRs/pi);
ptotS=(rcR*ptotL+rcL*ptotR-rcL*rcR*(u1R-u1L))/rcs;
/*ptotS=_C.energy ? (rcR*ptotL+rcL*ptotR-rcL*rcR*(u1R-u1L))/rcs
: (sR[i]*(u1L*rcL+ptotL)+sL[i]*(u1R*rcR-ptotR))
/(sR[i]-sL[i])-rhoLs*sM*sM;*/
u2L=uL[i2][i]/rhoL;
u3L=uL[i3][i]/rhoL;
u2R=uR[i2][i]/rhoR;
u3R=uR[i3][i]/rhoR;
if(HLLD_PRESSURE_CORR)
{
cuL=sqrt(2.*(ptotL-uL[ipr][i])/rhoL+SPQR(u1L,u2L,u3L));
cuR=sqrt(2.*(ptotR-uR[ipr][i])/rhoR+SPQR(u1R,u2R,u3R));
phi=2.*MAX(cuL,cuR)/(sR[i]-sL[i]-MAX(u1R,u1L)+MIN(u1R,u1L));
phi=MIN(1.,phi);
phi*=(2.-phi);
}
else phi=1.;
ptotS=(rcR*ptotL+rcL*ptotR-phi*rcL*rcR*(u1R-u1L))/rcs;
/*ptotS=_C.energy ? (rcR*ptotL+rcL*ptotR-rcL*rcR*(u1R-u1L))/rcs
: (sR[i]*(u1L*rcL+ptotL)+sL[i]*(u1R*rcR-ptotR))
/(sR[i]-sL[i])-rhoLs*sM*sM;*/
if(sM>0.)
{
e_nom=rhoL*SQR(sL[i]-u1L)-pi*SQR(b1);
e_den=rhoL*(sL[i]-sM)*(sL[i]-u1L)-pi*SQR(b1);
if(fabs(e_den)>DEPS*rhoL*smax*smax)
//if(fabs(e_den)>DEPS*rhoL*smax*smax)
if(e_den!=0.)
{
u2Ls=u2L-pi*b1*b2L*(sM-u1L)/e_den;
b2Ls=b2L*e_nom/e_den;
......@@ -215,7 +228,7 @@ for(i=is; i<=ie; i++)
uI[0][i]=rhoLs;
if(_C.energy) uI[4][i]=((sL[i]-u1L)*uL[4][i]-ptotL*u1L+ptotS*sM+b1*pi
*(u1L*b1+u2L*b2L+u3L*b3L-vdotbLs))/(sL[i]-sM);
if(sLs>0.)
if(sLs>=0.)
{
uI[i2][i]=rhoLs*u2Ls;
uI[i3][i]=rhoLs*u3Ls;
......@@ -231,7 +244,8 @@ for(i=is; i<=ie; i++)
e_nom=rhoR*SQR(sR[i]-u1R)-pi*SQR(b1);
e_den=rhoR*(sR[i]-sM)*(sR[i]-u1R)-pi*SQR(b1);
if(fabs(e_den)>DEPS*rhoR*smax*smax)
//if(fabs(e_den)>DEPS*rhoR*smax*smax)
if(e_den!=0.)
{
u2Rs=u2R-pi*b1*b2R*(sM-u1R)/e_den;
b2Rs=b2R*e_nom/e_den;
......@@ -271,7 +285,8 @@ for(i=is; i<=ie; i++)
e_nom=rhoR*SQR(sR[i]-u1R)-pi*SQR(b1);
e_den=rhoR*(sR[i]-sM)*(sR[i]-u1R)-pi*SQR(b1);
if(fabs(e_den)>DEPS*rhoR*smax*smax)
// if(fabs(e_den)>DEPS*rhoR*smax*smax)
if(e_den!=0.)
{
u2Rs=u2R-pi*b1*b2R*(sM-u1R)/e_den;
b2Rs=b2R*e_nom/e_den;
......@@ -297,7 +312,8 @@ for(i=is; i<=ie; i++)
e_nom=rhoL*SQR(sL[i]-u1L)-pi*SQR(b1);
e_den=rhoL*(sL[i]-sM)*(sL[i]-u1L)-pi*SQR(b1);
if(fabs(e_den)>DEPS*rhoL*smax*smax)
//if(fabs(e_den)>DEPS*rhoL*smax*smax)
if(e_den!=0.)
{
u2Ls=u2L-pi*b1*b2L*(sM-u1L)/e_den;
b2Ls=b2L*e_nom/e_den;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment