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 5 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
138 changes: 85 additions & 53 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 *Dens_TOV = new double [NBin]; // empirical TOV density Eq. (4)
double *Gamma_TOV = new double [NBin]; // metric function Eq. (5)
ChunYen-Chen marked this conversation as resolved.
Show resolved Hide resolved
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;

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] );
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 );
}

// 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 );

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];
// --> 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 );

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 @@ -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 @@ -204,28 +228,36 @@ void CPU_ComputeGREP( const int lv, const double Time, const Profile_t *DensAve,
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",
fprintf( File, "%5s %9s %22s %22s %22s %22s %22s %22s %22s %22s %22s %23s\n",
"# Bin", "NCell", "Bin_Center", "Bin_Edge", "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] );
// --> Dens, Engy, Vr, Pres are defined at the bin center
// Mass_NW, Mass_TOV, Gamma_TOV are defined at the right edge
// EffPot is defined at the left edge
ChunYen-Chen marked this conversation as resolved.
Show resolved Hide resolved
fprintf( File, "%5d %9ld %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %23.15e\n",
0, 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.15e %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], 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