From 11977eb82f1791c7eb8c4275f401dc25d9f52460 Mon Sep 17 00:00:00 2001 From: CGH0S7 <776459475@qq.com> Date: Mon, 13 Apr 2026 13:17:36 +0800 Subject: [PATCH] Merge wave and mass extraction interpolation (cherry picked from commit f3988ac8cac9f92200b22eee2b2a793ffe29cecc) --- AMSS_NCKU_source/bssn_class.C | 39 +- AMSS_NCKU_source/surface_integral.C | 962 +++++++++++++++++++++++----- AMSS_NCKU_source/surface_integral.h | 22 +- 3 files changed, 831 insertions(+), 192 deletions(-) diff --git a/AMSS_NCKU_source/bssn_class.C b/AMSS_NCKU_source/bssn_class.C index cc1f83f..1e9c79d 100644 --- a/AMSS_NCKU_source/bssn_class.C +++ b/AMSS_NCKU_source/bssn_class.C @@ -7678,32 +7678,35 @@ 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, - 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); + 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, + 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); + 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, + 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); + 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, + 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]); diff --git a/AMSS_NCKU_source/surface_integral.C b/AMSS_NCKU_source/surface_integral.C index 1b92167..1e58dce 100644 --- a/AMSS_NCKU_source/surface_integral.C +++ b/AMSS_NCKU_source/surface_integral.C @@ -45,8 +45,8 @@ surface_integral::surface_integral(int iSymmetry) : Symmetry(iSymmetry), wave_phi_cos(0), wave_phi_sin(0) { - MPI_Comm_rank(MPI_COMM_WORLD, &myrank); - MPI_Comm_size(MPI_COMM_WORLD, &cpusize); + MPI_Comm_rank(MPI_COMM_WORLD, &myrank); + MPI_Comm_size(MPI_COMM_WORLD, &cpusize); int N = 40; // read parameter from file { @@ -271,24 +271,24 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var * else cout << "WARNING: surface integral on multipatches" << endl; - const int InList = 2; - - MyList *DG_List = new MyList(Rpsi4); - DG_List->insert(Ipsi4); - - int n; - double *pox[3]; - for (int i = 0; i < 3; i++) - pox[i] = 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; + const int InList = 2; + + MyList *DG_List = new MyList(Rpsi4); + DG_List->insert(Ipsi4); + + int n; + double *pox[3]; + for (int i = 0; i < 3; i++) + pox[i] = 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) { @@ -301,16 +301,16 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var * Nmax = Nmin + mp - 1; } - double *shellf; - shellf = new double[n_tot * InList]; - - GH->PatL[lev]->data->Interp_Points(DG_List, n_tot, pox, shellf, Symmetry, Nmin, Nmax); + double *shellf; + shellf = new double[n_tot * InList]; + + GH->PatL[lev]->data->Interp_Points(DG_List, n_tot, pox, shellf, Symmetry, Nmin, Nmax); //|~~~~~> Integrate the dot product of Dphi with the surface normal. - double *RP_out, *IP_out; - RP_out = new double[NN]; - IP_out = new double[NN]; + double *RP_out, *IP_out; + RP_out = new double[NN]; + IP_out = new double[NN]; for (int ii = 0; ii < NN; ii++) { @@ -462,28 +462,28 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var * } //|------+ Communicate and sum the results from each processor. - { - double *RPIP_out = new double[2 * NN]; - double *RPIP = new double[2 * NN]; - memcpy(RPIP_out, RP_out, NN * sizeof(double)); - memcpy(RPIP_out + NN, IP_out, NN * sizeof(double)); - MPI_Allreduce(RPIP_out, RPIP, 2 * NN, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - memcpy(RP, RPIP, NN * sizeof(double)); - memcpy(IP, RPIP + NN, NN * sizeof(double)); - delete[] RPIP_out; - delete[] RPIP; - } - - //|------= Free memory. - - delete[] pox[0]; - delete[] pox[1]; - delete[] pox[2]; - delete[] shellf; - delete[] RP_out; - delete[] IP_out; - DG_List->clearList(); -} + { + double *RPIP_out = new double[2 * NN]; + double *RPIP = new double[2 * NN]; + memcpy(RPIP_out, RP_out, NN * sizeof(double)); + memcpy(RPIP_out + NN, IP_out, NN * sizeof(double)); + MPI_Allreduce(RPIP_out, RPIP, 2 * NN, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + memcpy(RP, RPIP, NN * sizeof(double)); + memcpy(IP, RPIP + NN, NN * sizeof(double)); + delete[] RPIP_out; + delete[] RPIP; + } + + //|------= Free memory. + + delete[] pox[0]; + delete[] pox[1]; + delete[] pox[2]; + delete[] shellf; + delete[] RP_out; + delete[] IP_out; + DG_List->clearList(); +} void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *Ipsi4, int spinw, int maxl, int NN, double *RP, double *IP, monitor *Monitor, MPI_Comm Comm_here) // NN is the length of RP and IP @@ -498,24 +498,24 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var * else cout << "WARNING: surface integral on multipatches" << endl; - const int InList = 2; - - MyList *DG_List = new MyList(Rpsi4); - DG_List->insert(Ipsi4); - - int n; - double *pox[3]; - for (int i = 0; i < 3; i++) - pox[i] = 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; - shellf = new double[n_tot * InList]; + const int InList = 2; + + MyList *DG_List = new MyList(Rpsi4); + DG_List->insert(Ipsi4); + + int n; + double *pox[3]; + for (int i = 0; i < 3; i++) + pox[i] = 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; + shellf = new double[n_tot * InList]; // misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Interp_Points"); @@ -544,9 +544,9 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var * //|~~~~~> Integrate the dot product of Dphi with the surface normal. - double *RP_out, *IP_out; - RP_out = new double[NN]; - IP_out = new double[NN]; + double *RP_out, *IP_out; + RP_out = new double[NN]; + IP_out = new double[NN]; for (int ii = 0; ii < NN; ii++) { @@ -698,28 +698,28 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var * } //|------+ Communicate and sum the results from each processor. - { - double *RPIP_out = new double[2 * NN]; - double *RPIP = new double[2 * NN]; - memcpy(RPIP_out, RP_out, NN * sizeof(double)); - memcpy(RPIP_out + NN, IP_out, NN * sizeof(double)); - MPI_Allreduce(RPIP_out, RPIP, 2 * NN, MPI_DOUBLE, MPI_SUM, Comm_here); - memcpy(RP, RPIP, NN * sizeof(double)); - memcpy(IP, RPIP + NN, NN * sizeof(double)); - delete[] RPIP_out; - delete[] RPIP; - } - - //|------= Free memory. - - delete[] pox[0]; - delete[] pox[1]; - delete[] pox[2]; - delete[] shellf; - delete[] RP_out; - delete[] IP_out; - DG_List->clearList(); -} + { + double *RPIP_out = new double[2 * NN]; + double *RPIP = new double[2 * NN]; + memcpy(RPIP_out, RP_out, NN * sizeof(double)); + memcpy(RPIP_out + NN, IP_out, NN * sizeof(double)); + MPI_Allreduce(RPIP_out, RPIP, 2 * NN, MPI_DOUBLE, MPI_SUM, Comm_here); + memcpy(RP, RPIP, NN * sizeof(double)); + memcpy(IP, RPIP + NN, NN * sizeof(double)); + delete[] RPIP_out; + delete[] RPIP; + } + + //|------= Free memory. + + delete[] pox[0]; + delete[] pox[1]; + delete[] pox[2]; + delete[] shellf; + delete[] RP_out; + delete[] IP_out; + DG_List->clearList(); +} //|---------------------------------------------------------------- // for shell patch //|---------------------------------------------------------------- @@ -2507,23 +2507,23 @@ void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, var *chi, var 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 i = 0; i < 3; i++) - pox[i] = 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; + DG_List->insert(Ayy); + DG_List->insert(Ayz); + DG_List->insert(Azz); + + int n; + double *pox[3]; + for (int i = 0; i < 3; i++) + pox[i] = 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) { @@ -2536,12 +2536,12 @@ void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, var *chi, var Nmax = Nmin + mp - 1; } - double *shellf; - shellf = new double[n_tot * InList]; - - // we have assumed there is only one box on this level, - // so we do not need loop boxes - GH->PatL[lev]->data->Interp_Points(DG_List, n_tot, pox, shellf, Symmetry, Nmin, Nmax); + double *shellf; + shellf = new double[n_tot * InList]; + + // we have assumed there is only one box on this level, + // so we do not need loop boxes + GH->PatL[lev]->data->Interp_Points(DG_List, n_tot, pox, shellf, Symmetry, Nmin, Nmax); double Mass_out = 0; double ang_outx, ang_outy, ang_outz; @@ -2611,7 +2611,7 @@ void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, var *chi, var #ifdef GaussInt // wtcostheta is even function respect costheta // 1/8\pi \int \psi^6 (y A^m_z - zA^m_y) dS_m - ang_outx = 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)) * wtcostheta[i]; + ang_outx = 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)) * wtcostheta[i]; // 1/8\pi \int \psi^6 (z A^m_x - xA^m_z) dS_m ang_outy = 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)) * wtcostheta[i]; // 1/8\pi \int \psi^6 (x A^m_y - yA^m_x) dS_m @@ -2708,12 +2708,12 @@ void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, var *chi, var Rout[5] = sy; Rout[6] = sz; - delete[] pox[0]; - delete[] pox[1]; - delete[] pox[2]; - delete[] shellf; - DG_List->clearList(); -} + delete[] pox[0]; + delete[] pox[1]; + delete[] pox[2]; + delete[] shellf; + DG_List->clearList(); +} void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, 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, @@ -2774,23 +2774,23 @@ void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, var *chi, var 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 i = 0; i < 3; i++) - pox[i] = 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; - shellf = new double[n_tot * InList]; + DG_List->insert(Ayy); + DG_List->insert(Ayz); + DG_List->insert(Azz); + + int n; + double *pox[3]; + for (int i = 0; i < 3; i++) + pox[i] = 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; + shellf = new double[n_tot * InList]; // we have assumed there is only one box on this level, // so we do not need loop boxes @@ -2980,12 +2980,12 @@ void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, var *chi, var Rout[5] = sy; Rout[6] = sz; - delete[] pox[0]; - delete[] pox[1]; - delete[] pox[2]; - delete[] shellf; - DG_List->clearList(); -} + delete[] pox[0]; + delete[] pox[1]; + delete[] pox[2]; + delete[] shellf; + DG_List->clearList(); +} //|---------------------------------------------------------------- // for shell patch //|---------------------------------------------------------------- @@ -3064,23 +3064,23 @@ void surface_integral::surf_MassPAng(double rex, int lev, ShellPatch *GH, var *c 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 i = 0; i < 3; i++) - pox[i] = 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; - shellf = new double[n_tot * InList]; + DG_List->insert(Ayy); + DG_List->insert(Ayz); + DG_List->insert(Azz); + + int n; + double *pox[3]; + for (int i = 0; i < 3; i++) + pox[i] = 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; + shellf = new double[n_tot * InList]; // we have assumed there is only one box on this level, // so we do not need loop boxes @@ -3266,15 +3266,635 @@ void surface_integral::surf_MassPAng(double rex, int lev, ShellPatch *GH, var *c Rout[5] = sy; Rout[6] = sz; - delete[] pox[0]; - delete[] pox[1]; - delete[] pox[2]; - delete[] shellf; - DG_List->clearList(); -} -//|---------------------------------------------------------------- -// do not discriminate box and shell -// for Gravitational wave specially symmetric case + delete[] pox[0]; + delete[] pox[1]; + delete[] pox[2]; + 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 *Pp = GH->PatL[lev]; + while (Pp) + { + MyList *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 *DG_List = new MyList(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 *Pp = GH->PatL; + while (Pp) + { + MyList *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 *DG_List = new MyList(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 //|---------------------------------------------------------------- void surface_integral::surf_Wave(double rex, cgh *GH, ShellPatch *SH, var *chi, var *trK, diff --git a/AMSS_NCKU_source/surface_integral.h b/AMSS_NCKU_source/surface_integral.h index e1e05e8..1b2d287 100644 --- a/AMSS_NCKU_source/surface_integral.h +++ b/AMSS_NCKU_source/surface_integral.h @@ -94,9 +94,25 @@ 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_Wave(double rex, cgh *GH, ShellPatch *SH, - var *chi, var *trK, - var *gxx, var *gxy, var *gxz, var *gyy, var *gyz, var *gzz, + 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, var *Axx, var *Axy, var *Axz, var *Ayy, var *Ayz, var *Azz, var *chix, var *chiy, var *chiz, var *trKx, var *trKy, var *trKz,