Reuse mass integrand across detector radii

This commit is contained in:
2026-04-13 11:55:41 +08:00
parent 3a58273501
commit 4b10519876
3 changed files with 226 additions and 206 deletions

View File

@@ -7114,13 +7114,17 @@ void bssn_class::AnalysisStuff(int lev, double dT_lev)
for (int pl = 2; pl < maxl + 1; pl++)
for (int pm = -pl; pm < pl + 1; pm++)
NN++;
RP = new double[NN];
IP = new double[NN];
RoutMAP = new double[7];
double Rex = maxrex;
for (int i = 0; i < decn; i++)
{
#ifdef Point_Psi4
RP = new double[NN];
IP = new double[NN];
RoutMAP = new double[7];
double Rex = maxrex;
bool patch_mass_prepared = false;
#ifdef WithShell
bool shell_mass_prepared = false;
#endif
for (int i = 0; i < decn; i++)
{
#ifdef Point_Psi4
Waveshell->surf_Wave(Rex, GH, SH,
phi, trK,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
@@ -7139,69 +7143,76 @@ void bssn_class::AnalysisStuff(int lev, double dT_lev)
Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
2, maxl, NN, RP, IP, ErrorMonitor);
#ifdef WithShell
if (lev > 0 || Rex < GH->bbox[0][0][3])
{
Waveshell->surf_MassPAng(Rex, lev, GH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor);
}
else
{
Waveshell->surf_MassPAng(Rex, lev, SH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor);
}
#else
Waveshell->surf_MassPAng(Rex, lev, GH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor);
#endif
#else
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before surface integral");
#ifdef WithShell
if (lev > 0 || Rex < GH->bbox[0][0][3])
if (lev > 0 || Rex < GH->bbox[0][0][3])
{
Waveshell->surf_MassPAng(Rex, lev, GH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor, !patch_mass_prepared);
patch_mass_prepared = true;
}
else
{
Waveshell->surf_MassPAng(Rex, lev, SH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor, !shell_mass_prepared);
shell_mass_prepared = true;
}
#else
Waveshell->surf_MassPAng(Rex, lev, GH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor, !patch_mass_prepared);
patch_mass_prepared = true;
#endif
#else
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before surface integral");
#ifdef WithShell
if (lev > 0 || Rex < GH->bbox[0][0][3])
{
Waveshell->surf_Wave(Rex, lev, GH, Rpsi4, Ipsi4, 2, maxl, NN, RP, IP, ErrorMonitor);
Waveshell->surf_MassPAng(Rex, lev, GH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor);
}
else
{
Waveshell->surf_Wave(Rex, lev, SH, Rpsi4, Ipsi4, 2, maxl, NN, RP, IP, ErrorMonitor);
Waveshell->surf_MassPAng(Rex, lev, SH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor);
}
#else
#if (PSTR == 0)
Waveshell->surf_Wave(Rex, lev, GH, Rpsi4, Ipsi4, 2, maxl, NN, RP, IP, ErrorMonitor);
Waveshell->surf_MassPAng(Rex, lev, GH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor);
#elif (PSTR == 1 || PSTR == 2)
Waveshell->surf_Wave(Rex, lev, GH, Rpsi4, Ipsi4, 2, maxl, NN, RP, IP, ErrorMonitor, GH->Commlev[lev]);
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"after surf_Wave");
Waveshell->surf_MassPAng(Rex, lev, GH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor, GH->Commlev[lev]);
#endif
#endif
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"end surface integral");
Waveshell->surf_MassPAng(Rex, lev, GH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor, !patch_mass_prepared);
patch_mass_prepared = true;
}
else
{
Waveshell->surf_Wave(Rex, lev, SH, Rpsi4, Ipsi4, 2, maxl, NN, RP, IP, ErrorMonitor);
Waveshell->surf_MassPAng(Rex, lev, SH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor, !shell_mass_prepared);
shell_mass_prepared = true;
}
#else
#if (PSTR == 0)
Waveshell->surf_Wave(Rex, lev, GH, Rpsi4, Ipsi4, 2, maxl, NN, RP, IP, ErrorMonitor);
Waveshell->surf_MassPAng(Rex, lev, GH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor, !patch_mass_prepared);
patch_mass_prepared = true;
#elif (PSTR == 1 || PSTR == 2)
Waveshell->surf_Wave(Rex, lev, GH, Rpsi4, Ipsi4, 2, maxl, NN, RP, IP, ErrorMonitor, GH->Commlev[lev]);
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"after surf_Wave");
Waveshell->surf_MassPAng(Rex, lev, GH, phi0, trK0,
gxx0, gxy0, gxz0, gyy0, gyz0, gzz0,
Axx0, Axy0, Axz0, Ayy0, Ayz0, Azz0,
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1, // here we can not touch rhs variables, but 1 variables
RoutMAP, ErrorMonitor, GH->Commlev[lev], !patch_mass_prepared);
patch_mass_prepared = true;
#endif
#endif
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"end surface integral");
#endif
if (i == 0)
{