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

Cache wave extraction angular kernels

Browse Source 
(cherry picked from commit e4c25eb21f)
This commit is contained in:
CGH0S7
2026-04-13 12:40:20 +08:00
committed by ianchb
parent c589097618
commit cce8a44fc4
2 changed files with 357 additions and 217 deletions
Download Patch File Download Diff File Expand all files Collapse all files

543
AMSS_NCKU_source/surface_integral.C
View File

@@ -36,8 +36,15 @@ using namespace std;
//| Constructor //| Constructor
//|============================================================================ //|============================================================================
surface_integral::surface_integral(int iSymmetry) : Symmetry(iSymmetry) surface_integral::surface_integral(int iSymmetry) : Symmetry(iSymmetry),
{ wave_cache_spinw(-1),
wave_cache_maxl(-1),
wave_cache_modes(0),
wave_theta_pos(0),
wave_theta_neg(0),
wave_phi_cos(0),
wave_phi_sin(0)
{
MPI_Comm_rank(MPI_COMM_WORLD, &myrank); MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
MPI_Comm_size(MPI_COMM_WORLD, &cpusize); MPI_Comm_size(MPI_COMM_WORLD, &cpusize);
int N = 40; int N = 40;
@@ -180,20 +187,80 @@ surface_integral::surface_integral(int iSymmetry) : Symmetry(iSymmetry)
//|============================================================================ //|============================================================================
//| Destructor //| Destructor
//|============================================================================ //|============================================================================
surface_integral::~surface_integral() surface_integral::~surface_integral()
{ {
delete[] nx_g; clear_wave_cache();
delete[] ny_g; delete[] nx_g;
delete[] nz_g; delete[] ny_g;
delete[] arcostheta; delete[] nz_g;
#ifdef GaussInt delete[] arcostheta;
delete[] wtcostheta; #ifdef GaussInt
#endif delete[] wtcostheta;
} #endif
//|---------------------------------------------------------------- }
// spin weighted spinw component of psi4, general routine void surface_integral::clear_wave_cache()
// l takes from spinw to maxl; m takes from -l to l {
//|---------------------------------------------------------------- delete[] wave_theta_pos;
delete[] wave_theta_neg;
delete[] wave_phi_cos;
delete[] wave_phi_sin;
wave_theta_pos = 0;
wave_theta_neg = 0;
wave_phi_cos = 0;
wave_phi_sin = 0;
wave_cache_spinw = -1;
wave_cache_maxl = -1;
wave_cache_modes = 0;
}
void surface_integral::build_wave_cache(int spinw, int maxl)
{
if (wave_cache_spinw == spinw && wave_cache_maxl == maxl && wave_theta_pos && wave_theta_neg && wave_phi_cos && wave_phi_sin)
return;
clear_wave_cache();
int modes = 0;
for (int pl = spinw; pl < maxl + 1; pl++)
for (int pm = -pl; pm < pl + 1; pm++)
modes++;
wave_theta_pos = new double[modes * N_theta];
wave_theta_neg = new double[modes * N_theta];
wave_phi_cos = new double[modes * N_phi];
wave_phi_sin = new double[modes * N_phi];
int countlm = 0;
for (int pl = spinw; pl < maxl + 1; pl++)
for (int pm = -pl; pm < pl + 1; pm++)
{
const double prefactor = sqrt((2 * pl + 1.0) / 4.0 / PI);
for (int i = 0; i < N_theta; i++)
{
#ifdef GaussInt
const double weight = wtcostheta[i];
#else
const double weight = 1.0;
#endif
wave_theta_pos[countlm * N_theta + i] = prefactor * misc::Wigner_d_function(pl, pm, spinw, arcostheta[i]) * weight;
wave_theta_neg[countlm * N_theta + i] = prefactor * misc::Wigner_d_function(pl, pm, spinw, -arcostheta[i]) * weight;
}
for (int j = 0; j < N_phi; j++)
{
const double phase = pm * (j + 0.5) * dphi;
wave_phi_cos[countlm * N_phi + j] = cos(phase);
wave_phi_sin[countlm * N_phi + j] = sin(phase);
}
countlm++;
}
wave_cache_spinw = spinw;
wave_cache_maxl = maxl;
wave_cache_modes = modes;
}
//|----------------------------------------------------------------
// spin weighted spinw component of psi4, general routine
// l takes from spinw to maxl; m takes from -l to l
//|----------------------------------------------------------------
void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *Ipsi4, 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, int spinw, int maxl, int NN, double *RP, double *IP,
monitor *Monitor) // NN is the length of RP and IP monitor *Monitor) // NN is the length of RP and IP
@@ -254,100 +321,134 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *
double costheta, thetap; double costheta, thetap;
double cosmphi, sinmphi; double cosmphi, sinmphi;
int i, j; int i, j;
int lpsy = 0; int lpsy = 0;
if (Symmetry == 0) if (Symmetry == 0)
lpsy = 1; lpsy = 1;
else if (Symmetry == 1) else if (Symmetry == 1)
lpsy = 2; lpsy = 2;
else if (Symmetry == 2) else if (Symmetry == 2)
lpsy = 8; lpsy = 8;
double psi4RR, psi4II; double psi4RR, psi4II;
for (n = Nmin; n <= Nmax; n++) if (Symmetry == 0 || Symmetry == 1)
{ {
// need round off always build_wave_cache(spinw, maxl);
i = int(n / N_phi); // int(1.723) = 1, int(-1.732) = -1 for (n = Nmin; n <= Nmax; n++)
j = n - i * N_phi; {
i = int(n / N_phi);
int countlm = 0; j = n - i * N_phi;
for (int pl = spinw; pl < maxl + 1; pl++)
for (int pm = -pl; pm < pl + 1; pm++) const double psi4RR0 = shellf[InList * n];
{ const double psi4II0 = shellf[InList * n + 1];
for (int lp = 0; lp < lpsy; lp++) const double psi4RR1 = Rpsi4->SoA[2] * psi4RR0;
{ const double psi4II1 = Ipsi4->SoA[2] * psi4II0;
switch (lp) for (int countlm = 0; countlm < wave_cache_modes; countlm++)
{ {
case 0: //+++ (theta, phi) const int theta_idx = countlm * N_theta + i;
costheta = arcostheta[i]; const int phi_idx = countlm * N_phi + j;
cosmphi = cos(pm * (j + 0.5) * dphi); const double theta_pos = wave_theta_pos[theta_idx];
sinmphi = sin(pm * (j + 0.5) * dphi); const double cosmphi_here = wave_phi_cos[phi_idx];
psi4RR = shellf[InList * n]; const double sinmphi_here = wave_phi_sin[phi_idx];
psi4II = shellf[InList * n + 1];
break; RP_out[countlm] += theta_pos * (psi4RR0 * cosmphi_here + psi4II0 * sinmphi_here);
case 1: //++- (pi-theta, phi) IP_out[countlm] += theta_pos * (psi4II0 * cosmphi_here - psi4RR0 * sinmphi_here);
costheta = -arcostheta[i]; if (Symmetry == 1)
cosmphi = cos(pm * (j + 0.5) * dphi); {
sinmphi = sin(pm * (j + 0.5) * dphi); const double theta_neg = wave_theta_neg[theta_idx];
psi4RR = Rpsi4->SoA[2] * shellf[InList * n]; RP_out[countlm] += theta_neg * (psi4RR1 * cosmphi_here + psi4II1 * sinmphi_here);
psi4II = Ipsi4->SoA[2] * shellf[InList * n + 1]; IP_out[countlm] += theta_neg * (psi4II1 * cosmphi_here - psi4RR1 * sinmphi_here);
break; }
case 2: //+-+ (theta, 2*pi-phi) }
costheta = arcostheta[i]; }
cosmphi = cos(pm * (j + 0.5) * dphi); }
sinmphi = -sin(pm * (j + 0.5) * dphi); else
psi4RR = Rpsi4->SoA[1] * shellf[InList * n]; {
psi4II = Ipsi4->SoA[1] * shellf[InList * n + 1]; for (n = Nmin; n <= Nmax; n++)
break; {
case 3: //+-- (pi-theta, 2*pi-phi) // need round off always
costheta = -arcostheta[i]; i = int(n / N_phi); // int(1.723) = 1, int(-1.732) = -1
cosmphi = cos(pm * (j + 0.5) * dphi); j = n - i * N_phi;
sinmphi = -sin(pm * (j + 0.5) * dphi);
psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * shellf[InList * n]; int countlm = 0;
psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * shellf[InList * n + 1]; for (int pl = spinw; pl < maxl + 1; pl++)
break; for (int pm = -pl; pm < pl + 1; pm++)
case 4: //-++ (theta, pi-phi) {
costheta = arcostheta[i]; for (int lp = 0; lp < lpsy; lp++)
cosmphi = cos(pm * (PI - (j + 0.5) * dphi)); {
sinmphi = sin(pm * (PI - (j + 0.5) * dphi)); switch (lp)
psi4RR = Rpsi4->SoA[0] * shellf[InList * n]; {
psi4II = Ipsi4->SoA[0] * shellf[InList * n + 1]; case 0: //+++ (theta, phi)
break; costheta = arcostheta[i];
case 5: //-+- (pi-theta, pi-phi) cosmphi = cos(pm * (j + 0.5) * dphi);
costheta = -arcostheta[i]; sinmphi = sin(pm * (j + 0.5) * dphi);
cosmphi = cos(pm * (PI - (j + 0.5) * dphi)); psi4RR = shellf[InList * n];
sinmphi = sin(pm * (PI - (j + 0.5) * dphi)); psi4II = shellf[InList * n + 1];
psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[0] * shellf[InList * n]; break;
psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[0] * shellf[InList * n + 1]; case 1: //++- (pi-theta, phi)
break; costheta = -arcostheta[i];
case 6: //--+ (theta, pi+phi) cosmphi = cos(pm * (j + 0.5) * dphi);
costheta = arcostheta[i]; sinmphi = sin(pm * (j + 0.5) * dphi);
cosmphi = cos(pm * (PI + (j + 0.5) * dphi)); psi4RR = Rpsi4->SoA[2] * shellf[InList * n];
sinmphi = sin(pm * (PI + (j + 0.5) * dphi)); psi4II = Ipsi4->SoA[2] * shellf[InList * n + 1];
psi4RR = Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n]; break;
psi4II = Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1]; case 2: //+-+ (theta, 2*pi-phi)
break; costheta = arcostheta[i];
case 7: //--- (pi-theta, pi+phi) cosmphi = cos(pm * (j + 0.5) * dphi);
costheta = -arcostheta[i]; sinmphi = -sin(pm * (j + 0.5) * dphi);
cosmphi = cos(pm * (PI + (j + 0.5) * dphi)); psi4RR = Rpsi4->SoA[1] * shellf[InList * n];
sinmphi = sin(pm * (PI + (j + 0.5) * dphi)); psi4II = Ipsi4->SoA[1] * shellf[InList * n + 1];
psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n]; break;
psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1]; case 3: //+-- (pi-theta, 2*pi-phi)
} costheta = -arcostheta[i];
cosmphi = cos(pm * (j + 0.5) * dphi);
thetap = sqrt((2 * pl + 1.0) / 4.0 / PI) * misc::Wigner_d_function(pl, pm, spinw, costheta); // note the variation from -2 to 2 sinmphi = -sin(pm * (j + 0.5) * dphi);
#ifdef GaussInt psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * shellf[InList * n];
// wtcostheta is even function respect costheta psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * shellf[InList * n + 1];
RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi) * wtcostheta[i]; break;
IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi) * wtcostheta[i]; case 4: //-++ (theta, pi-phi)
#else costheta = arcostheta[i];
RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi); cosmphi = cos(pm * (PI - (j + 0.5) * dphi));
IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi); sinmphi = sin(pm * (PI - (j + 0.5) * dphi));
#endif psi4RR = Rpsi4->SoA[0] * shellf[InList * n];
} psi4II = Ipsi4->SoA[0] * shellf[InList * n + 1];
countlm++; // no sanity check for countlm and NN which should be noted in the input parameters break;
} case 5: //-+- (pi-theta, pi-phi)
} costheta = -arcostheta[i];
cosmphi = cos(pm * (PI - (j + 0.5) * dphi));
sinmphi = sin(pm * (PI - (j + 0.5) * dphi));
psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[0] * shellf[InList * n];
psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[0] * shellf[InList * n + 1];
break;
case 6: //--+ (theta, pi+phi)
costheta = arcostheta[i];
cosmphi = cos(pm * (PI + (j + 0.5) * dphi));
sinmphi = sin(pm * (PI + (j + 0.5) * dphi));
psi4RR = Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n];
psi4II = Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1];
break;
case 7: //--- (pi-theta, pi+phi)
costheta = -arcostheta[i];
cosmphi = cos(pm * (PI + (j + 0.5) * dphi));
sinmphi = sin(pm * (PI + (j + 0.5) * dphi));
psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n];
psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1];
}
thetap = sqrt((2 * pl + 1.0) / 4.0 / PI) * misc::Wigner_d_function(pl, pm, spinw, costheta); // note the variation from -2 to 2
#ifdef GaussInt
// wtcostheta is even function respect costheta
RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi) * wtcostheta[i];
IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi) * wtcostheta[i];
#else
RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi);
IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi);
#endif
}
countlm++; // no sanity check for countlm and NN which should be noted in the input parameters
}
}
}
for (int ii = 0; ii < NN; ii++) for (int ii = 0; ii < NN; ii++)
{ {
@@ -456,100 +557,134 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *
double costheta, thetap; double costheta, thetap;
double cosmphi, sinmphi; double cosmphi, sinmphi;
int i, j; int i, j;
int lpsy = 0; int lpsy = 0;
if (Symmetry == 0) if (Symmetry == 0)
lpsy = 1; lpsy = 1;
else if (Symmetry == 1) else if (Symmetry == 1)
lpsy = 2; lpsy = 2;
else if (Symmetry == 2) else if (Symmetry == 2)
lpsy = 8; lpsy = 8;
double psi4RR, psi4II; double psi4RR, psi4II;
for (n = Nmin; n <= Nmax; n++) if (Symmetry == 0 || Symmetry == 1)
{ {
// need round off always build_wave_cache(spinw, maxl);
i = int(n / N_phi); // int(1.723) = 1, int(-1.732) = -1 for (n = Nmin; n <= Nmax; n++)
j = n - i * N_phi; {
i = int(n / N_phi);
int countlm = 0; j = n - i * N_phi;
for (int pl = spinw; pl < maxl + 1; pl++)
for (int pm = -pl; pm < pl + 1; pm++) const double psi4RR0 = shellf[InList * n];
{ const double psi4II0 = shellf[InList * n + 1];
for (int lp = 0; lp < lpsy; lp++) const double psi4RR1 = Rpsi4->SoA[2] * psi4RR0;
{ const double psi4II1 = Ipsi4->SoA[2] * psi4II0;
switch (lp) for (int countlm = 0; countlm < wave_cache_modes; countlm++)
{ {
case 0: //+++ (theta, phi) const int theta_idx = countlm * N_theta + i;
costheta = arcostheta[i]; const int phi_idx = countlm * N_phi + j;
cosmphi = cos(pm * (j + 0.5) * dphi); const double theta_pos = wave_theta_pos[theta_idx];
sinmphi = sin(pm * (j + 0.5) * dphi); const double cosmphi_here = wave_phi_cos[phi_idx];
psi4RR = shellf[InList * n]; const double sinmphi_here = wave_phi_sin[phi_idx];
psi4II = shellf[InList * n + 1];
break; RP_out[countlm] += theta_pos * (psi4RR0 * cosmphi_here + psi4II0 * sinmphi_here);
case 1: //++- (pi-theta, phi) IP_out[countlm] += theta_pos * (psi4II0 * cosmphi_here - psi4RR0 * sinmphi_here);
costheta = -arcostheta[i]; if (Symmetry == 1)
cosmphi = cos(pm * (j + 0.5) * dphi); {
sinmphi = sin(pm * (j + 0.5) * dphi); const double theta_neg = wave_theta_neg[theta_idx];
psi4RR = Rpsi4->SoA[2] * shellf[InList * n]; RP_out[countlm] += theta_neg * (psi4RR1 * cosmphi_here + psi4II1 * sinmphi_here);
psi4II = Ipsi4->SoA[2] * shellf[InList * n + 1]; IP_out[countlm] += theta_neg * (psi4II1 * cosmphi_here - psi4RR1 * sinmphi_here);
break; }
case 2: //+-+ (theta, 2*pi-phi) }
costheta = arcostheta[i]; }
cosmphi = cos(pm * (j + 0.5) * dphi); }
sinmphi = -sin(pm * (j + 0.5) * dphi); else
psi4RR = Rpsi4->SoA[1] * shellf[InList * n]; {
psi4II = Ipsi4->SoA[1] * shellf[InList * n + 1]; for (n = Nmin; n <= Nmax; n++)
break; {
case 3: //+-- (pi-theta, 2*pi-phi) // need round off always
costheta = -arcostheta[i]; i = int(n / N_phi); // int(1.723) = 1, int(-1.732) = -1
cosmphi = cos(pm * (j + 0.5) * dphi); j = n - i * N_phi;
sinmphi = -sin(pm * (j + 0.5) * dphi);
psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * shellf[InList * n]; int countlm = 0;
psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * shellf[InList * n + 1]; for (int pl = spinw; pl < maxl + 1; pl++)
break; for (int pm = -pl; pm < pl + 1; pm++)
case 4: //-++ (theta, pi-phi) {
costheta = arcostheta[i]; for (int lp = 0; lp < lpsy; lp++)
cosmphi = cos(pm * (PI - (j + 0.5) * dphi)); {
sinmphi = sin(pm * (PI - (j + 0.5) * dphi)); switch (lp)
psi4RR = Rpsi4->SoA[0] * shellf[InList * n]; {
psi4II = Ipsi4->SoA[0] * shellf[InList * n + 1]; case 0: //+++ (theta, phi)
break; costheta = arcostheta[i];
case 5: //-+- (pi-theta, pi-phi) cosmphi = cos(pm * (j + 0.5) * dphi);
costheta = -arcostheta[i]; sinmphi = sin(pm * (j + 0.5) * dphi);
cosmphi = cos(pm * (PI - (j + 0.5) * dphi)); psi4RR = shellf[InList * n];
sinmphi = sin(pm * (PI - (j + 0.5) * dphi)); psi4II = shellf[InList * n + 1];
psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[0] * shellf[InList * n]; break;
psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[0] * shellf[InList * n + 1]; case 1: //++- (pi-theta, phi)
break; costheta = -arcostheta[i];
case 6: //--+ (theta, pi+phi) cosmphi = cos(pm * (j + 0.5) * dphi);
costheta = arcostheta[i]; sinmphi = sin(pm * (j + 0.5) * dphi);
cosmphi = cos(pm * (PI + (j + 0.5) * dphi)); psi4RR = Rpsi4->SoA[2] * shellf[InList * n];
sinmphi = sin(pm * (PI + (j + 0.5) * dphi)); psi4II = Ipsi4->SoA[2] * shellf[InList * n + 1];
psi4RR = Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n]; break;
psi4II = Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1]; case 2: //+-+ (theta, 2*pi-phi)
break; costheta = arcostheta[i];
case 7: //--- (pi-theta, pi+phi) cosmphi = cos(pm * (j + 0.5) * dphi);
costheta = -arcostheta[i]; sinmphi = -sin(pm * (j + 0.5) * dphi);
cosmphi = cos(pm * (PI + (j + 0.5) * dphi)); psi4RR = Rpsi4->SoA[1] * shellf[InList * n];
sinmphi = sin(pm * (PI + (j + 0.5) * dphi)); psi4II = Ipsi4->SoA[1] * shellf[InList * n + 1];
psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n]; break;
psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1]; case 3: //+-- (pi-theta, 2*pi-phi)
} costheta = -arcostheta[i];
cosmphi = cos(pm * (j + 0.5) * dphi);
thetap = sqrt((2 * pl + 1.0) / 4.0 / PI) * misc::Wigner_d_function(pl, pm, spinw, costheta); // note the variation from -2 to 2 sinmphi = -sin(pm * (j + 0.5) * dphi);
#ifdef GaussInt psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * shellf[InList * n];
// wtcostheta is even function respect costheta psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * shellf[InList * n + 1];
RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi) * wtcostheta[i]; break;
IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi) * wtcostheta[i]; case 4: //-++ (theta, pi-phi)
#else costheta = arcostheta[i];
RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi); cosmphi = cos(pm * (PI - (j + 0.5) * dphi));
IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi); sinmphi = sin(pm * (PI - (j + 0.5) * dphi));
#endif psi4RR = Rpsi4->SoA[0] * shellf[InList * n];
} psi4II = Ipsi4->SoA[0] * shellf[InList * n + 1];
countlm++; // no sanity check for countlm and NN which should be noted in the input parameters break;
} case 5: //-+- (pi-theta, pi-phi)
} costheta = -arcostheta[i];
cosmphi = cos(pm * (PI - (j + 0.5) * dphi));
sinmphi = sin(pm * (PI - (j + 0.5) * dphi));
psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[0] * shellf[InList * n];
psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[0] * shellf[InList * n + 1];
break;
case 6: //--+ (theta, pi+phi)
costheta = arcostheta[i];
cosmphi = cos(pm * (PI + (j + 0.5) * dphi));
sinmphi = sin(pm * (PI + (j + 0.5) * dphi));
psi4RR = Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n];
psi4II = Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1];
break;
case 7: //--- (pi-theta, pi+phi)
costheta = -arcostheta[i];
cosmphi = cos(pm * (PI + (j + 0.5) * dphi));
sinmphi = sin(pm * (PI + (j + 0.5) * dphi));
psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n];
psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1];
}
thetap = sqrt((2 * pl + 1.0) / 4.0 / PI) * misc::Wigner_d_function(pl, pm, spinw, costheta); // note the variation from -2 to 2
#ifdef GaussInt
// wtcostheta is even function respect costheta
RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi) * wtcostheta[i];
IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi) * wtcostheta[i];
#else
RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi);
IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi);
#endif
}
countlm++; // no sanity check for countlm and NN which should be noted in the input parameters
}
}
}
for (int ii = 0; ii < NN; ii++) for (int ii = 0; ii < NN; ii++)
{ {

31
AMSS_NCKU_source/surface_integral.h
View File

@@ -27,19 +27,24 @@ using namespace std;
class surface_integral class surface_integral
{ {
private: private:
int Symmetry, factor; int Symmetry, factor;
int N_theta, N_phi; // Number of points in Theta & Phi directions int N_theta, N_phi; // Number of points in Theta & Phi directions
double dphi, dcostheta; double dphi, dcostheta;
double *arcostheta, *wtcostheta; double *arcostheta, *wtcostheta;
int n_tot; // size of arrays int n_tot; // size of arrays
double *nx_g, *ny_g, *nz_g; // global list of unit normals double *nx_g, *ny_g, *nz_g; // global list of unit normals
int myrank, cpusize; int myrank, cpusize;
int wave_cache_spinw, wave_cache_maxl, wave_cache_modes;
public: double *wave_theta_pos, *wave_theta_neg;
surface_integral(int iSymmetry); double *wave_phi_cos, *wave_phi_sin;
~surface_integral(); void clear_wave_cache();
void build_wave_cache(int spinw, int maxl);
public:
surface_integral(int iSymmetry);
~surface_integral();
void surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *Ipsi4, void surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *Ipsi4,
int spinw, int maxl, int NN, double *RP, double *IP, int spinw, int maxl, int NN, double *RP, double *IP,