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

Jump to bottom

Improve the accuracy and stability of the GREP routines #123

Open
wants to merge 7 commits into
base: ccsn
Choose a base branch
from
Open

Improve the accuracy and stability of the GREP routines #123

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
4642889
Revise the GREP routine
hfhsieh May 14, 2024
6adf12c
Subtract energy shift from the energy profile for computing GREP
hfhsieh May 14, 2024
1802746
Minor
hfhsieh May 14, 2024
76eb2b8
Minor
hfhsieh May 14, 2024
548084f
Merge branch 'ccsn' into ccsn_grep_rev
hfhsieh Oct 23, 2024
c32f44f
Minor
hfhsieh Nov 5, 2024
ae645a3
[Bugfix] Fix compilation error if GRAVITY is disabled
hfhsieh Nov 19, 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
164 changes: 104 additions & 60 deletions src/SelfGravity/CPU_Poisson/CPU_ComputeGREP.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,9 @@
#include "GAMER.h"


#define LinearInterp( x, xa, xb, ya, yb ) ( (ya) + ((yb) - (ya)) * ((x) - (xa)) / ((xb) - (xa)) )




//-------------------------------------------------------------------------------------------------------
Expand All @@ -13,7 +16,7 @@
// 4. Ref: Marek et al., 2006, A&A, 445, 273 (arXiv: 0502161), sec. 2, case A
//
// Parameter : lv : Current AMR level
// Time : Current physical time at the refinement level=lv
// Time : Current physical time at the specified level
// DensAve : Profile_t object storing the spherically averaged profile of density
// EngyAve : Profile_t object storing the spherically averaged profile of internal energy density
// VrAve : Profile_t object storing the spherically averaged profile of radial velocity
Expand All @@ -30,7 +33,7 @@ void CPU_ComputeGREP( const int lv, const double Time, const Profile_t *DensAve,
const double _c2 = 1.0 / c2;
const double FourPI = 4.0 * M_PI;
const double FourThirdPI = FourPI / 3.0;
const double Tolerance = 1.0e-5;
const double Tolerance = 1.0e-10;

int NIter = 0;
int NBin = DensAve->NBin;
Expand All @@ -41,59 +44,77 @@ void CPU_ComputeGREP( const int lv, const double Time, const Profile_t *DensAve,
double *Vr = VrAve ->Data;
double *Pres = PresAve->Data;

double *Weight = new double [NBin]; // volume of each bin
double *Edge_L = new double [NBin]; // left edge of each bin
double *Vr_L = new double [NBin]; // radial velocity at the left edge of each bin
double *Mass_NW = new double [NBin]; // enclosed Newtonian mass for \bar_Phi(r) in Eq. (7)
double *Mass_TOV = new double [NBin]; // enclosed TOV mass for \bar_Phi(r)_TOV in Eq. (7)
double *Mass_TOV_USG = new double [NBin]; // enclosed TOV mass in the previous iteration
double *Mass_TOV = new double [NBin]; // enclosed TOV mass for \bar_Phi(r)_TOV in Eq. (7)
double *Mass_TOV_USG = new double [NBin]; // enclosed TOV mass in the previous iteration
double *Dens_TOV = new double [NBin]; // empirical TOV density Eq. (4)
double *Gamma_TOV = new double [NBin]; // metric function Eq. (5)
double *dVol = new double [NBin]; // volume of each bin
double *EdgeL = new double [NBin]; // left edge of each bin

for (int i=0; i<NBin; i++) Mass_TOV_USG[i] = 0.0;
for (int i=0; i<NBin; i++) Gamma_TOV [i] = 1.0;
// (1) preparation stage
for (int i=0; i<NBin; i++) Mass_TOV_USG[i] = 0.0;
for (int i=0; i<NBin; i++) Gamma_TOV [i] = 1.0;

Edge_L[0] = 0.0;
for (int i=1; i<NBin; i++) Edge_L[i] = ( GREP_LOGBIN ) ? sqrt( Radius[i] * Radius[i-1] )
: 0.5*( Radius[i] + Radius[i-1] );

// find the maximum radius that the stored weights are reliable
for (int d=0; d<3; d++)
{
MaxRadius = MIN( MIN( amr->BoxSize[d] - DensAve->Center[d], DensAve->Center[d] ),
MaxRadius );
}

EdgeL [0] = 0.0;
for (int i=1; i<NBin; i++) EdgeL [i] = ( GREP_LOGBIN ) ? sqrt( Radius[i-1] * Radius[i] )
: 0.5*( Radius[i-1] + Radius[i] );
// update the weight
for (int i=0; i<NBin-1; i++)
{
Weight[i] = ( Edge_L[i+1] > MaxRadius )
? FourThirdPI * ( CUBE(Edge_L[i+1]) - CUBE(Edge_L[i]) )
: DensAve->Weight[i];
}

for (int i=0; i<NBin-1; i++) dVol [i] = FourThirdPI * ( CUBE( EdgeL[i+1] ) - CUBE( EdgeL[i] ) );
dVol [NBin-1] = FourThirdPI * ( CUBE( MaxRadius ) - CUBE( EdgeL[NBin-1] ) );
Weight[NBin-1] = FourThirdPI * ( CUBE(DensAve->MaxRadius) - CUBE(Edge_L[NBin-1]) );

// compute the mass of each bin
for (int i=0; i<NBin; i++) Mass_NW[i] = Weight[i] * Dens[i];

// 1. compute Mass_NW at the outer edge of each bin
Mass_NW [0] = dVol [0] * Dens[0];
for (int i=1; i<NBin-1; i++) Mass_NW [i] = dVol [i] * ( 0.5 * Dens[i] + 0.25 * ( Dens[i-1] + Dens[i+1] ) );
Mass_NW[NBin-1] = dVol[NBin-1] * Dens[NBin-1];
// interpolate the radial velocity at the bin edge
Vr_L[0] = 0.0;
for (int i=1; i<NBin; i++) Vr_L[i] = LinearInterp( Edge_L[i], Radius[i-1], Radius[i],
Mass_NW[i-1] * Vr[i-1], Mass_NW[i] * Vr[i] )
/ LinearInterp( Edge_L[i], Radius[i-1], Radius[i],
Mass_NW[i-1], Mass_NW[i] );
ChunYen-Chen marked this conversation as resolved.
Show resolved Hide resolved

for (int i=1; i<NBin; i++) Mass_NW[i] += Mass_NW[i-1];
// compute the enclosed Newtonian mass, defined at the right edge of each bin
for (int i=1; i<NBin; i++) Mass_NW[i] += Mass_NW[i-1];


// 2. iteratively compute Mass_TOV and Gamma_TOV at the outer edge of each bin
// (2) iteratively compute Mass_TOV and Gamma_TOV, defined at the right edge of each bin
// --> ignore the last bin since it is not used later
bool IsConverged = false;

while ( ! IsConverged && ( NIter++ < GREP_MAXITER ) )
{
// update Mass_TOV
for (int i=0; i<NBin; i++) Dens_TOV[i] = Gamma_TOV[i] * ( Dens[i] + Engy[i] * _c2 );
// --> include the last bin here for troubleshooting information
for (int i=0; i<NBin; i++) Dens_TOV[i] = Gamma_TOV[i] * ( Dens[i] + Engy[i] * _c2 );

Mass_TOV [0] = dVol [0] * Dens_TOV[0];
for (int i=1; i<NBin-1; i++) Mass_TOV [i] = dVol [i] * ( 0.5 * Dens_TOV[i] + 0.25 * ( Dens_TOV[i-1] + Dens_TOV[i+1] ) );
Mass_TOV[NBin-1] = dVol[NBin-1] * Dens_TOV[NBin-1];

for (int i=1; i<NBin ; i++) Mass_TOV[i] += Mass_TOV[i-1];
for (int i=0; i<NBin; i++) Mass_TOV[i] = Weight[i] * Dens_TOV[i];
for (int i=1; i<NBin; i++) Mass_TOV[i] += Mass_TOV[i-1];

// update Gamma_TOV
for (int i=0; i<NBin-1; i++) Gamma_TOV [i] = 1.0 + ( 0.25 * SQR( Vr[i] + Vr[i+1] )
- 2.0 * NEWTON_G * Mass_TOV[i] / EdgeL[i+1] ) * _c2;
Gamma_TOV[NBin-1] = 1.0 + ( SQR( Vr[NBin-1] )
- 2.0 * NEWTON_G * Mass_TOV[NBin-1] / MaxRadius ) * _c2;
for (int i=0; i<NBin-1; i++) Gamma_TOV[i] = 1.0 + ( SQR( Vr_L[i+1] )
- 2.0 * NEWTON_G * Mass_TOV[i] / Edge_L[i+1] ) * _c2;

for (int i=0; i<NBin; i++) Gamma_TOV[i] = sqrt( MAX( TINY_NUMBER, Gamma_TOV[i] ) );
for (int i=0; i<NBin-1; i++) Gamma_TOV[i] = sqrt( MAX( TINY_NUMBER, Gamma_TOV[i] ) );

// check whether the Mass_TOV is converged
IsConverged = true;

for (int i=0; i<NBin; i++)
for (int i=0; i<NBin-1; i++)
{
double RelErr = fabs( Mass_TOV_USG[i] - Mass_TOV[i] ) / Mass_TOV[i];

Expand Down Expand Up @@ -138,7 +159,7 @@ void CPU_ComputeGREP( const int lv, const double Time, const Profile_t *DensAve,

// data
for (int i=0; i<NBin; i++)
fprintf( File, "%5d %9ld %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e\n",
fprintf( File, "%5d %9ld %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e\n",
i, DensAve->NCell[i], Radius[i], Dens[i], Engy[i], Vr[i], Pres[i],
Mass_NW[i], Mass_TOV[i], Mass_TOV_USG[i], Gamma_TOV[i],
fabs( Mass_TOV_USG[i] - Mass_TOV[i] ) / Mass_TOV[i] );
Expand All @@ -151,37 +172,40 @@ void CPU_ComputeGREP( const int lv, const double Time, const Profile_t *DensAve,
} // if ( ! IsConverged )


// 3. compute the effective GR potential at the bin center
// (3) compute the effective GR potential, defined at the left edge of each bin
EffPot->NBin = DensAve->NBin;
EffPot->LogBin = DensAve->LogBin;
EffPot->LogBinRatio = DensAve->LogBinRatio;
EffPot->MaxRadius = DensAve->MaxRadius;
EffPot->AllocateMemory();

for (int d=0; d<3; d++) EffPot->Center[d] = DensAve->Center[d];
for (int b=0; b<NBin; b++) EffPot->Radius[b] = DensAve->Radius[b];
for (int b=0; b<NBin; b++) EffPot->Radius[b] = Edge_L[b];

// set the outer boundary condition of the potential to -G * M_out / r_out
EffPot->Data[NBin-1] = -NEWTON_G * ( Mass_TOV[NBin-1] - Mass_NW[NBin-1] ) / Radius[NBin-1];
// set the outer boundary condition of the potential to -GM/r at the left edge of last bin
EffPot->Data[NBin-1] = -NEWTON_G * ( Mass_TOV[NBin-2] - Mass_NW[NBin-2] ) / EffPot->Radius[NBin-1];

for (int i=NBin-2; i>0; i--)
{
double dr = Radius[i] - Radius[i+1];
double Mass_NW_C = 0.5 * ( Mass_NW[i] + Mass_NW[i-1] );
double Mass_TOV_C = ( 0.5 * ( Mass_TOV[i] + Mass_TOV[i-1] ) + FourPI * CUBE( Radius[i] ) * Pres[i] * _c2 )
* ( 1.0 + ( Engy[i] + Pres[i] ) / ( Dens[i] * c2 ) )
/ SQR( 0.5 * ( Gamma_TOV[i] + Gamma_TOV[i-1] ) );
const double dr = EffPot->Radius[i] - EffPot->Radius[i+1];

// compute the integrand using data at the bin center
// --> approximate the Mass_NW, Mass_TOV, and Gamma_TOV at the bin center via linear interpolation
const double ratio = ( Radius[i] - Edge_L[i] ) / ( Edge_L[i+1] - Edge_L[i] );

double Mass_NW_C = Mass_NW [i-1] + ratio * ( Mass_NW [i] - Mass_NW [i-1] );
double Mass_TOV_C = Mass_TOV [i-1] + ratio * ( Mass_TOV [i] - Mass_TOV [i-1] );
double Gamma_TOV_C = Gamma_TOV[i-1] + ratio * ( Gamma_TOV[i] - Gamma_TOV[i-1] );

Mass_TOV_C = ( Mass_TOV_C + FourPI * CUBE( Radius[i] ) * Pres[i] * _c2 )
* ( 1.0 + ( Engy[i] + Pres[i] ) / ( Dens[i] * c2 ) ) / SQR( Gamma_TOV_C );


EffPot->Data[i] = EffPot->Data[i+1] - dr * NEWTON_G * ( Mass_NW_C - Mass_TOV_C ) / SQR( Radius[i] );
}

double dr = Radius[0] - Radius[1];
double Mass_NW_C = Mass_NW[0];
double Mass_TOV_C = ( Mass_TOV[0] + FourPI * CUBE( Radius[0] ) * Pres[0] * _c2 )
* ( 1.0 + ( Engy[0] + Pres[0] ) / ( Dens[0] * c2 ) )
/ SQR( Gamma_TOV[0] );

EffPot->Data[0] = EffPot->Data[1] - dr * NEWTON_G * ( Mass_NW_C - Mass_TOV_C ) / SQR( Radius[0] ) ;
// set the potential for the innermost bin
EffPot->Data[0] = EffPot->Data[1];


// troubleshooting information
Expand All @@ -195,37 +219,57 @@ void CPU_ComputeGREP( const int lv, const double Time, const Profile_t *DensAve,

// metadata
fprintf( File, "# Step : %ld\n", Step );
fprintf( File, "# Time : %.7e\n", Time );
fprintf( File, "# Time : %13.7e\n", Time );
fprintf( File, "# GREP_CENTER_METHOD : %d\n", GREP_CENTER_METHOD );
fprintf( File, "# Center : %13.7e %13.7e %13.7e\n", EffPot->Center[0], EffPot->Center[1], EffPot->Center[2] );
fprintf( File, "# Maximum Radius : %13.7e\n", EffPot->MaxRadius );
fprintf( File, "# LogBin : %d\n", EffPot->LogBin );
fprintf( File, "# LogBinRatio : %13.7e\n", EffPot->LogBinRatio );
fprintf( File, "# Number of Iterations : %d\n", NIter );
fprintf( File, "# NBin : %d\n", NBin );
fprintf( File, "# -------------------------------------------------\n" );
fprintf( File, "%5s %9s %22s %22s %22s %22s %22s %22s %22s %22s %23s\n",
"# Bin", "NCell", "Radius", "Density", "Energy", "Vr", "Pressure",
"Mass_NW", "Mass_TOV", "Gamma_TOV", "Eff_Pot");
fprintf( File, "# NBin : %d\n", NBin );
fprintf( File, "#\n" );
fprintf( File, "# Bin : bin index\n" );
fprintf( File, "# NCell : number of cell\n" );
fprintf( File, "# Bin_Center : bin center\n" );
fprintf( File, "# Bin_Edge_L : left bin edge\n" );
fprintf( File, "# Density : density defined at the center\n" );
fprintf( File, "# Energy : internal energy density defined at the center\n" );
fprintf( File, "# Vr : radial velocity defined at the center\n" );
fprintf( File, "# Pressure : pressure defined at the center\n" );
fprintf( File, "# Mass_NW : enclosed Newtonian mass defined at the right edge\n" );
fprintf( File, "# Mass_TOV : enclosed TOV mass defined at the right edge\n" );
fprintf( File, "# Gamma_TOV : metric function defined at the right edge\n" );
fprintf( File, "# Eff_Pot : effective potential correction defined at the left edge\n" );
fprintf( File, "# -------------------------------------------------------------------------------\n" );
fprintf( File, "%5s %9s %22s %22s %22s %22s %22s %22s %22s %22s %22s %22s\n",
"# [1]", "[2]", "[3]", "[4]", "[5]", "[6]", "[7]", "[8]", "[9]", "[10]", "[11]", "[12]" );
fprintf( File, "%5s %9s %22s %22s %22s %22s %22s %22s %22s %22s %22s %22s\n",
"# Bin", "NCell", "Bin_Center", "Bin_Edge_L", "Density", "Energy", "Vr", "Pressure",
"Mass_NW", "Mass_TOV", "Gamma_TOV", "Eff_Pot" );

// data
for (int i=0; i<NBin; i++)
fprintf( File, "%5d %9ld %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %23.15e\n",
i, DensAve->NCell[i], Radius[i], Dens[i], Engy[i], Vr[i], Pres[i],
Mass_NW[i], Mass_TOV[i], Gamma_TOV[i], EffPot->Data[i] );
fprintf( File, "%5d %9ld %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e\n",
1, DensAve->NCell[0], Radius[0], Edge_L[0], Dens[0], Engy[0], Vr[0], Pres[0],
NULL_REAL, NULL_REAL, NULL_REAL, EffPot->Data[0] );

for (int i=1; i<NBin; i++)
fprintf( File, "%5d %9ld %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e %22.14e\n",
i+1, DensAve->NCell[i], Radius[i], Edge_L[i], Dens[i], Engy[i], Vr[i], Pres[i],
Mass_NW[i-1], Mass_TOV[i-1], Gamma_TOV[i-1], EffPot->Data[i] );

fclose( File );
}
# endif


// free memory
delete [] Weight;
delete [] Edge_L;
delete [] Vr_L;
delete [] Mass_NW;
delete [] Mass_TOV;
delete [] Mass_TOV_USG;
delete [] Dens_TOV;
delete [] Gamma_TOV;
delete [] dVol;
delete [] EdgeL;

} // FUNCTION : CPU_ComputeGREP
2 changes: 1 addition & 1 deletion src/SelfGravity/CPU_Poisson/CPU_ExtPot_GREP.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#endif


#define LinearInterp( x, xa, xb, ya, yb ) ( ( ((x) - (xa)) * (yb) + ((xb) - (x)) * (ya) ) / ((xb) - (xa)) )
#define LinearInterp( x, xa, xb, ya, yb ) ( (ya) + ((yb) - (ya)) * ((x) - (xa)) / ((xb) - (xa)) )


#ifdef __CUDACC__
Expand Down
10 changes: 10 additions & 0 deletions src/SelfGravity/Poi_UserWorkBeforePoisson_GREP.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -308,6 +308,16 @@ void Poi_Prepare_GREP( const double Time, const int lv )
} // if ( GREP_OPT_PRES == GREP_PRES_BINDATA )


// subtract the energy shift from the profile of internal energy density
# if ( EOS == EOS_NUCLEAR )
const real sEint2CGS = EoS_AuxArray_Flt[NUC_AUX_VSQR2CGS];
const real EnergyShift_Code = EoS_AuxArray_Flt[NUC_AUX_ESHIFT ] / sEint2CGS;

for (int b=0; b<Engy_Tot->NBin; b++)
Engy_Tot->Data[b] -= EnergyShift_Code * Dens_Tot->Data[b];
# endif


// check the profiles before computing the effective GR potential
Profile_t *GREP_Check_List[] = { Dens_Tot, Engy_Tot, Vr_Tot, Pres_Tot };

Expand Down
27 changes: 21 additions & 6 deletions src/TestProblem/Hydro/CCSN/Record_CCSN.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -387,6 +387,14 @@ void Detect_CoreBounce()


// (2) criterion 2: any cells within 30km has entropy larger than 3
double Center[3];

# ifdef GRAVITY
for (int i=0; i<3; i++) Center[i] = GREP_Center[i]
# else
for (int i=0; i<3; i++) Center[i] = amr->BoxCenter[i];
# endif


// allocate memory for per-thread arrays
# ifdef OPENMP
Expand Down Expand Up @@ -423,9 +431,9 @@ void Detect_CoreBounce()
for (int j=0; j<PS1; j++) { const double y = amr->patch[0][lv][PID]->EdgeL[1] + (j+0.5)*dh;
for (int i=0; i<PS1; i++) { const double x = amr->patch[0][lv][PID]->EdgeL[0] + (i+0.5)*dh;

const double x0 = x - GREP_Center[0];
const double y0 = y - GREP_Center[1];
const double z0 = z - GREP_Center[2];
const double x0 = x - Center[0];
const double y0 = y - Center[1];
const double z0 = z - Center[2];
const double r = sqrt( SQR( x0 ) + SQR( y0 ) + SQR( z0 ) );

// ignore cells outside 30km
Expand Down Expand Up @@ -498,6 +506,7 @@ void Detect_Shock()
const int PRES = NCOMP_PASSIVE + 4;


double Center[3];
double Shock_Min = HUGE_NUMBER;
double Shock_Max = -HUGE_NUMBER;
double Shock_Ave = 0.0;
Expand All @@ -523,6 +532,12 @@ void Detect_Shock()
OMP_Shock_Found [t] = false;
}

# ifdef GRAVITY
for (int i=0; i<3; i++) Center[i] = GREP_Center[i]
# else
for (int i=0; i<3; i++) Center[i] = amr->BoxCenter[i];
# endif


for (int lv=0; lv<NLEVEL; lv++)
{
Expand Down Expand Up @@ -614,9 +629,9 @@ void Detect_Shock()
for (int j=0; j<PS1; j++) { const double y = amr->patch[0][lv][PID]->EdgeL[1] + (j+0.5)*dh; const int jj = j+NGhost;
for (int i=0; i<PS1; i++) { const double x = amr->patch[0][lv][PID]->EdgeL[0] + (i+0.5)*dh; const int ii = i+NGhost;

const double dx = x - GREP_Center[0];
const double dy = y - GREP_Center[1];
const double dz = z - GREP_Center[2];
const double dx = x - Center[0];
const double dy = y - Center[1];
const double dz = z - Center[2];
const double r = sqrt( SQR( dx ) + SQR( dy ) + SQR( dz ) );

// (3) evaluate the undivided gradient of pressure and the undivided divergence of velocity
Expand Down