64-BitBrainstorm_2026/AMSS-NCKU
8
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge wave and mass extraction interpolation

Browse Source 
(cherry picked from commit f3988ac8ca)
This commit is contained in:
CGH0S7
2026-04-13 13:17:36 +08:00
committed by ianchb
parent cce8a44fc4
commit 11977eb82f
3 changed files with 831 additions and 192 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
AMSS_NCKU_source/bssn_class.C
View File

@@ -7678,31 +7678,34 @@ void bssn_class::AnalysisStuff(int lev, double dT_lev)
#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,
Waveshell->surf_WaveMassPAng(Rex, lev, GH,
Rpsi4, Ipsi4, 2, maxl, NN, RP, IP,
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
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1,
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,
Waveshell->surf_WaveMassPAng(Rex, lev, SH,
Rpsi4, Ipsi4, 2, maxl, NN, RP, IP,
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
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1,
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,
Waveshell->surf_WaveMassPAng(Rex, lev, GH,
Rpsi4, Ipsi4, 2, maxl, NN, RP, IP,
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
Gmx0, Gmy0, Gmz0, Sfx1, Sfy1, Sfz1,
RoutMAP, ErrorMonitor, !patch_mass_prepared);
patch_mass_prepared = true;
#elif (PSTR == 1 || PSTR == 2)

620
AMSS_NCKU_source/surface_integral.C
View File

@@ -3272,6 +3272,626 @@ void surface_integral::surf_MassPAng(double rex, int lev, ShellPatch *GH, var *c
delete[] shellf;
DG_List->clearList();
}
void surface_integral::surf_WaveMassPAng(double rex, int lev, cgh *GH,
var *Rpsi4, var *Ipsi4, int spinw, int maxl, int NN, double *RP, double *IP,
var *chi, var *trK,
var *gxx, var *gxy, var *gxz, var *gyy, var *gyz, var *gzz,
var *Axx, var *Axy, var *Axz, var *Ayy, var *Ayz, var *Azz,
var *Gmx, var *Gmy, var *Gmz,
var *Sfx_rhs, var *Sfy_rhs, var *Sfz_rhs,
double *Rout, monitor *Monitor, bool refresh_mass_fields)
{
if (Symmetry != 0 && Symmetry != 1)
{
surf_Wave(rex, lev, GH, Rpsi4, Ipsi4, spinw, maxl, NN, RP, IP, Monitor);
surf_MassPAng(rex, lev, GH, chi, trK,
gxx, gxy, gxz, gyy, gyz, gzz,
Axx, Axy, Axz, Ayy, Ayz, Azz,
Gmx, Gmy, Gmz,
Sfx_rhs, Sfy_rhs, Sfz_rhs,
Rout, Monitor, refresh_mass_fields);
return;
}
if (myrank == 0 && GH->grids[lev] != 1)
if (Monitor && Monitor->outfile)
Monitor->outfile << "WARNING: surface integral on multipatches" << endl;
else
cout << "WARNING: surface integral on multipatches" << endl;
if (refresh_mass_fields)
{
MyList<Patch> *Pp = GH->PatL[lev];
while (Pp)
{
MyList<Block> *BP = Pp->data->blb;
while (BP)
{
Block *cg = BP->data;
if (myrank == cg->rank)
{
f_admmass_bssn(cg->shape, cg->X[0], cg->X[1], cg->X[2],
cg->fgfs[chi->sgfn], cg->fgfs[trK->sgfn],
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn],
cg->fgfs[Gmx->sgfn], cg->fgfs[Gmy->sgfn], cg->fgfs[Gmz->sgfn],
cg->fgfs[Sfx_rhs->sgfn], cg->fgfs[Sfy_rhs->sgfn], cg->fgfs[Sfz_rhs->sgfn],
Symmetry);
}
if (BP == Pp->data->ble)
break;
BP = BP->next;
}
Pp = Pp->next;
}
}
const int InList = 19;
const int idx_rpsi4 = 0, idx_ipsi4 = 1;
const int idx_sfx = 2, idx_sfy = 3, idx_sfz = 4;
const int idx_chi = 5, idx_trk = 6;
const int idx_gxx = 7, idx_gxy = 8, idx_gxz = 9, idx_gyy = 10, idx_gyz = 11, idx_gzz = 12;
const int idx_axx = 13, idx_axy = 14, idx_axz = 15, idx_ayy = 16, idx_ayz = 17, idx_azz = 18;
MyList<var> *DG_List = new MyList<var>(Rpsi4);
DG_List->insert(Ipsi4);
DG_List->insert(Sfx_rhs);
DG_List->insert(Sfy_rhs);
DG_List->insert(Sfz_rhs);
DG_List->insert(chi);
DG_List->insert(trK);
DG_List->insert(gxx);
DG_List->insert(gxy);
DG_List->insert(gxz);
DG_List->insert(gyy);
DG_List->insert(gyz);
DG_List->insert(gzz);
DG_List->insert(Axx);
DG_List->insert(Axy);
DG_List->insert(Axz);
DG_List->insert(Ayy);
DG_List->insert(Ayz);
DG_List->insert(Azz);
int n;
double *pox[3];
for (int ia = 0; ia < 3; ia++)
pox[ia] = new double[n_tot];
for (n = 0; n < n_tot; n++)
{
pox[0][n] = rex * nx_g[n];
pox[1][n] = rex * ny_g[n];
pox[2][n] = rex * nz_g[n];
}
int mp, Lp, Nmin, Nmax;
mp = n_tot / cpusize;
Lp = n_tot - cpusize * mp;
if (Lp > myrank)
{
Nmin = myrank * mp + myrank;
Nmax = Nmin + mp;
}
else
{
Nmin = myrank * mp + Lp;
Nmax = Nmin + mp - 1;
}
double *shellf = new double[n_tot * InList];
GH->PatL[lev]->data->Interp_Points(DG_List, n_tot, pox, shellf, Symmetry, Nmin, Nmax);
double *RP_out = new double[NN];
double *IP_out = new double[NN];
for (int ii = 0; ii < NN; ii++)
{
RP_out[ii] = 0.0;
IP_out[ii] = 0.0;
}
double Mass_out = 0.0;
double ang_outx = 0.0, ang_outy = 0.0, ang_outz = 0.0;
double p_outx = 0.0, p_outy = 0.0, p_outz = 0.0;
const double f1o8 = 0.125;
build_wave_cache(spinw, maxl);
for (n = Nmin; n <= Nmax; n++)
{
const int base = InList * n;
const int i = int(n / N_phi);
const int j = n - i * N_phi;
const double psi4RR0 = shellf[base + idx_rpsi4];
const double psi4II0 = shellf[base + idx_ipsi4];
const double psi4RR1 = Rpsi4->SoA[2] * psi4RR0;
const double psi4II1 = Ipsi4->SoA[2] * psi4II0;
for (int countlm = 0; countlm < wave_cache_modes; countlm++)
{
const int theta_idx = countlm * N_theta + i;
const int phi_idx = countlm * N_phi + j;
const double theta_pos = wave_theta_pos[theta_idx];
const double cosmphi_here = wave_phi_cos[phi_idx];
const double sinmphi_here = wave_phi_sin[phi_idx];
RP_out[countlm] += theta_pos * (psi4RR0 * cosmphi_here + psi4II0 * sinmphi_here);
IP_out[countlm] += theta_pos * (psi4II0 * cosmphi_here - psi4RR0 * sinmphi_here);
if (Symmetry == 1)
{
const double theta_neg = wave_theta_neg[theta_idx];
RP_out[countlm] += theta_neg * (psi4RR1 * cosmphi_here + psi4II1 * sinmphi_here);
IP_out[countlm] += theta_neg * (psi4II1 * cosmphi_here - psi4RR1 * sinmphi_here);
}
}
double Chi = shellf[base + idx_chi];
double TRK = shellf[base + idx_trk];
double Gxx = shellf[base + idx_gxx] + 1.0;
double Gxy = shellf[base + idx_gxy];
double Gxz = shellf[base + idx_gxz];
double Gyy = shellf[base + idx_gyy] + 1.0;
double Gyz = shellf[base + idx_gyz];
double Gzz = shellf[base + idx_gzz] + 1.0;
double axx = shellf[base + idx_axx];
double axy = shellf[base + idx_axy];
double axz = shellf[base + idx_axz];
double ayy = shellf[base + idx_ayy];
double ayz = shellf[base + idx_ayz];
double azz = shellf[base + idx_azz];
Chi = 1.0 / (1.0 + Chi);
const double Psi = Chi * sqrt(Chi);
#ifdef GaussInt
const double theta_weight = wtcostheta[i];
Mass_out += (shellf[base + idx_sfx] * nx_g[n] + shellf[base + idx_sfy] * ny_g[n] + shellf[base + idx_sfz] * nz_g[n]) * theta_weight;
#else
const double theta_weight = 1.0;
Mass_out += shellf[base + idx_sfx] * nx_g[n] + shellf[base + idx_sfy] * ny_g[n] + shellf[base + idx_sfz] * nz_g[n];
#endif
double detg = Gxx * Gyy * Gzz + Gxy * Gyz * Gxz + Gxz * Gxy * Gyz -
Gxz * Gyy * Gxz - Gxy * Gxy * Gzz - Gxx * Gyz * Gyz;
const double gupxx = (Gyy * Gzz - Gyz * Gyz) / detg;
const double gupxy = -(Gxy * Gzz - Gyz * Gxz) / detg;
const double gupxz = (Gxy * Gyz - Gyy * Gxz) / detg;
const double gupyy = (Gxx * Gzz - Gxz * Gxz) / detg;
const double gupyz = -(Gxx * Gyz - Gxy * Gxz) / detg;
const double gupzz = (Gxx * Gyy - Gxy * Gxy) / detg;
const double aupxx = gupxx * axx + gupxy * axy + gupxz * axz;
const double aupxy = gupxx * axy + gupxy * ayy + gupxz * ayz;
const double aupxz = gupxx * axz + gupxy * ayz + gupxz * azz;
const double aupyx = gupxy * axx + gupyy * axy + gupyz * axz;
const double aupyy = gupxy * axy + gupyy * ayy + gupyz * ayz;
const double aupyz = gupxy * axz + gupyy * ayz + gupyz * azz;
const double aupzx = gupxz * axx + gupyz * axy + gupzz * axz;
const double aupzy = gupxz * axy + gupyz * ayy + gupzz * ayz;
const double aupzz = gupxz * axz + gupyz * ayz + gupzz * azz;
if (Symmetry == 0)
{
ang_outx += f1o8 * Psi * (nx_g[n] * (pox[1][n] * aupxz - pox[2][n] * aupxy) + ny_g[n] * (pox[1][n] * aupyz - pox[2][n] * aupyy) + nz_g[n] * (pox[1][n] * aupzz - pox[2][n] * aupzy)) * theta_weight;
ang_outy += f1o8 * Psi * (nx_g[n] * (pox[2][n] * aupxx - pox[0][n] * aupxz) + ny_g[n] * (pox[2][n] * aupyx - pox[0][n] * aupyz) + nz_g[n] * (pox[2][n] * aupzx - pox[0][n] * aupzz)) * theta_weight;
ang_outz += f1o8 * Psi * (nx_g[n] * (pox[0][n] * aupxy - pox[1][n] * aupxx) + ny_g[n] * (pox[0][n] * aupyy - pox[1][n] * aupyx) + nz_g[n] * (pox[0][n] * aupzy - pox[1][n] * aupzx)) * theta_weight;
}
else
{
ang_outz += f1o8 * Psi * (nx_g[n] * (pox[0][n] * aupxy - pox[1][n] * aupxx) + ny_g[n] * (pox[0][n] * aupyy - pox[1][n] * aupyx) + nz_g[n] * (pox[0][n] * aupzy - pox[1][n] * aupzx)) * theta_weight;
}
axx = Chi * (axx + Gxx * TRK / 3.0);
axy = Chi * (axy + Gxy * TRK / 3.0);
axz = Chi * (axz + Gxz * TRK / 3.0);
ayy = Chi * (ayy + Gyy * TRK / 3.0);
ayz = Chi * (ayz + Gyz * TRK / 3.0);
azz = Chi * (azz + Gzz * TRK / 3.0);
axx -= TRK;
ayy -= TRK;
azz -= TRK;
p_outx += f1o8 * Psi * (nx_g[n] * axx + ny_g[n] * axy + nz_g[n] * axz) * theta_weight;
p_outy += f1o8 * Psi * (nx_g[n] * axy + ny_g[n] * ayy + nz_g[n] * ayz) * theta_weight;
if (Symmetry == 0)
p_outz += f1o8 * Psi * (nx_g[n] * axz + ny_g[n] * ayz + nz_g[n] * azz) * theta_weight;
}
for (int ii = 0; ii < NN; ii++)
{
#ifdef GaussInt
RP_out[ii] = RP_out[ii] * rex * dphi;
IP_out[ii] = IP_out[ii] * rex * dphi;
#else
RP_out[ii] = RP_out[ii] * rex * dphi * dcostheta;
IP_out[ii] = IP_out[ii] * rex * dphi * dcostheta;
#endif
}
double mass, px, py, pz, sx, sy, sz;
{
double *reduce_out = new double[2 * NN + 7];
double *reduce_in = new double[2 * NN + 7];
memcpy(reduce_out, RP_out, NN * sizeof(double));
memcpy(reduce_out + NN, IP_out, NN * sizeof(double));
reduce_out[2 * NN + 0] = Mass_out;
reduce_out[2 * NN + 1] = ang_outx;
reduce_out[2 * NN + 2] = ang_outy;
reduce_out[2 * NN + 3] = ang_outz;
reduce_out[2 * NN + 4] = p_outx;
reduce_out[2 * NN + 5] = p_outy;
reduce_out[2 * NN + 6] = p_outz;
MPI_Allreduce(reduce_out, reduce_in, 2 * NN + 7, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
memcpy(RP, reduce_in, NN * sizeof(double));
memcpy(IP, reduce_in + NN, NN * sizeof(double));
mass = reduce_in[2 * NN + 0];
sx = reduce_in[2 * NN + 1];
sy = reduce_in[2 * NN + 2];
sz = reduce_in[2 * NN + 3];
px = reduce_in[2 * NN + 4];
py = reduce_in[2 * NN + 5];
pz = reduce_in[2 * NN + 6];
delete[] reduce_out;
delete[] reduce_in;
}
#ifdef GaussInt
mass = mass * rex * rex * dphi * factor;
sx = sx * rex * rex * dphi * (1.0 / PI) * factor;
sy = sy * rex * rex * dphi * (1.0 / PI) * factor;
sz = sz * rex * rex * dphi * (1.0 / PI) * factor;
px = px * rex * rex * dphi * (1.0 / PI) * factor;
py = py * rex * rex * dphi * (1.0 / PI) * factor;
pz = pz * rex * rex * dphi * (1.0 / PI) * factor;
#else
mass = mass * rex * rex * dphi * dcostheta * factor;
sx = sx * rex * rex * dphi * dcostheta * (1.0 / PI) * factor;
sy = sy * rex * rex * dphi * dcostheta * (1.0 / PI) * factor;
sz = sz * rex * rex * dphi * dcostheta * (1.0 / PI) * factor;
px = px * rex * rex * dphi * dcostheta * (1.0 / PI) * factor;
py = py * rex * rex * dphi * dcostheta * (1.0 / PI) * factor;
pz = pz * rex * rex * dphi * dcostheta * (1.0 / PI) * factor;
#endif
Rout[0] = mass;
Rout[1] = px;
Rout[2] = py;
Rout[3] = pz;
Rout[4] = sx;
Rout[5] = sy;
Rout[6] = sz;
delete[] pox[0];
delete[] pox[1];
delete[] pox[2];
delete[] shellf;
delete[] RP_out;
delete[] IP_out;
DG_List->clearList();
}
void surface_integral::surf_WaveMassPAng(double rex, int lev, ShellPatch *GH,
var *Rpsi4, var *Ipsi4, int spinw, int maxl, int NN, double *RP, double *IP,
var *chi, var *trK,
var *gxx, var *gxy, var *gxz, var *gyy, var *gyz, var *gzz,
var *Axx, var *Axy, var *Axz, var *Ayy, var *Ayz, var *Azz,
var *Gmx, var *Gmy, var *Gmz,
var *Sfx_rhs, var *Sfy_rhs, var *Sfz_rhs,
double *Rout, monitor *Monitor, bool refresh_mass_fields)
{
if (Symmetry != 0 && Symmetry != 1)
{
surf_Wave(rex, lev, GH, Rpsi4, Ipsi4, spinw, maxl, NN, RP, IP, Monitor);
surf_MassPAng(rex, lev, GH, chi, trK,
gxx, gxy, gxz, gyy, gyz, gzz,
Axx, Axy, Axz, Ayy, Ayz, Azz,
Gmx, Gmy, Gmz,
Sfx_rhs, Sfy_rhs, Sfz_rhs,
Rout, Monitor, refresh_mass_fields);
return;
}
if (lev != 0)
{
if (myrank == 0)
{
if (Monitor && Monitor->outfile)
Monitor->outfile << "WARNING: shell surface integral not on level 0" << endl;
else
cout << "WARNING: shell surface integral not on level 0" << endl;
}
return;
}
if (refresh_mass_fields)
{
MyList<ss_patch> *Pp = GH->PatL;
while (Pp)
{
MyList<Block> *BL = Pp->data->blb;
int fngfs = Pp->data->fngfs;
while (BL)
{
Block *cg = BL->data;
if (myrank == cg->rank)
{
f_admmass_bssn_ss(cg->shape, cg->X[0], cg->X[1], cg->X[2],
cg->fgfs[fngfs + ShellPatch::gx], cg->fgfs[fngfs + ShellPatch::gy], cg->fgfs[fngfs + ShellPatch::gz],
cg->fgfs[fngfs + ShellPatch::drhodx], cg->fgfs[fngfs + ShellPatch::drhody], cg->fgfs[fngfs + ShellPatch::drhodz],
cg->fgfs[fngfs + ShellPatch::dsigmadx], cg->fgfs[fngfs + ShellPatch::dsigmady], cg->fgfs[fngfs + ShellPatch::dsigmadz],
cg->fgfs[fngfs + ShellPatch::dRdx], cg->fgfs[fngfs + ShellPatch::dRdy], cg->fgfs[fngfs + ShellPatch::dRdz],
cg->fgfs[fngfs + ShellPatch::drhodxx], cg->fgfs[fngfs + ShellPatch::drhodxy], cg->fgfs[fngfs + ShellPatch::drhodxz],
cg->fgfs[fngfs + ShellPatch::drhodyy], cg->fgfs[fngfs + ShellPatch::drhodyz], cg->fgfs[fngfs + ShellPatch::drhodzz],
cg->fgfs[fngfs + ShellPatch::dsigmadxx], cg->fgfs[fngfs + ShellPatch::dsigmadxy], cg->fgfs[fngfs + ShellPatch::dsigmadxz],
cg->fgfs[fngfs + ShellPatch::dsigmadyy], cg->fgfs[fngfs + ShellPatch::dsigmadyz], cg->fgfs[fngfs + ShellPatch::dsigmadzz],
cg->fgfs[fngfs + ShellPatch::dRdxx], cg->fgfs[fngfs + ShellPatch::dRdxy], cg->fgfs[fngfs + ShellPatch::dRdxz],
cg->fgfs[fngfs + ShellPatch::dRdyy], cg->fgfs[fngfs + ShellPatch::dRdyz], cg->fgfs[fngfs + ShellPatch::dRdzz],
cg->fgfs[chi->sgfn], cg->fgfs[trK->sgfn],
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn],
cg->fgfs[Gmx->sgfn], cg->fgfs[Gmy->sgfn], cg->fgfs[Gmz->sgfn],
cg->fgfs[Sfx_rhs->sgfn], cg->fgfs[Sfy_rhs->sgfn], cg->fgfs[Sfz_rhs->sgfn],
Symmetry, Pp->data->sst);
}
if (BL == Pp->data->ble)
break;
BL = BL->next;
}
Pp = Pp->next;
}
}
const int InList = 19;
const int idx_rpsi4 = 0, idx_ipsi4 = 1;
const int idx_sfx = 2, idx_sfy = 3, idx_sfz = 4;
const int idx_chi = 5, idx_trk = 6;
const int idx_gxx = 7, idx_gxy = 8, idx_gxz = 9, idx_gyy = 10, idx_gyz = 11, idx_gzz = 12;
const int idx_axx = 13, idx_axy = 14, idx_axz = 15, idx_ayy = 16, idx_ayz = 17, idx_azz = 18;
MyList<var> *DG_List = new MyList<var>(Rpsi4);
DG_List->insert(Ipsi4);
DG_List->insert(Sfx_rhs);
DG_List->insert(Sfy_rhs);
DG_List->insert(Sfz_rhs);
DG_List->insert(chi);
DG_List->insert(trK);
DG_List->insert(gxx);
DG_List->insert(gxy);
DG_List->insert(gxz);
DG_List->insert(gyy);
DG_List->insert(gyz);
DG_List->insert(gzz);
DG_List->insert(Axx);
DG_List->insert(Axy);
DG_List->insert(Axz);
DG_List->insert(Ayy);
DG_List->insert(Ayz);
DG_List->insert(Azz);
int n;
double *pox[3];
for (int ia = 0; ia < 3; ia++)
pox[ia] = new double[n_tot];
for (n = 0; n < n_tot; n++)
{
pox[0][n] = rex * nx_g[n];
pox[1][n] = rex * ny_g[n];
pox[2][n] = rex * nz_g[n];
}
double *shellf = new double[n_tot * InList];
GH->Interp_Points(DG_List, n_tot, pox, shellf, Symmetry);
int mp, Lp, Nmin, Nmax;
mp = n_tot / cpusize;
Lp = n_tot - cpusize * mp;
if (Lp > myrank)
{
Nmin = myrank * mp + myrank;
Nmax = Nmin + mp;
}
else
{
Nmin = myrank * mp + Lp;
Nmax = Nmin + mp - 1;
}
double *RP_out = new double[NN];
double *IP_out = new double[NN];
for (int ii = 0; ii < NN; ii++)
{
RP_out[ii] = 0.0;
IP_out[ii] = 0.0;
}
double Mass_out = 0.0;
double ang_outx = 0.0, ang_outy = 0.0, ang_outz = 0.0;
double p_outx = 0.0, p_outy = 0.0, p_outz = 0.0;
const double f1o8 = 0.125;
build_wave_cache(spinw, maxl);
for (n = Nmin; n <= Nmax; n++)
{
const int base = InList * n;
const int i = int(n / N_phi);
const int j = n - i * N_phi;
const double psi4RR0 = shellf[base + idx_rpsi4];
const double psi4II0 = shellf[base + idx_ipsi4];
const double psi4RR1 = Rpsi4->SoA[2] * psi4RR0;
const double psi4II1 = Ipsi4->SoA[2] * psi4II0;
for (int countlm = 0; countlm < wave_cache_modes; countlm++)
{
const int theta_idx = countlm * N_theta + i;
const int phi_idx = countlm * N_phi + j;
const double theta_pos = wave_theta_pos[theta_idx];
const double cosmphi_here = wave_phi_cos[phi_idx];
const double sinmphi_here = wave_phi_sin[phi_idx];
RP_out[countlm] += theta_pos * (psi4RR0 * cosmphi_here + psi4II0 * sinmphi_here);
IP_out[countlm] += theta_pos * (psi4II0 * cosmphi_here - psi4RR0 * sinmphi_here);
if (Symmetry == 1)
{
const double theta_neg = wave_theta_neg[theta_idx];
RP_out[countlm] += theta_neg * (psi4RR1 * cosmphi_here + psi4II1 * sinmphi_here);
IP_out[countlm] += theta_neg * (psi4II1 * cosmphi_here - psi4RR1 * sinmphi_here);
}
}
double Chi = shellf[base + idx_chi];
double TRK = shellf[base + idx_trk];
double Gxx = shellf[base + idx_gxx] + 1.0;
double Gxy = shellf[base + idx_gxy];
double Gxz = shellf[base + idx_gxz];
double Gyy = shellf[base + idx_gyy] + 1.0;
double Gyz = shellf[base + idx_gyz];
double Gzz = shellf[base + idx_gzz] + 1.0;
double axx = shellf[base + idx_axx];
double axy = shellf[base + idx_axy];
double axz = shellf[base + idx_axz];
double ayy = shellf[base + idx_ayy];
double ayz = shellf[base + idx_ayz];
double azz = shellf[base + idx_azz];
Chi = 1.0 / (1.0 + Chi);
const double Psi = Chi * sqrt(Chi);
#ifdef GaussInt
const double theta_weight = wtcostheta[i];
Mass_out += (shellf[base + idx_sfx] * nx_g[n] + shellf[base + idx_sfy] * ny_g[n] + shellf[base + idx_sfz] * nz_g[n]) * theta_weight;
#else
const double theta_weight = 1.0;
Mass_out += shellf[base + idx_sfx] * nx_g[n] + shellf[base + idx_sfy] * ny_g[n] + shellf[base + idx_sfz] * nz_g[n];
#endif
double detg = Gxx * Gyy * Gzz + Gxy * Gyz * Gxz + Gxz * Gxy * Gyz -
Gxz * Gyy * Gxz - Gxy * Gxy * Gzz - Gxx * Gyz * Gyz;
const double gupxx = (Gyy * Gzz - Gyz * Gyz) / detg;
const double gupxy = -(Gxy * Gzz - Gyz * Gxz) / detg;
const double gupxz = (Gxy * Gyz - Gyy * Gxz) / detg;
const double gupyy = (Gxx * Gzz - Gxz * Gxz) / detg;
const double gupyz = -(Gxx * Gyz - Gxy * Gxz) / detg;
const double gupzz = (Gxx * Gyy - Gxy * Gxy) / detg;
const double aupxx = gupxx * axx + gupxy * axy + gupxz * axz;
const double aupxy = gupxx * axy + gupxy * ayy + gupxz * ayz;
const double aupxz = gupxx * axz + gupxy * ayz + gupxz * azz;
const double aupyx = gupxy * axx + gupyy * axy + gupyz * axz;
const double aupyy = gupxy * axy + gupyy * ayy + gupyz * ayz;
const double aupyz = gupxy * axz + gupyy * ayz + gupyz * azz;
const double aupzx = gupxz * axx + gupyz * axy + gupzz * axz;
const double aupzy = gupxz * axy + gupyz * ayy + gupzz * ayz;
const double aupzz = gupxz * axz + gupyz * ayz + gupzz * azz;
if (Symmetry == 0)
{
ang_outx += f1o8 * Psi * (nx_g[n] * (pox[1][n] * aupxz - pox[2][n] * aupxy) + ny_g[n] * (pox[1][n] * aupyz - pox[2][n] * aupyy) + nz_g[n] * (pox[1][n] * aupzz - pox[2][n] * aupzy)) * theta_weight;
ang_outy += f1o8 * Psi * (nx_g[n] * (pox[2][n] * aupxx - pox[0][n] * aupxz) + ny_g[n] * (pox[2][n] * aupyx - pox[0][n] * aupyz) + nz_g[n] * (pox[2][n] * aupzx - pox[0][n] * aupzz)) * theta_weight;
ang_outz += f1o8 * Psi * (nx_g[n] * (pox[0][n] * aupxy - pox[1][n] * aupxx) + ny_g[n] * (pox[0][n] * aupyy - pox[1][n] * aupyx) + nz_g[n] * (pox[0][n] * aupzy - pox[1][n] * aupzx)) * theta_weight;
}
else
{
ang_outz += f1o8 * Psi * (nx_g[n] * (pox[0][n] * aupxy - pox[1][n] * aupxx) + ny_g[n] * (pox[0][n] * aupyy - pox[1][n] * aupyx) + nz_g[n] * (pox[0][n] * aupzy - pox[1][n] * aupzx)) * theta_weight;
}
axx = Chi * (axx + Gxx * TRK / 3.0);
axy = Chi * (axy + Gxy * TRK / 3.0);
axz = Chi * (axz + Gxz * TRK / 3.0);
ayy = Chi * (ayy + Gyy * TRK / 3.0);
ayz = Chi * (ayz + Gyz * TRK / 3.0);
azz = Chi * (azz + Gzz * TRK / 3.0);
axx -= TRK;
ayy -= TRK;
azz -= TRK;
p_outx += f1o8 * Psi * (nx_g[n] * axx + ny_g[n] * axy + nz_g[n] * axz) * theta_weight;
p_outy += f1o8 * Psi * (nx_g[n] * axy + ny_g[n] * ayy + nz_g[n] * ayz) * theta_weight;
if (Symmetry == 0)
p_outz += f1o8 * Psi * (nx_g[n] * axz + ny_g[n] * ayz + nz_g[n] * azz) * theta_weight;
}
for (int ii = 0; ii < NN; ii++)
{
#ifdef GaussInt
RP_out[ii] = RP_out[ii] * rex * dphi;
IP_out[ii] = IP_out[ii] * rex * dphi;
#else
RP_out[ii] = RP_out[ii] * rex * dphi * dcostheta;
IP_out[ii] = IP_out[ii] * rex * dphi * dcostheta;
#endif
}
double mass, px, py, pz, sx, sy, sz;
{
double *reduce_out = new double[2 * NN + 7];
double *reduce_in = new double[2 * NN + 7];
memcpy(reduce_out, RP_out, NN * sizeof(double));
memcpy(reduce_out + NN, IP_out, NN * sizeof(double));
reduce_out[2 * NN + 0] = Mass_out;
reduce_out[2 * NN + 1] = ang_outx;
reduce_out[2 * NN + 2] = ang_outy;
reduce_out[2 * NN + 3] = ang_outz;
reduce_out[2 * NN + 4] = p_outx;
reduce_out[2 * NN + 5] = p_outy;
reduce_out[2 * NN + 6] = p_outz;
MPI_Allreduce(reduce_out, reduce_in, 2 * NN + 7, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
memcpy(RP, reduce_in, NN * sizeof(double));
memcpy(IP, reduce_in + NN, NN * sizeof(double));
mass = reduce_in[2 * NN + 0];
sx = reduce_in[2 * NN + 1];
sy = reduce_in[2 * NN + 2];
sz = reduce_in[2 * NN + 3];
px = reduce_in[2 * NN + 4];
py = reduce_in[2 * NN + 5];
pz = reduce_in[2 * NN + 6];
delete[] reduce_out;
delete[] reduce_in;
}
#ifdef GaussInt
mass = mass * rex * rex * dphi * factor;
sx = sx * rex * rex * dphi * (1.0 / PI) * factor;
sy = sy * rex * rex * dphi * (1.0 / PI) * factor;
sz = sz * rex * rex * dphi * (1.0 / PI) * factor;
px = px * rex * rex * dphi * (1.0 / PI) * factor;
py = py * rex * rex * dphi * (1.0 / PI) * factor;
pz = pz * rex * rex * dphi * (1.0 / PI) * factor;
#else
mass = mass * rex * rex * dphi * dcostheta * factor;
sx = sx * rex * rex * dphi * dcostheta * (1.0 / PI) * factor;
sy = sy * rex * rex * dphi * dcostheta * (1.0 / PI) * factor;
sz = sz * rex * rex * dphi * dcostheta * (1.0 / PI) * factor;
px = px * rex * rex * dphi * dcostheta * (1.0 / PI) * factor;
py = py * rex * rex * dphi * dcostheta * (1.0 / PI) * factor;
pz = pz * rex * rex * dphi * dcostheta * (1.0 / PI) * factor;
#endif
Rout[0] = mass;
Rout[1] = px;
Rout[2] = py;
Rout[3] = pz;
Rout[4] = sx;
Rout[5] = sy;
Rout[6] = sz;
delete[] pox[0];
delete[] pox[1];
delete[] pox[2];
delete[] shellf;
delete[] RP_out;
delete[] IP_out;
DG_List->clearList();
}
//|----------------------------------------------------------------
// do not discriminate box and shell
// for Gravitational wave specially symmetric case

16
AMSS_NCKU_source/surface_integral.h
View File

@@ -94,6 +94,22 @@ public:
var *Gmx, var *Gmy, var *Gmz,
var *Sfx_rhs, var *Sfy_rhs, var *Sfz_rhs,
double *Rout, monitor *Monitor, bool refresh_mass_fields = true);
void surf_WaveMassPAng(double rex, int lev, cgh *GH,
var *Rpsi4, var *Ipsi4, int spinw, int maxl, int NN, double *RP, double *IP,
var *chi, var *trK,
var *gxx, var *gxy, var *gxz, var *gyy, var *gyz, var *gzz,
var *Axx, var *Axy, var *Axz, var *Ayy, var *Ayz, var *Azz,
var *Gmx, var *Gmy, var *Gmz,
var *Sfx_rhs, var *Sfy_rhs, var *Sfz_rhs,
double *Rout, monitor *Monitor, bool refresh_mass_fields = true);
void surf_WaveMassPAng(double rex, int lev, ShellPatch *GH,
var *Rpsi4, var *Ipsi4, int spinw, int maxl, int NN, double *RP, double *IP,
var *chi, var *trK,
var *gxx, var *gxy, var *gxz, var *gyy, var *gyz, var *gzz,
var *Axx, var *Axy, var *Axz, var *Ayy, var *Ayz, var *Azz,
var *Gmx, var *Gmy, var *Gmz,
var *Sfx_rhs, var *Sfy_rhs, var *Sfz_rhs,
double *Rout, monitor *Monitor, bool refresh_mass_fields = true);
void surf_Wave(double rex, cgh *GH, ShellPatch *SH,
var *chi, var *trK,
var *gxx, var *gxy, var *gxz, var *gyy, var *gyz, var *gzz,