diff --git a/AMSS_NCKU_Input.py b/AMSS_NCKU_Input.py index 40d252b..e3cee28 100755 --- a/AMSS_NCKU_Input.py +++ b/AMSS_NCKU_Input.py @@ -31,7 +31,7 @@ GPU_Part = 0.0 ## Setting the physical system and numerical method Symmetry = "equatorial-symmetry" ## Symmetry of System: choose equatorial-symmetry、no-symmetry、octant-symmetry -Equation_Class = "Z4C" ## Evolution Equation: choose "BSSN", "BSSN-EScalar", "BSSN-EM", "Z4C" +Equation_Class = "BSSN-EScalar" ## Evolution Equation: choose "BSSN", "BSSN-EScalar", "BSSN-EM", "Z4C" ## If "BSSN-EScalar" is chosen, it is necessary to set other parameters below Initial_Data_Method = "Ansorg-TwoPuncture" ## initial data method: choose "Ansorg-TwoPuncture", "Lousto-Analytical", "Cao-Analytical", "KerrSchild-Analytical" Time_Evolution_Method = "runge-kutta-45" ## time evolution method: choose "runge-kutta-45" diff --git a/AMSS_NCKU_source/Parallel.C b/AMSS_NCKU_source/Parallel.C index 239574e..555acdd 100644 --- a/AMSS_NCKU_source/Parallel.C +++ b/AMSS_NCKU_source/Parallel.C @@ -341,7 +341,7 @@ bool cuda_state_count_direct_supported(int state_count) #if USE_CUDA_Z4C && (ABEtype == 2) return state_count == Z4C_CUDA_STATE_COUNT; #elif USE_CUDA_BSSN - return state_count > 0 && state_count <= BSSN_CUDA_STATE_COUNT; + return state_count == BSSN_CUDA_STATE_COUNT; #else (void)state_count; return false; @@ -393,6 +393,14 @@ bool cuda_can_direct_pack(const Parallel::gridseg *src, const Parallel::gridseg #elif USE_CUDA_BSSN if (bssn_cuda_has_resident_state(src->Bg) == 0) return false; + if (VarLists) + { + double *view_ptrs[BSSN_CUDA_STATE_COUNT]; + if (!cuda_build_bssn_host_views(src->Bg, VarLists, BSSN_CUDA_STATE_COUNT, view_ptrs)) + return false; + if (bssn_cuda_resident_state_matches(src->Bg, view_ptrs) == 0) + return false; + } if (type == 1) return true; int a[3], b[3]; @@ -427,7 +435,17 @@ bool cuda_can_direct_unpack(const Parallel::gridseg *dst, int type, MyList (void)VarListd; return true; #elif USE_CUDA_BSSN - return bssn_cuda_has_resident_state(dst->Bg) != 0; + if (bssn_cuda_has_resident_state(dst->Bg) == 0) + return false; + if (VarListd) + { + double *view_ptrs[BSSN_CUDA_STATE_COUNT]; + if (!cuda_build_bssn_host_views(dst->Bg, VarListd, BSSN_CUDA_STATE_COUNT, view_ptrs)) + return false; + if (bssn_cuda_resident_state_matches(dst->Bg, view_ptrs) == 0) + return false; + } + return true; #else return false; #endif @@ -443,7 +461,7 @@ bool cuda_direct_pack_segment(double *buffer, if (state_count != Z4C_CUDA_STATE_COUNT) return false; #elif USE_CUDA_BSSN - if (state_count <= 0 || state_count > BSSN_CUDA_STATE_COUNT) + if (state_count != BSSN_CUDA_STATE_COUNT) return false; #else return false; @@ -490,7 +508,7 @@ bool cuda_direct_unpack_segment(double *buffer, if (state_count != Z4C_CUDA_STATE_COUNT) return false; #elif USE_CUDA_BSSN - if (state_count <= 0 || state_count > BSSN_CUDA_STATE_COUNT) + if (state_count != BSSN_CUDA_STATE_COUNT) return false; #else return false; @@ -771,7 +789,7 @@ bool cuda_direct_pack_segment_to_device(double *buffer, } #endif #if USE_CUDA_BSSN - if (state_count <= 0 || state_count > BSSN_CUDA_STATE_COUNT) + if (state_count != BSSN_CUDA_STATE_COUNT) return false; const double t0 = sync_profile_enabled() ? MPI_Wtime() : 0.0; bool ok = false; @@ -963,7 +981,7 @@ bool cuda_direct_unpack_segment_from_device(double *buffer, } #endif #if USE_CUDA_BSSN - if (state_count <= 0 || state_count > BSSN_CUDA_STATE_COUNT) + if (state_count != BSSN_CUDA_STATE_COUNT) return false; const double t0 = sync_profile_enabled() ? MPI_Wtime() : 0.0; const int i0 = cuda_seg_begin(dst, dst->Bg, 0); @@ -1017,7 +1035,7 @@ bool cuda_download_resident_subset_to_host(Block *block, } #endif #if USE_CUDA_BSSN - if (!block || state_count <= 0 || state_count > BSSN_CUDA_STATE_COUNT) + if (!block || state_count != BSSN_CUDA_STATE_COUNT) return false; if (bssn_cuda_has_resident_state(block) == 0) return true; @@ -1032,6 +1050,8 @@ bool cuda_download_resident_subset_to_host(Block *block, views[i] = block->fgfs[v->data->sgfn]; v = v->next; } + if (bssn_cuda_resident_state_matches(block, views) == 0) + return false; return bssn_cuda_download_state_subset(block, block->shape, state_count, indices, views) == 0; #else (void)block; (void)vars; (void)state_count; @@ -1085,7 +1105,7 @@ bool cuda_device_state_count_supported(int state_count) return true; #endif #if USE_CUDA_BSSN - return state_count > 0 && state_count <= BSSN_CUDA_STATE_COUNT; + return state_count == BSSN_CUDA_STATE_COUNT; #else (void)state_count; return false; @@ -7259,6 +7279,8 @@ void Parallel::prepare_inter_time_level(Patch *Pat, cuda_build_bssn_host_views(cg, VarList2, state_count, src2_views) && cuda_build_bssn_host_views(cg, VarList3, state_count, dst_views) && bssn_cuda_has_resident_state(cg) && + bssn_cuda_resident_state_matches(cg, src1_views) && + bssn_cuda_resident_state_matches(cg, src2_views) && bssn_cuda_prepare_inter_time_level(cg, cg->shape, src1_views, src2_views, 0, dst_views, 2, tindex) == 0) @@ -7336,6 +7358,9 @@ void Parallel::prepare_inter_time_level(Patch *Pat, cuda_build_bssn_host_views(cg, VarList3, state_count, src3_views) && cuda_build_bssn_host_views(cg, VarList4, state_count, dst_views) && bssn_cuda_has_resident_state(cg) && + bssn_cuda_resident_state_matches(cg, src1_views) && + bssn_cuda_resident_state_matches(cg, src2_views) && + bssn_cuda_resident_state_matches(cg, src3_views) && bssn_cuda_prepare_inter_time_level(cg, cg->shape, src1_views, src2_views, src3_views, dst_views, 3, tindex) == 0) diff --git a/AMSS_NCKU_source/bssnEM_class.C b/AMSS_NCKU_source/bssnEM_class.C index 97bfbbd..8ebd9ac 100644 --- a/AMSS_NCKU_source/bssnEM_class.C +++ b/AMSS_NCKU_source/bssnEM_class.C @@ -15,10 +15,13 @@ using namespace std; #include "misc.h" #include "Ansorg.h" #include "fmisc.h" -#include "Parallel.h" -#include "bssnEM_class.h" -#include "bssn_rhs.h" -#include "empart.h" +#include "Parallel.h" +#include "bssnEM_class.h" +#include "bssn_rhs.h" +#if USE_CUDA_BSSN +#include "bssn_rhs_cuda.h" +#endif +#include "empart.h" #include "initial_puncture.h" #include "initial_maxwell.h" #include "enforce_algebra.h" @@ -32,11 +35,111 @@ using namespace std; #ifdef With_AHF #include "derivatives.h" #include "myglobal.h" -#endif - -//================================================================================================ - -// Define bssnEM_class +#endif + +//================================================================================================ + +#if USE_CUDA_BSSN +namespace { + +bool fill_bssn_cuda_views_prefix(Block *cg, MyList *vars, + double **host_views, + double *propspeeds = nullptr, + double *soa_flat = nullptr) +{ + int idx = 0; + while (vars && idx < BSSN_CUDA_STATE_COUNT) + { + host_views[idx] = cg->fgfs[vars->data->sgfn]; + if (propspeeds) + propspeeds[idx] = vars->data->propspeed; + if (soa_flat) + { + soa_flat[3 * idx + 0] = vars->data->SoA[0]; + soa_flat[3 * idx + 1] = vars->data->SoA[1]; + soa_flat[3 * idx + 2] = vars->data->SoA[2]; + } + vars = vars->next; + ++idx; + } + return idx == BSSN_CUDA_STATE_COUNT; +} + +void skip_bssn_cuda_prefix(MyList *&a, MyList *&b, MyList *&c) +{ + for (int i = 0; i < BSSN_CUDA_STATE_COUNT && a && b && c; ++i) + { + a = a->next; + b = b->next; + c = c->next; + } +} + +void skip_bssn_cuda_prefix(MyList *&a, MyList *&b, + MyList *&c, MyList *&d) +{ + for (int i = 0; i < BSSN_CUDA_STATE_COUNT && a && b && c && d; ++i) + { + a = a->next; + b = b->next; + c = c->next; + d = d->next; + } +} + +int run_bssn_em_cuda_substep(Block *cg, + MyList *state_in_list, + MyList *state_out_list, + Patch *patch, + double &dT_lev, + double &TRK4, + int &iter_count, + int &Symmetry, + int lev, + double &ndeps, + int &co, + double &chitiny, + var *rho, var *Sx, var *Sy, var *Sz, + var *Sxx, var *Sxy, var *Sxz, + var *Syy, var *Syz, var *Szz) +{ + double *state_in[BSSN_CUDA_STATE_COUNT]; + double *state_out[BSSN_CUDA_STATE_COUNT]; + double *matter[BSSN_CUDA_MATTER_COUNT] = { + cg->fgfs[rho->sgfn], cg->fgfs[Sx->sgfn], cg->fgfs[Sy->sgfn], cg->fgfs[Sz->sgfn], + cg->fgfs[Sxx->sgfn], cg->fgfs[Sxy->sgfn], cg->fgfs[Sxz->sgfn], + cg->fgfs[Syy->sgfn], cg->fgfs[Syz->sgfn], cg->fgfs[Szz->sgfn]}; + double propspeed[BSSN_CUDA_STATE_COUNT]; + double soa_flat[3 * BSSN_CUDA_STATE_COUNT]; + if (!fill_bssn_cuda_views_prefix(cg, state_in_list, state_in, propspeed, soa_flat) || + !fill_bssn_cuda_views_prefix(cg, state_out_list, state_out)) + return 1; + + int apply_bam_bc = 0; +#if (SommerType == 0) +#ifndef WithShell + apply_bam_bc = (lev == 0) ? 1 : 0; +#endif +#endif + int use_zero_matter = 0; + int keep_resident_state = 0; + int apply_enforce_ga = 0; + return bssn_cuda_rk4_substep(cg, + cg->shape, cg->X[0], cg->X[1], cg->X[2], + state_in, state_out, matter, + propspeed, soa_flat, patch->bbox, + dT_lev, TRK4, iter_count, apply_bam_bc, + Symmetry, lev, ndeps, co, + use_zero_matter, + keep_resident_state, apply_enforce_ga, chitiny); +} + +} +#endif + +//================================================================================================ + +// Define bssnEM_class // It inherits some members and methods from the parent class bssn_class and modifies others. // The modified members and methods are defined below (and in the header bssnEM_class.h). @@ -853,10 +956,11 @@ void bssnEM_class::Step(int lev, int YN) cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn], cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn]); -#endif - - if ( - f_compute_rhs_empart(cg->shape, cg->X[0], cg->X[1], cg->X[2], +#endif + + bool used_gpu_substep = false; + if ( + f_compute_rhs_empart(cg->shape, cg->X[0], cg->X[1], cg->X[2], cg->fgfs[phi0->sgfn], cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn], @@ -873,11 +977,20 @@ void bssnEM_class::Step(int lev, int YN) cg->fgfs[Kpsi_rhs->sgfn], cg->fgfs[Kphi_rhs->sgfn], cg->fgfs[rho->sgfn], cg->fgfs[Sx->sgfn], cg->fgfs[Sy->sgfn], cg->fgfs[Sz->sgfn], - cg->fgfs[Sxx->sgfn], cg->fgfs[Sxy->sgfn], cg->fgfs[Sxz->sgfn], - cg->fgfs[Syy->sgfn], cg->fgfs[Syz->sgfn], cg->fgfs[Szz->sgfn], - Symmetry, lev, ndeps) || - f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2], - cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn], + cg->fgfs[Sxx->sgfn], cg->fgfs[Sxy->sgfn], cg->fgfs[Sxz->sgfn], + cg->fgfs[Syy->sgfn], cg->fgfs[Syz->sgfn], cg->fgfs[Szz->sgfn], + Symmetry, lev, ndeps) || +#if USE_CUDA_BSSN + ((used_gpu_substep = + (run_bssn_em_cuda_substep(cg, StateList, SynchList_pre, Pp->data, + dT_lev, TRK4, iter_count, Symmetry, lev, + ndeps, pre, chitiny, + rho, Sx, Sy, Sz, Sxx, Sxy, Sxz, Syy, Syz, Szz) == 0)) + ? 0 + : 1) || +#endif + (!used_gpu_substep && f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2], + cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn], cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn], cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], @@ -906,10 +1019,10 @@ void bssnEM_class::Step(int lev, int YN) cg->fgfs[Gamzyy->sgfn], cg->fgfs[Gamzyz->sgfn], cg->fgfs[Gamzzz->sgfn], cg->fgfs[Rxx->sgfn], cg->fgfs[Rxy->sgfn], cg->fgfs[Rxz->sgfn], cg->fgfs[Ryy->sgfn], cg->fgfs[Ryz->sgfn], cg->fgfs[Rzz->sgfn], - cg->fgfs[Cons_Ham->sgfn], - cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn], - cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn], - Symmetry, lev, ndeps, pre)) + cg->fgfs[Cons_Ham->sgfn], + cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn], + cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn], + Symmetry, lev, ndeps, pre))) { cout << "find NaN in domain: (" << cg->bbox[0] << ":" << cg->bbox[3] << "," @@ -919,11 +1032,15 @@ void bssnEM_class::Step(int lev, int YN) } // rk4 substep and boundary - { - MyList *varl0 = StateList, *varl = SynchList_pre, *varlrhs = RHSList; - // we do not check the correspondence here - - while (varl0) + { + MyList *varl0 = StateList, *varl = SynchList_pre, *varlrhs = RHSList; + // we do not check the correspondence here +#if USE_CUDA_BSSN + if (used_gpu_substep) + skip_bssn_cuda_prefix(varl0, varl, varlrhs); +#endif + + while (varl0) { #ifndef WithShell if (lev == 0) // sommerfeld indeed @@ -1309,10 +1426,11 @@ void bssnEM_class::Step(int lev, int YN) 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]); -#endif - - if ( - f_compute_rhs_empart(cg->shape, cg->X[0], cg->X[1], cg->X[2], +#endif + + bool used_gpu_substep = false; + if ( + f_compute_rhs_empart(cg->shape, cg->X[0], cg->X[1], cg->X[2], cg->fgfs[phi->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], @@ -1329,11 +1447,20 @@ void bssnEM_class::Step(int lev, int YN) cg->fgfs[Kpsi1->sgfn], cg->fgfs[Kphi1->sgfn], cg->fgfs[rho->sgfn], cg->fgfs[Sx->sgfn], cg->fgfs[Sy->sgfn], cg->fgfs[Sz->sgfn], - cg->fgfs[Sxx->sgfn], cg->fgfs[Sxy->sgfn], cg->fgfs[Sxz->sgfn], - cg->fgfs[Syy->sgfn], cg->fgfs[Syz->sgfn], cg->fgfs[Szz->sgfn], - Symmetry, lev, ndeps) || - f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2], - cg->fgfs[phi->sgfn], cg->fgfs[trK->sgfn], + cg->fgfs[Sxx->sgfn], cg->fgfs[Sxy->sgfn], cg->fgfs[Sxz->sgfn], + cg->fgfs[Syy->sgfn], cg->fgfs[Syz->sgfn], cg->fgfs[Szz->sgfn], + Symmetry, lev, ndeps) || +#if USE_CUDA_BSSN + ((used_gpu_substep = + (run_bssn_em_cuda_substep(cg, SynchList_pre, SynchList_cor, Pp->data, + dT_lev, TRK4, iter_count, Symmetry, lev, + ndeps, cor, chitiny, + rho, Sx, Sy, Sz, Sxx, Sxy, Sxz, Syy, Syz, Szz) == 0)) + ? 0 + : 1) || +#endif + (!used_gpu_substep && f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2], + cg->fgfs[phi->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], @@ -1361,10 +1488,10 @@ void bssnEM_class::Step(int lev, int YN) cg->fgfs[Gamzyy->sgfn], cg->fgfs[Gamzyz->sgfn], cg->fgfs[Gamzzz->sgfn], cg->fgfs[Rxx->sgfn], cg->fgfs[Rxy->sgfn], cg->fgfs[Rxz->sgfn], cg->fgfs[Ryy->sgfn], cg->fgfs[Ryz->sgfn], cg->fgfs[Rzz->sgfn], - cg->fgfs[Cons_Ham->sgfn], - cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn], - cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn], - Symmetry, lev, ndeps, cor)) + cg->fgfs[Cons_Ham->sgfn], + cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn], + cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn], + Symmetry, lev, ndeps, cor))) { cout << "find NaN in domain: (" << cg->bbox[0] << ":" << cg->bbox[3] << "," @@ -1373,11 +1500,15 @@ void bssnEM_class::Step(int lev, int YN) ERROR = 1; } // rk4 substep and boundary - { - MyList *varl0 = StateList, *varl = SynchList_pre, *varl1 = SynchList_cor, *varlrhs = RHSList; - // we do not check the correspondence here - - while (varl0) + { + MyList *varl0 = StateList, *varl = SynchList_pre, *varl1 = SynchList_cor, *varlrhs = RHSList; + // we do not check the correspondence here +#if USE_CUDA_BSSN + if (used_gpu_substep) + skip_bssn_cuda_prefix(varl0, varl, varl1, varlrhs); +#endif + + while (varl0) { #ifndef WithShell if (lev == 0) // sommerfeld indeed diff --git a/AMSS_NCKU_source/bssnEScalar_class.C b/AMSS_NCKU_source/bssnEScalar_class.C index 4d05827..927787a 100644 --- a/AMSS_NCKU_source/bssnEScalar_class.C +++ b/AMSS_NCKU_source/bssnEScalar_class.C @@ -15,10 +15,13 @@ using namespace std; #include "misc.h" #include "Ansorg.h" #include "fmisc.h" -#include "Parallel.h" -#include "bssnEScalar_class.h" -#include "bssn_rhs.h" -#include "initial_puncture.h" +#include "Parallel.h" +#include "bssnEScalar_class.h" +#include "bssn_rhs.h" +#if USE_CUDA_BSSN +#include "bssn_rhs_cuda.h" +#endif +#include "initial_puncture.h" #include "enforce_algebra.h" #include "rungekutta4_rout.h" #include "sommerfeld_rout.h" @@ -26,14 +29,216 @@ using namespace std; #include "shellfunctions.h" #include "parameters.h" -#ifdef With_AHF -#include "derivatives.h" -#include "myglobal.h" -#endif - -//================================================================================================ - -// Define bssnEScalar_class +#ifdef With_AHF +#include "derivatives.h" +#include "myglobal.h" +#endif + +//================================================================================================ + +namespace { + +int amss_escalar_analysis_map_every() +{ + static int every = -1; + if (every < 0) + { + const char *env = getenv("AMSS_ANALYSIS_MAP_EVERY"); + every = (env && atoi(env) > 0) ? atoi(env) : 1; + } + return every; +} + +} + +//================================================================================================ + +#if USE_CUDA_BSSN +extern "C" { +#ifdef fortran1 +void set_escalar_parameter(double &, double &, double &, double &, double &); +#endif +#ifdef fortran2 +void SET_ESCALAR_PARAMETER(double &, double &, double &, double &, double &); +#endif +#ifdef fortran3 +void set_escalar_parameter_(double &, double &, double &, double &, double &); +#endif +} + +namespace { + +bool fill_bssn_cuda_views_prefix(Block *cg, MyList *vars, + double **host_views, + double *propspeeds = nullptr, + double *soa_flat = nullptr) +{ + int idx = 0; + while (vars && idx < BSSN_CUDA_STATE_COUNT) + { + host_views[idx] = cg->fgfs[vars->data->sgfn]; + if (propspeeds) + propspeeds[idx] = vars->data->propspeed; + if (soa_flat) + { + soa_flat[3 * idx + 0] = vars->data->SoA[0]; + soa_flat[3 * idx + 1] = vars->data->SoA[1]; + soa_flat[3 * idx + 2] = vars->data->SoA[2]; + } + vars = vars->next; + ++idx; + } + return idx == BSSN_CUDA_STATE_COUNT; +} + +void skip_bssn_cuda_prefix(MyList *&a, MyList *&b, MyList *&c) +{ + for (int i = 0; i < BSSN_CUDA_STATE_COUNT && a && b && c; ++i) + { + a = a->next; + b = b->next; + c = c->next; + } +} + +void skip_bssn_cuda_prefix(MyList *&a, MyList *&b, + MyList *&c, MyList *&d) +{ + for (int i = 0; i < BSSN_CUDA_STATE_COUNT && a && b && c && d; ++i) + { + a = a->next; + b = b->next; + c = c->next; + d = d->next; + } +} + +MyList *clone_var_list_prefix(MyList *src, int count) +{ + MyList *dst = nullptr; + MyList *tail = nullptr; + for (int i = 0; i < count && src; ++i, src = src->next) + { + MyList *node = new MyList(src->data); + if (!dst) + dst = node; + else + tail->next = node; + tail = node; + } + return dst; +} + +void clear_var_list(MyList *&list) +{ + if (list) + { + list->clearList(); + list = nullptr; + } +} + +void download_bssn_cuda_prefix_if_present(MyList *PatL, + MyList *vars, + int myrank) +{ + while (PatL) + { + MyList *BP = PatL->data->blb; + while (BP) + { + Block *cg = BP->data; + if (myrank == cg->rank) + { + double *views[BSSN_CUDA_STATE_COUNT]; + if (fill_bssn_cuda_views_prefix(cg, vars, views)) + bssn_cuda_download_resident_state_if_present(cg, cg->shape, views); + } + if (BP == PatL->data->ble) + break; + BP = BP->next; + } + PatL = PatL->next; + } +} + +int run_bssn_escalar_cuda_substep(Block *cg, + MyList *state_in_list, + MyList *state_out_list, + Patch *patch, + double &dT_lev, + double &TRK4, + int &iter_count, + int &Symmetry, + int lev, + double &ndeps, + int &co, + double &chitiny, + var *Sphi_in, var *Spi_in, + var *Sphi_rhs, var *Spi_rhs, + var *rho, var *Sx, var *Sy, var *Sz, + var *Sxx, var *Sxy, var *Sxz, + var *Syy, var *Syz, var *Szz) +{ + double *state_in[BSSN_CUDA_STATE_COUNT]; + double *state_out[BSSN_CUDA_STATE_COUNT]; + double propspeed[BSSN_CUDA_STATE_COUNT]; + double soa_flat[3 * BSSN_CUDA_STATE_COUNT]; + if (!fill_bssn_cuda_views_prefix(cg, state_in_list, state_in, propspeed, soa_flat) || + !fill_bssn_cuda_views_prefix(cg, state_out_list, state_out)) + return 1; + + double a2 = 0.0, phi0 = 0.0, r0 = 0.0, sigma0 = 0.0, l2 = 0.0; +#ifdef fortran1 + set_escalar_parameter(a2, phi0, r0, sigma0, l2); +#endif +#ifdef fortran2 + SET_ESCALAR_PARAMETER(a2, phi0, r0, sigma0, l2); +#endif +#ifdef fortran3 + set_escalar_parameter_(a2, phi0, r0, sigma0, l2); +#endif + int apply_enforce_ga = 0; +#if (AGM == 0) + apply_enforce_ga = 1; +#elif (AGM == 1) + apply_enforce_ga = (iter_count == 3) ? 1 : 0; +#endif + if (bssn_cuda_compute_escalar_matter(cg, + cg->shape, cg->X[0], cg->X[1], cg->X[2], + state_in, + cg->fgfs[Sphi_in->sgfn], + cg->fgfs[Spi_in->sgfn], + cg->fgfs[Sphi_rhs->sgfn], + cg->fgfs[Spi_rhs->sgfn], + a2, Symmetry, lev, ndeps, co, apply_enforce_ga)) + return 1; + + int apply_bam_bc = 0; +#if (SommerType == 0) +#ifndef WithShell + apply_bam_bc = (lev == 0) ? 1 : 0; +#endif +#endif + int use_zero_matter = 0; + int keep_resident_state = 1; + double **matter_precomputed = nullptr; + return bssn_cuda_rk4_substep(cg, + cg->shape, cg->X[0], cg->X[1], cg->X[2], + state_in, state_out, matter_precomputed, + propspeed, soa_flat, patch->bbox, + dT_lev, TRK4, iter_count, apply_bam_bc, + Symmetry, lev, ndeps, co, + use_zero_matter, + keep_resident_state, apply_enforce_ga, chitiny); +} + +} +#endif + +//================================================================================================ + +// Define bssnEScalar_class // It inherits some members and methods from the parent class bssn_class and modifies others. // The modified members and methods are defined below (and in the header bssnEScalar_class.h). @@ -41,19 +246,27 @@ using namespace std; //================================================================================================ -bssnEScalar_class::bssnEScalar_class(double Couranti, double StartTimei, double TotalTimei, +bssnEScalar_class::bssnEScalar_class(double Couranti, double StartTimei, double TotalTimei, double DumpTimei, double d2DumpTimei, double CheckTimei, double AnasTimei, - int Symmetryi, int checkruni, char *checkfilenamei, + int Symmetryi, int checkruni, char *checkfilenamei, double numepssi, double numepsbi, double numepshi, - int a_levi, int maxli, int decni, double maxrexi, double drexi) - : bssn_class(Couranti, StartTimei, TotalTimei, - DumpTimei, d2DumpTimei, CheckTimei, AnasTimei, - Symmetryi, checkruni, checkfilenamei, numepssi, numepsbi, numepshi, - a_levi, maxli, decni, maxrexi, drexi) -{ - // setup Monitors - { + int a_levi, int maxli, int decni, double maxrexi, double drexi) + : bssn_class(Couranti, StartTimei, TotalTimei, + DumpTimei, d2DumpTimei, CheckTimei, AnasTimei, + Symmetryi, checkruni, checkfilenamei, numepssi, numepsbi, numepshi, + a_levi, maxli, decni, maxrexi, drexi) +{ + BSSNStateList = nullptr; + BSSNSynchList_pre = nullptr; + BSSNSynchList_cor = nullptr; + ScalarSynchList_pre = nullptr; + ScalarSynchList_cor = nullptr; + sync_cache_scalar_pre = nullptr; + sync_cache_scalar_cor = nullptr; + + // setup Monitors + { char str[50]; stringstream a_stream; a_stream.setf(ios::left); @@ -106,12 +319,22 @@ void bssnEScalar_class::Initialize() ConstraintList->insert(Cons_Gz); - DumpList->insert(Sphi0); - DumpList->insert(Spi0); - DumpList->insert(Cons_fR); - - CheckPoint->addvariablelist(StateList); - CheckPoint->addvariablelist(OldStateList); + DumpList->insert(Sphi0); + DumpList->insert(Spi0); + DumpList->insert(Cons_fR); + +#if USE_CUDA_BSSN + BSSNStateList = clone_var_list_prefix(StateList, BSSN_CUDA_STATE_COUNT); + BSSNSynchList_pre = clone_var_list_prefix(SynchList_pre, BSSN_CUDA_STATE_COUNT); + BSSNSynchList_cor = clone_var_list_prefix(SynchList_cor, BSSN_CUDA_STATE_COUNT); + ScalarSynchList_pre = new MyList(Sphi); + ScalarSynchList_pre->insert(Spi); + ScalarSynchList_cor = new MyList(Sphi1); + ScalarSynchList_cor->insert(Spi1); +#endif + + CheckPoint->addvariablelist(StateList); + CheckPoint->addvariablelist(OldStateList); int myrank = 0; @@ -152,6 +375,12 @@ void bssnEScalar_class::Initialize() #endif Initialize_Level_Runtime(); +#if USE_CUDA_BSSN + if (!sync_cache_scalar_pre) + sync_cache_scalar_pre = new Parallel::SyncCache[GH->levels]; + if (!sync_cache_scalar_cor) + sync_cache_scalar_cor = new Parallel::SyncCache[GH->levels]; +#endif double h = GH->PatL[0]->data->blb->data->getdX(0); for (int i = 1; i < dim; i++) @@ -179,10 +408,34 @@ void bssnEScalar_class::Initialize() //================================================================================================ -bssnEScalar_class::~bssnEScalar_class() -{ - delete Sphio; - delete Spio; +bssnEScalar_class::~bssnEScalar_class() +{ +#if USE_CUDA_BSSN + clear_var_list(BSSNStateList); + clear_var_list(BSSNSynchList_pre); + clear_var_list(BSSNSynchList_cor); + clear_var_list(ScalarSynchList_pre); + clear_var_list(ScalarSynchList_cor); + if (sync_cache_scalar_pre) + { + const int levels = GH ? GH->levels : 0; + for (int i = 0; i < levels; ++i) + sync_cache_scalar_pre[i].destroy(); + delete[] sync_cache_scalar_pre; + sync_cache_scalar_pre = nullptr; + } + if (sync_cache_scalar_cor) + { + const int levels = GH ? GH->levels : 0; + for (int i = 0; i < levels; ++i) + sync_cache_scalar_cor[i].destroy(); + delete[] sync_cache_scalar_cor; + sync_cache_scalar_cor = nullptr; + } +#endif + + delete Sphio; + delete Spio; delete Sphi0; delete Spi0; delete Sphi; @@ -729,20 +982,34 @@ void bssnEScalar_class::Step(int lev, int YN) { MyList *BP = Pp->data->blb; while (BP) - { - Block *cg = BP->data; - if (myrank == cg->rank) - { -#if (AGM == 0) - f_enforce_ga(cg->shape, - cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], + { + Block *cg = BP->data; + if (myrank == cg->rank) + { +#if !USE_CUDA_BSSN +#if (AGM == 0) + f_enforce_ga(cg->shape, + cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn], - cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], - cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn]); -#endif - - if (f_compute_rhs_bssn_escalar(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2], - cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn], + cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], + cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn]); +#endif +#endif + + bool used_gpu_substep = false; + if ( +#if USE_CUDA_BSSN + ((used_gpu_substep = + (run_bssn_escalar_cuda_substep(cg, StateList, SynchList_pre, Pp->data, + dT_lev, TRK4, iter_count, Symmetry, lev, + ndeps, pre, chitiny, + Sphi0, Spi0, Sphi_rhs, Spi_rhs, + rho, Sx, Sy, Sz, Sxx, Sxy, Sxz, Syy, Syz, Szz) == 0)) + ? 0 + : 1) || +#endif + (!used_gpu_substep && f_compute_rhs_bssn_escalar(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2], + cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn], cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn], cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], @@ -773,10 +1040,10 @@ void bssnEScalar_class::Step(int lev, int YN) cg->fgfs[Gamzyy->sgfn], cg->fgfs[Gamzyz->sgfn], cg->fgfs[Gamzzz->sgfn], cg->fgfs[Rxx->sgfn], cg->fgfs[Rxy->sgfn], cg->fgfs[Rxz->sgfn], cg->fgfs[Ryy->sgfn], cg->fgfs[Ryz->sgfn], cg->fgfs[Rzz->sgfn], - cg->fgfs[Cons_Ham->sgfn], - cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn], - cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn], - Symmetry, lev, ndeps, pre)) + cg->fgfs[Cons_Ham->sgfn], + cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn], + cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn], + Symmetry, lev, ndeps, pre))) { cout << "find NaN in domain: (" << cg->bbox[0] << ":" << cg->bbox[3] << "," @@ -786,9 +1053,13 @@ void bssnEScalar_class::Step(int lev, int YN) } // rk4 substep and boundary - { - MyList *varl0 = StateList, *varl = SynchList_pre, *varlrhs = RHSList; // we do not check the correspondence here - while (varl0) + { + MyList *varl0 = StateList, *varl = SynchList_pre, *varlrhs = RHSList; // we do not check the correspondence here +#if USE_CUDA_BSSN + if (used_gpu_substep) + skip_bssn_cuda_prefix(varl0, varl, varlrhs); +#endif + while (varl0) { #ifndef WithShell if (lev == 0) // sommerfeld indeed @@ -823,8 +1094,9 @@ void bssnEScalar_class::Step(int lev, int YN) varlrhs = varlrhs->next; } } - f_lowerboundset(cg->shape, cg->fgfs[phi->sgfn], chitiny); - } + if (!used_gpu_substep) + f_lowerboundset(cg->shape, cg->fgfs[phi->sgfn], chitiny); + } if (BP == Pp->data->ble) break; BP = BP->next; @@ -995,7 +1267,12 @@ void bssnEScalar_class::Step(int lev, int YN) } #endif +#if USE_CUDA_BSSN + Parallel::Sync_cached(GH->PatL[lev], BSSNSynchList_pre, Symmetry, sync_cache_pre[lev]); + Parallel::Sync_cached(GH->PatL[lev], ScalarSynchList_pre, Symmetry, sync_cache_scalar_pre[lev]); +#else Parallel::Sync_cached(GH->PatL[lev], SynchList_pre, Symmetry, sync_cache_pre[lev]); +#endif #ifdef WithShell if (lev == 0) @@ -1065,26 +1342,40 @@ void bssnEScalar_class::Step(int lev, int YN) MyList *BP = Pp->data->blb; while (BP) { - Block *cg = BP->data; - if (myrank == cg->rank) - { -#if (AGM == 0) - f_enforce_ga(cg->shape, - cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], + Block *cg = BP->data; + if (myrank == cg->rank) + { +#if !USE_CUDA_BSSN +#if (AGM == 0) + f_enforce_ga(cg->shape, + 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[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]); #elif (AGM == 1) if (iter_count == 3) f_enforce_ga(cg->shape, - cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->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]); -#endif - - if (f_compute_rhs_bssn_escalar(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2], - cg->fgfs[phi->sgfn], cg->fgfs[trK->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]); +#endif +#endif + + bool used_gpu_substep = false; + if ( +#if USE_CUDA_BSSN + ((used_gpu_substep = + (run_bssn_escalar_cuda_substep(cg, SynchList_pre, SynchList_cor, Pp->data, + dT_lev, TRK4, iter_count, Symmetry, lev, + ndeps, cor, chitiny, + Sphi, Spi, Sphi_rhs, Spi_rhs, + rho, Sx, Sy, Sz, Sxx, Sxy, Sxz, Syy, Syz, Szz) == 0)) + ? 0 + : 1) || +#endif + (!used_gpu_substep && f_compute_rhs_bssn_escalar(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2], + cg->fgfs[phi->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], @@ -1116,10 +1407,10 @@ void bssnEScalar_class::Step(int lev, int YN) cg->fgfs[Gamzyy->sgfn], cg->fgfs[Gamzyz->sgfn], cg->fgfs[Gamzzz->sgfn], cg->fgfs[Rxx->sgfn], cg->fgfs[Rxy->sgfn], cg->fgfs[Rxz->sgfn], cg->fgfs[Ryy->sgfn], cg->fgfs[Ryz->sgfn], cg->fgfs[Rzz->sgfn], - cg->fgfs[Cons_Ham->sgfn], - cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn], - cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn], - Symmetry, lev, ndeps, cor)) + cg->fgfs[Cons_Ham->sgfn], + cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn], + cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn], + Symmetry, lev, ndeps, cor))) { cout << "find NaN in domain: (" << cg->bbox[0] << ":" << cg->bbox[3] << "," @@ -1128,11 +1419,15 @@ void bssnEScalar_class::Step(int lev, int YN) ERROR = 1; } // rk4 substep and boundary - { - MyList *varl0 = StateList, *varl = SynchList_pre, *varl1 = SynchList_cor, *varlrhs = RHSList; - // we do not check the correspondence here - - while (varl0) + { + MyList *varl0 = StateList, *varl = SynchList_pre, *varl1 = SynchList_cor, *varlrhs = RHSList; + // we do not check the correspondence here +#if USE_CUDA_BSSN + if (used_gpu_substep) + skip_bssn_cuda_prefix(varl0, varl, varl1, varlrhs); +#endif + + while (varl0) { #ifndef WithShell if (lev == 0) // sommerfeld indeed @@ -1168,8 +1463,9 @@ void bssnEScalar_class::Step(int lev, int YN) varlrhs = varlrhs->next; } } - f_lowerboundset(cg->shape, cg->fgfs[phi1->sgfn], chitiny); - } + if (!used_gpu_substep) + f_lowerboundset(cg->shape, cg->fgfs[phi1->sgfn], chitiny); + } if (BP == Pp->data->ble) break; BP = BP->next; @@ -1351,7 +1647,12 @@ void bssnEScalar_class::Step(int lev, int YN) } #endif +#if USE_CUDA_BSSN + Parallel::Sync_cached(GH->PatL[lev], BSSNSynchList_cor, Symmetry, sync_cache_cor[lev]); + Parallel::Sync_cached(GH->PatL[lev], ScalarSynchList_cor, Symmetry, sync_cache_scalar_cor[lev]); +#else Parallel::Sync_cached(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_cor[lev]); +#endif #ifdef WithShell if (lev == 0) @@ -1451,9 +1752,13 @@ void bssnEScalar_class::Step(int lev, int YN) } } -#if (RPS == 0) - // mesh refinement boundary part - RestrictProlong(lev, YN, BB); +#if (RPS == 0) + // mesh refinement boundary part +#if USE_CUDA_BSSN + if (!getenv("AMSS_ESCALAR_SPLIT_RP") || atoi(getenv("AMSS_ESCALAR_SPLIT_RP")) == 0) + download_bssn_cuda_prefix_if_present(GH->PatL[lev], SynchList_cor, myrank); +#endif + RestrictProlong(lev, YN, BB); #ifdef WithShell if (lev == 0) @@ -1472,15 +1777,15 @@ void bssnEScalar_class::Step(int lev, int YN) } #endif -#endif - // note the data structure before update +#endif + // note the data structure before update // SynchList_cor 1 ----------- // // StateList 0 ----------- // // OldStateList old ----------- - // update - Pp = GH->PatL[lev]; + // update + Pp = GH->PatL[lev]; while (Pp) { MyList *BP = Pp->data->blb; @@ -2053,17 +2358,34 @@ void bssnEScalar_class::Interp_Constraint() //================================================================================================ -void bssnEScalar_class::Constraint_Out() -{ - // Use the same variables as in the parent class here - // Otherwise the correct time will not be passed - LastConsOut += dT * pow(0.5, Mymax(0, trfls)); - - if (LastConsOut >= AnasTime) - // Constraint violation - { - // recompute least the constraint data lost for moved new grid - for (int lev = 0; lev < GH->levels; lev++) +void bssnEScalar_class::Constraint_Out() +{ + // Use the same variables as in the parent class here + // Otherwise the correct time will not be passed + LastConsOut += dT * pow(0.5, Mymax(0, trfls)); + + if (LastConsOut >= AnasTime) + // Constraint violation + { + const int constraint_map_every = amss_escalar_analysis_map_every(); + static long long constraint_map_counter = 0; + const bool refresh_constraints = + constraint_map_every <= 1 || + (constraint_map_counter % constraint_map_every) == 0; + constraint_map_counter++; + if (!refresh_constraints) + { + LastConsOut = 0; + return; + } + +#if USE_CUDA_BSSN + for (int lev = 0; lev < GH->levels; lev++) + download_bssn_cuda_prefix_if_present(GH->PatL[lev], StateList, myrank); +#endif + + // recompute least the constraint data lost for moved new grid + for (int lev = 0; lev < GH->levels; lev++) { // make sure the data consistent for higher levels { diff --git a/AMSS_NCKU_source/bssnEScalar_class.h b/AMSS_NCKU_source/bssnEScalar_class.h index 3e26005..2e4fac5 100644 --- a/AMSS_NCKU_source/bssnEScalar_class.h +++ b/AMSS_NCKU_source/bssnEScalar_class.h @@ -54,17 +54,21 @@ public: void Interp_Constraint(); void Constraint_Out(); -protected: - var *Sphio, *Spio; - var *Sphi0, *Spi0; +protected: + var *Sphio, *Spio; + var *Sphi0, *Spi0; var *Sphi, *Spi; var *Sphi1, *Spi1; var *Sphi_rhs, *Spi_rhs; - - var *Cons_fR; - - monitor *MaxScalar_Monitor; -}; + + var *Cons_fR; + + MyList *BSSNStateList, *BSSNSynchList_pre, *BSSNSynchList_cor; + MyList *ScalarSynchList_pre, *ScalarSynchList_cor; + Parallel::SyncCache *sync_cache_scalar_pre, *sync_cache_scalar_cor; + + monitor *MaxScalar_Monitor; +}; #endif /* BSSNESCALAR_CLASS_H */ diff --git a/AMSS_NCKU_source/bssnEScalar_rhs.f90 b/AMSS_NCKU_source/bssnEScalar_rhs.f90 index 79327c7..7a60ae7 100644 --- a/AMSS_NCKU_source/bssnEScalar_rhs.f90 +++ b/AMSS_NCKU_source/bssnEScalar_rhs.f90 @@ -3,11 +3,143 @@ !! note that the potential for scalar field in F(R) gravity !! is defined in the file Set_Rho_ADM.f90 -#include "macrodef.fh" - -! rhs for scalar and GR variables -! here we consider vacuum spacetime only - function compute_rhs_bssn_escalar(ex, T,X, Y, Z, & +#include "macrodef.fh" + +! scalar RHS and stress-energy only; BSSN RHS can be supplied by CUDA. + function compute_rhs_bssn_escalar_matter(ex, T, X, Y, Z, & + chi , trK , & + dxx , gxy , gxz , dyy , gyz , dzz, & + Axx , Axy , Axz , Ayy , Ayz , Azz, & + Gamx , Gamy , Gamz , & + Lap , betax , betay , betaz , & + dtSfx , dtSfy , dtSfz , & + Sphi , Spi , & + Sphi_rhs , Spi_rhs , & + rho,Sx,Sy,Sz,Sxx,Sxy,Sxz,Syy,Syz,Szz, & + Symmetry,Lev,eps) result(gont) + implicit none + + integer,intent(in ):: ex(1:3), Symmetry,Lev + real*8, intent(in ):: T + real*8, intent(in ):: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3)) + real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: chi,dxx,dyy,dzz + real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: trK + real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gxy,gxz,gyz + real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Axx,Axy,Axz,Ayy,Ayz,Azz + real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Gamx,Gamy,Gamz + real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Lap, betax, betay, betaz + real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: dtSfx, dtSfy, dtSfz + real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Sphi,Spi + real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Sphi_rhs,Spi_rhs + real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: rho,Sx,Sy,Sz + real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Sxx,Sxy,Sxz,Syy,Syz,Szz + real*8,intent(in) :: eps + integer::gont + + real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz + real*8, dimension(ex(1),ex(2),ex(3)) :: chix,chiy,chiz + real*8, dimension(ex(1),ex(2),ex(3)) :: Lapx,Lapy,Lapz + real*8, dimension(ex(1),ex(2),ex(3)) :: Kx,Ky,Kz,S + real*8, dimension(ex(1),ex(2),ex(3)) :: f,fxx,fxy,fxz,fyy,fyz,fzz + real*8, dimension(ex(1),ex(2),ex(3)) :: alpn1,chin1 + real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz + real*8, dimension(ex(1),ex(2),ex(3)) :: gupyy,gupyz,gupzz + + real*8 :: dX + real*8, parameter :: ZEO=0.d0, ONE = 1.D0, TWO = 2.D0, HALF = 0.5D0 + real*8, parameter :: SYM = 1.D0 + + dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) & + +sum(Gamx)+sum(Gamy)+sum(Gamz) & + +sum(Lap)+sum(Sphi)+sum(Spi) + if(dX.ne.dX) then + if(sum(chi).ne.sum(chi))write(*,*)"bssn_escalar_matter: find NaN in chi" + if(sum(trK).ne.sum(trK))write(*,*)"bssn_escalar_matter: find NaN in trk" + if(sum(dxx).ne.sum(dxx))write(*,*)"bssn_escalar_matter: find NaN in dxx" + if(sum(gxy).ne.sum(gxy))write(*,*)"bssn_escalar_matter: find NaN in gxy" + if(sum(gxz).ne.sum(gxz))write(*,*)"bssn_escalar_matter: find NaN in gxz" + if(sum(dyy).ne.sum(dyy))write(*,*)"bssn_escalar_matter: find NaN in dyy" + if(sum(gyz).ne.sum(gyz))write(*,*)"bssn_escalar_matter: find NaN in gyz" + if(sum(dzz).ne.sum(dzz))write(*,*)"bssn_escalar_matter: find NaN in dzz" + if(sum(Gamx).ne.sum(Gamx))write(*,*)"bssn_escalar_matter: find NaN in Gamx" + if(sum(Gamy).ne.sum(Gamy))write(*,*)"bssn_escalar_matter: find NaN in Gamy" + if(sum(Gamz).ne.sum(Gamz))write(*,*)"bssn_escalar_matter: find NaN in Gamz" + if(sum(Lap).ne.sum(Lap))write(*,*)"bssn_escalar_matter: find NaN in Lap" + if(sum(Sphi).ne.sum(Sphi))write(*,*)"bssn_escalar_matter: find NaN in Sphi" + if(sum(Spi).ne.sum(Spi))write(*,*)"bssn_escalar_matter: find NaN in Spi" + gont = 1 + return + endif + + alpn1 = Lap + ONE + chin1 = chi + ONE + gxx = dxx + ONE + gyy = dyy + ONE + gzz = dzz + ONE + + call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) + call fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) + + gupzz = gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - & + gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz + gupxx = ( gyy * gzz - gyz * gyz ) / gupzz + gupxy = - ( gxy * gzz - gyz * gxz ) / gupzz + gupxz = ( gxy * gyz - gyy * gxz ) / gupzz + gupyy = ( gxx * gzz - gxz * gxz ) / gupzz + gupyz = - ( gxx * gyz - gxy * gxz ) / gupzz + gupzz = ( gxx * gyy - gxy * gxy ) / gupzz + +#if 1 + Sphi_rhs = alpn1 * Spi + call fderivs(ex,Sphi,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) + call fdderivs(ex,Sphi,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) + Spi_rhs = gupxx * fxx + gupyy * fyy + gupzz * fzz + & + ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO - & + ((Gamx+(gupxx*chix+gupxy*chiy+gupxz*chiz)/TWO/chin1)*Kx & + + (Gamy+(gupxy*chix+gupyy*chiy+gupyz*chiz)/TWO/chin1)*Ky & + + (Gamz+(gupxz*chix+gupyz*chiy+gupzz*chiz)/TWO/chin1)*Kz) + Spi_rhs = Spi_rhs*alpn1 + & + (gupxx*Lapx*Kx + gupxy*Lapx*Ky + gupxz*Lapx*Kz & + +gupxy*Lapy*Kx + gupyy*Lapy*Ky + gupyz*Lapy*Kz & + +gupxz*Lapz*Kx + gupyz*Lapz*Ky + gupzz*Lapz*Kz) + + call frpotential(ex,Sphi,f,S) + Spi_rhs = Spi_rhs*chin1 + alpn1*(trK*Spi - S) + rho = chin1*((gupxx * Kx * Kx + gupyy * Ky * Ky + gupzz * Kz * Kz)/TWO + & + gupxy * Kx * Ky + gupxz * Kx * Kz + gupyz * Ky * Kz ) & + + Spi*Spi/TWO+f + Sx = -Spi*Kx + Sy = -Spi*Ky + Sz = -Spi*Kz + f = (rho - Spi*Spi)/chin1 + Sxx = Kx*Kx-f*gxx + Sxy = Kx*Ky-f*gxy + Sxz = Kx*Kz-f*gxz + Syy = Ky*Ky-f*gyy + Syz = Ky*Kz-f*gyz + Szz = Kz*Kz-f*gzz +#else + Sphi_rhs = ZEO + Spi_rhs = ZEO + rho = ZEO + Sx = ZEO + Sy = ZEO + Sz = ZEO + Sxx = ZEO + Sxy = ZEO + Sxz = ZEO + Syy = ZEO + Syz = ZEO + Szz = ZEO +#endif + + gont = 0 + return +end function compute_rhs_bssn_escalar_matter + +! rhs for scalar and GR variables +! here we consider vacuum spacetime only + function compute_rhs_bssn_escalar(ex, T,X, Y, Z, & chi , trK , & dxx , gxy , gxz , dyy , gyz , dzz, & Axx , Axy , Axz , Ayy , Ayz , Azz, & diff --git a/AMSS_NCKU_source/bssn_class.C b/AMSS_NCKU_source/bssn_class.C index c4164b3..8b1860d 100644 --- a/AMSS_NCKU_source/bssn_class.C +++ b/AMSS_NCKU_source/bssn_class.C @@ -79,6 +79,111 @@ int amss_analysis_map_every() return every; } +#if USE_CUDA_BSSN +int amss_escalar_split_rp_enabled() +{ + static int enabled = -1; + if (enabled < 0) + { + const char *env = getenv("AMSS_ESCALAR_SPLIT_RP"); + enabled = (env && atoi(env) != 0) ? 1 : 0; + } + return enabled; +} + +int amss_escalar_split_rp_recursive_enabled() +{ + static int enabled = -1; + if (enabled < 0) + { + const char *env = getenv("AMSS_ESCALAR_SPLIT_RP_RECURSIVE"); + enabled = (env && atoi(env) != 0) ? 1 : 0; + } + return enabled; +} + +MyList *clone_var_sublist(MyList *src, int skip, int take) +{ + for (int i = 0; i < skip && src; ++i) + src = src->next; + MyList *dst = nullptr; + MyList *tail = nullptr; + int copied = 0; + while (src && (take < 0 || copied < take)) + { + MyList *node = new MyList(src->data); + if (!dst) + dst = node; + else + tail->next = node; + tail = node; + src = src->next; + ++copied; + } + return dst; +} + +void clear_tmp_var_list(MyList *&list) +{ + if (list) + { + list->clearList(); + list = nullptr; + } +} + +int var_list_count(MyList *vars) +{ + int count = 0; + while (vars) + { + ++count; + vars = vars->next; + } + return count; +} + +bool bssn_prefix_views(Block *cg, MyList *vars, double **views) +{ + if (!cg || !vars || !views) + return false; + for (int i = 0; i < BSSN_CUDA_STATE_COUNT; ++i) + { + if (!vars) + return false; + views[i] = cg->fgfs[vars->data->sgfn]; + if (!views[i]) + return false; + vars = vars->next; + } + return true; +} + +void download_bssn_prefix_for_list(MyList *PatL, + MyList *vars, + int myrank) +{ + while (PatL) + { + MyList *BP = PatL->data->blb; + while (BP) + { + Block *cg = BP->data; + if (myrank == cg->rank) + { + double *views[BSSN_CUDA_STATE_COUNT]; + if (bssn_prefix_views(cg, vars, views)) + bssn_cuda_download_resident_state_if_present(cg, cg->shape, views); + } + if (BP == PatL->data->ble) + break; + BP = BP->next; + } + PatL = PatL->next; + } +} +#endif + } // Compile-time switch for per-timestep memory usage collection/printing. @@ -7000,6 +7105,108 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB, // a_stream.setf(ios::left); #endif +#if USE_CUDA_BSSN && (ABEtype == 1) && (RPB == 0) && (MIXOUTB == 0) + if (lev > 0 && amss_escalar_split_rp_recursive_enabled() && var_list_count(SL) > BSSN_CUDA_STATE_COUNT) + { + MyList *SLb = clone_var_sublist(SL, 0, BSSN_CUDA_STATE_COUNT); + MyList *OLb = clone_var_sublist(OL, 0, BSSN_CUDA_STATE_COUNT); + MyList *corLb = clone_var_sublist(corL, 0, BSSN_CUDA_STATE_COUNT); + MyList *preb = clone_var_sublist(SynchList_pre, 0, BSSN_CUDA_STATE_COUNT); + MyList *SLs = clone_var_sublist(SL, BSSN_CUDA_STATE_COUNT, -1); + MyList *OLs = clone_var_sublist(OL, BSSN_CUDA_STATE_COUNT, -1); + MyList *corLs = clone_var_sublist(corL, BSSN_CUDA_STATE_COUNT, -1); + MyList *pres = clone_var_sublist(SynchList_pre, BSSN_CUDA_STATE_COUNT, -1); + + if (lev > trfls && YN == 0) + { + MyList *Pp = GH->PatL[lev - 1]; + while (Pp) + { + if (BB) + { + Parallel::prepare_inter_time_level(Pp->data, SLb, OLb, corLb, preb, 0); + Parallel::prepare_inter_time_level(Pp->data, SLs, OLs, corLs, pres, 0); + } + else + { + Parallel::prepare_inter_time_level(Pp->data, SLb, OLb, preb, 0); + Parallel::prepare_inter_time_level(Pp->data, SLs, OLs, pres, 0); + } + Pp = Pp->next; + } +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], SLb, preb, Symmetry); +#else + Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], SLb, preb, Symmetry, sync_cache_restrict[lev]); +#endif + Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], SLs, pres, Symmetry); +#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1) +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::Sync(GH->PatL[lev - 1], preb, Symmetry); +#else + Parallel::Sync_cached(GH->PatL[lev - 1], preb, Symmetry, sync_cache_rp_coarse[lev]); +#endif + Parallel::Sync(GH->PatL[lev - 1], pres, Symmetry); +#endif +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], preb, SLb, Symmetry); +#else + Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], preb, SLb, Symmetry, sync_cache_outbd[lev]); +#endif + Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], pres, SLs, Symmetry); + } + else + { +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], SLb, SLb, Symmetry); +#else + Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], SLb, SLb, Symmetry, sync_cache_restrict[lev]); +#endif + Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], SLs, SLs, Symmetry); +#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1) +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::Sync(GH->PatL[lev - 1], SLb, Symmetry); +#else + Parallel::Sync_cached(GH->PatL[lev - 1], SLb, Symmetry, sync_cache_rp_coarse[lev]); +#endif + Parallel::Sync(GH->PatL[lev - 1], SLs, Symmetry); +#endif +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], SLb, SLb, Symmetry); +#else + Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], SLb, SLb, Symmetry, sync_cache_outbd[lev]); +#endif + Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], SLs, SLs, Symmetry); + } + +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::Sync(GH->PatL[lev], SLb, Symmetry); +#else + Parallel::Sync_cached(GH->PatL[lev], SLb, Symmetry, sync_cache_rp_fine[lev]); +#endif + Parallel::Sync(GH->PatL[lev], SLs, Symmetry); + + clear_tmp_var_list(SLb); + clear_tmp_var_list(OLb); + clear_tmp_var_list(corLb); + clear_tmp_var_list(preb); + clear_tmp_var_list(SLs); + clear_tmp_var_list(OLs); + clear_tmp_var_list(corLs); + clear_tmp_var_list(pres); + STEP_TIMER_ADD(TB_RESTRICT_PROLONG, timer_restrict_prolong); + return; + } + if (lev > 0 && var_list_count(SL) > BSSN_CUDA_STATE_COUNT) + { + download_bssn_prefix_for_list(GH->PatL[lev], SL, myrank); + download_bssn_prefix_for_list(GH->PatL[lev - 1], SL, myrank); + download_bssn_prefix_for_list(GH->PatL[lev - 1], OL, myrank); + if (BB) + download_bssn_prefix_for_list(GH->PatL[lev - 1], corL, myrank); + } +#endif + if (lev > 0) { MyList *Pp, *Ppc; @@ -7355,6 +7562,117 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB) // OldStateList 0 ----------- // // SynchList_cor old ----------- +#if USE_CUDA_BSSN && (ABEtype == 1) && (RPB == 0) && (MIXOUTB == 0) + if (lev > 0 && amss_escalar_split_rp_enabled() && + var_list_count(StateList) > BSSN_CUDA_STATE_COUNT) + { + MyList *StateB = clone_var_sublist(StateList, 0, BSSN_CUDA_STATE_COUNT); + MyList *OldB = clone_var_sublist(OldStateList, 0, BSSN_CUDA_STATE_COUNT); + MyList *PreB = clone_var_sublist(SynchList_pre, 0, BSSN_CUDA_STATE_COUNT); + MyList *CorB = clone_var_sublist(SynchList_cor, 0, BSSN_CUDA_STATE_COUNT); + MyList *StateS = clone_var_sublist(StateList, BSSN_CUDA_STATE_COUNT, -1); + MyList *OldS = clone_var_sublist(OldStateList, BSSN_CUDA_STATE_COUNT, -1); + MyList *PreS = clone_var_sublist(SynchList_pre, BSSN_CUDA_STATE_COUNT, -1); + MyList *CorS = clone_var_sublist(SynchList_cor, BSSN_CUDA_STATE_COUNT, -1); + + if (lev > trfls && YN == 0) + { + if (myrank == 0) + cout << "/=: " << GH->Lt[lev - 1] << "," << GH->Lt[lev] + dT_lev << endl; + MyList *Pp = GH->PatL[lev - 1]; + while (Pp) + { + if (BB) + { + Parallel::prepare_inter_time_level(Pp->data, StateB, OldB, CorB, PreB, 0); + Parallel::prepare_inter_time_level(Pp->data, StateS, OldS, CorS, PreS, 0); + } + else + { + Parallel::prepare_inter_time_level(Pp->data, StateB, OldB, PreB, 0); + Parallel::prepare_inter_time_level(Pp->data, StateS, OldS, PreS, 0); + } + Pp = Pp->next; + } + +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], CorB, PreB, Symmetry); + Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], CorS, PreS, Symmetry); +#else + Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], CorB, PreB, Symmetry, sync_cache_restrict[lev]); + Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], CorS, PreS, Symmetry, sync_cache_restrict[lev]); +#endif + +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::Sync(GH->PatL[lev - 1], PreB, Symmetry); + Parallel::Sync(GH->PatL[lev - 1], PreS, Symmetry); +#else +#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1) + Parallel::Sync_cached(GH->PatL[lev - 1], PreB, Symmetry, sync_cache_rp_coarse[lev]); + Parallel::Sync_cached(GH->PatL[lev - 1], PreS, Symmetry, sync_cache_rp_coarse[lev]); +#endif +#endif + +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], PreB, CorB, Symmetry); + Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], PreS, CorS, Symmetry); +#else + Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], PreB, CorB, Symmetry, sync_cache_outbd[lev]); + Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], PreS, CorS, Symmetry, sync_cache_outbd[lev]); +#endif + } + else + { + if (myrank == 0) + cout << "===: " << GH->Lt[lev - 1] << "," << GH->Lt[lev] + dT_lev << endl; +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], CorB, StateB, Symmetry); + Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], CorS, StateS, Symmetry); +#else + Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], CorB, StateB, Symmetry, sync_cache_restrict[lev]); + Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], CorS, StateS, Symmetry, sync_cache_restrict[lev]); +#endif + +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::Sync(GH->PatL[lev - 1], StateB, Symmetry); + Parallel::Sync(GH->PatL[lev - 1], StateS, Symmetry); +#else +#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1) + Parallel::Sync_cached(GH->PatL[lev - 1], StateB, Symmetry, sync_cache_rp_coarse[lev]); + Parallel::Sync_cached(GH->PatL[lev - 1], StateS, Symmetry, sync_cache_rp_coarse[lev]); +#endif +#endif + +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], StateB, CorB, Symmetry); + Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], StateS, CorS, Symmetry); +#else + Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], StateB, CorB, Symmetry, sync_cache_outbd[lev]); + Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], StateS, CorS, Symmetry, sync_cache_outbd[lev]); +#endif + } + +#if AMSS_LEGACY_ABE_TRANSFER + Parallel::Sync(GH->PatL[lev], CorB, Symmetry); + Parallel::Sync(GH->PatL[lev], CorS, Symmetry); +#else + Parallel::Sync_cached(GH->PatL[lev], CorB, Symmetry, sync_cache_rp_fine[lev]); + Parallel::Sync_cached(GH->PatL[lev], CorS, Symmetry, sync_cache_rp_fine[lev]); +#endif + + clear_tmp_var_list(StateB); + clear_tmp_var_list(OldB); + clear_tmp_var_list(PreB); + clear_tmp_var_list(CorB); + clear_tmp_var_list(StateS); + clear_tmp_var_list(OldS); + clear_tmp_var_list(PreS); + clear_tmp_var_list(CorS); + STEP_TIMER_ADD(TB_RESTRICT_PROLONG, timer_restrict_prolong); + return; + } +#endif + if (lev > 0) { MyList *Pp, *Ppc; diff --git a/AMSS_NCKU_source/bssn_rhs.h b/AMSS_NCKU_source/bssn_rhs.h index 96a545a..cc0c520 100644 --- a/AMSS_NCKU_source/bssn_rhs.h +++ b/AMSS_NCKU_source/bssn_rhs.h @@ -5,8 +5,9 @@ #ifdef fortran1 #define f_compute_rhs_bssn compute_rhs_bssn #define f_compute_rhs_bssn_ss compute_rhs_bssn_ss -#define f_compute_rhs_bssn_escalar compute_rhs_bssn_escalar -#define f_compute_rhs_bssn_escalar_ss compute_rhs_bssn_escalar_ss +#define f_compute_rhs_bssn_escalar compute_rhs_bssn_escalar +#define f_compute_rhs_bssn_escalar_matter compute_rhs_bssn_escalar_matter +#define f_compute_rhs_bssn_escalar_ss compute_rhs_bssn_escalar_ss #define f_compute_rhs_Z4c compute_rhs_z4c #define f_compute_rhs_Z4cnot compute_rhs_z4cnot #define f_compute_rhs_Z4c_ss compute_rhs_z4c_ss @@ -15,8 +16,9 @@ #ifdef fortran2 #define f_compute_rhs_bssn COMPUTE_RHS_BSSN #define f_compute_rhs_bssn_ss COMPUTE_RHS_BSSN_SS -#define f_compute_rhs_bssn_escalar COMPUTE_RHS_BSSN_ESCALAR -#define f_compute_rhs_bssn_escalar_ss COMPUTE_RHS_BSSN_ESCALAR_SS +#define f_compute_rhs_bssn_escalar COMPUTE_RHS_BSSN_ESCALAR +#define f_compute_rhs_bssn_escalar_matter COMPUTE_RHS_BSSN_ESCALAR_MATTER +#define f_compute_rhs_bssn_escalar_ss COMPUTE_RHS_BSSN_ESCALAR_SS #define f_compute_rhs_Z4c COMPUTE_RHS_Z4C #define f_compute_rhs_Z4cnot COMPUTE_RHS_Z4CNOT #define f_compute_rhs_Z4c_ss COMPUTE_RHS_Z4C_SS @@ -25,8 +27,9 @@ #ifdef fortran3 #define f_compute_rhs_bssn compute_rhs_bssn_ #define f_compute_rhs_bssn_ss compute_rhs_bssn_ss_ -#define f_compute_rhs_bssn_escalar compute_rhs_bssn_escalar_ -#define f_compute_rhs_bssn_escalar_ss compute_rhs_bssn_escalar_ss_ +#define f_compute_rhs_bssn_escalar compute_rhs_bssn_escalar_ +#define f_compute_rhs_bssn_escalar_matter compute_rhs_bssn_escalar_matter_ +#define f_compute_rhs_bssn_escalar_ss compute_rhs_bssn_escalar_ss_ #define f_compute_rhs_Z4c compute_rhs_z4c_ #define f_compute_rhs_Z4cnot compute_rhs_z4cnot_ #define f_compute_rhs_Z4c_ss compute_rhs_z4c_ss_ @@ -96,10 +99,24 @@ extern "C" int &, int &, double &, int &, int &); } -extern "C" -{ - int f_compute_rhs_bssn_escalar(int *, double &, double *, double *, double *, // ex,T,X,Y,Z - double *, double *, // chi, trK +extern "C" +{ + int f_compute_rhs_bssn_escalar_matter(int *, double &, double *, double *, double *, // ex,T,X,Y,Z + double *, double *, // chi, trK + double *, double *, double *, double *, double *, double *, // gij + double *, double *, double *, double *, double *, double *, // Aij + double *, double *, double *, // Gam + double *, double *, double *, double *, double *, double *, double *, // Gauge + double *, double *, // Sphi, Spi + double *, double *, // Sphi, Spi rhs + double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, // stress-energy + int &, int &, double &); +} + +extern "C" +{ + int f_compute_rhs_bssn_escalar(int *, double &, double *, double *, double *, // ex,T,X,Y,Z + double *, double *, // chi, trK double *, double *, double *, double *, double *, double *, // gij double *, double *, double *, double *, double *, double *, // Aij double *, double *, double *, // Gam diff --git a/AMSS_NCKU_source/bssn_rhs_cuda.cu b/AMSS_NCKU_source/bssn_rhs_cuda.cu index 378c79e..61a793e 100644 --- a/AMSS_NCKU_source/bssn_rhs_cuda.cu +++ b/AMSS_NCKU_source/bssn_rhs_cuda.cu @@ -4947,6 +4947,118 @@ static void bind_state_output_slots(const std::array } } +__global__ __launch_bounds__(128, 4) +void kern_escalar_matter_rhs( + const double * __restrict__ chi, + const double * __restrict__ trK, + const double * __restrict__ dxx, + const double * __restrict__ gxy, + const double * __restrict__ gxz, + const double * __restrict__ dyy, + const double * __restrict__ gyz, + const double * __restrict__ dzz, + const double * __restrict__ Gamx, + const double * __restrict__ Gamy, + const double * __restrict__ Gamz, + const double * __restrict__ Lap, + const double * __restrict__ Sphi, + const double * __restrict__ Spi, + const double * __restrict__ chix, + const double * __restrict__ chiy, + const double * __restrict__ chiz, + const double * __restrict__ Lapx, + const double * __restrict__ Lapy, + const double * __restrict__ Lapz, + const double * __restrict__ Kx, + const double * __restrict__ Ky, + const double * __restrict__ Kz, + const double * __restrict__ fxx, + const double * __restrict__ fxy, + const double * __restrict__ fxz, + const double * __restrict__ fyy, + const double * __restrict__ fyz, + const double * __restrict__ fzz, + double * __restrict__ Sphi_rhs, + double * __restrict__ Spi_rhs, + double * __restrict__ rho, + double * __restrict__ Sx, + double * __restrict__ Sy, + double * __restrict__ Sz, + double * __restrict__ Sxx, + double * __restrict__ Sxy, + double * __restrict__ Sxz, + double * __restrict__ Syy, + double * __restrict__ Syz, + double * __restrict__ Szz, + double a2) +{ + const double TWO = 2.0; + const double HALF = 0.5; + + for (int i = blockIdx.x * blockDim.x + threadIdx.x; + i < d_gp.all; + i += blockDim.x * gridDim.x) + { + const double chin1 = chi[i] + 1.0; + const double alpn1 = Lap[i] + 1.0; + const double gxx = dxx[i] + 1.0; + const double gyy = dyy[i] + 1.0; + const double gzz = dzz[i] + 1.0; + const double lgxy = gxy[i]; + const double lgxz = gxz[i]; + const double lgyz = gyz[i]; + + double det = gxx * gyy * gzz + lgxy * lgyz * lgxz + lgxz * lgxy * lgyz + - lgxz * gyy * lgxz - lgxy * lgxy * gzz - gxx * lgyz * lgyz; + const double gupxx = (gyy * gzz - lgyz * lgyz) / det; + const double gupxy = -(lgxy * gzz - lgyz * lgxz) / det; + const double gupxz = (lgxy * lgyz - gyy * lgxz) / det; + const double gupyy = (gxx * gzz - lgxz * lgxz) / det; + const double gupyz = -(gxx * lgyz - lgxy * lgxz) / det; + const double gupzz = (gxx * gyy - lgxy * lgxy) / det; + + double V = 0.0; + double dVdSphi = 0.0; +#if (EScalar_CC == 2 || EScalar_CC == 3) + const double sqrt_pi_over_3 = sqrt(PI_VAL / 3.0); + const double e4 = exp(4.0 * sqrt_pi_over_3 * Sphi[i]); + const double e8n = exp(-8.0 * sqrt_pi_over_3 * Sphi[i]); + const double inv_a2 = 1.0 / a2; + V = e8n * (1.0 - e4) * (1.0 - e4) / (32.0 * PI_VAL) * inv_a2; + dVdSphi = inv_a2 / 12.0 * sqrt(3.0 / PI_VAL) * e8n * (-1.0 + e4); +#else + (void)a2; +#endif + + Sphi_rhs[i] = alpn1 * Spi[i]; + + double srhs = gupxx * fxx[i] + gupyy * fyy[i] + gupzz * fzz[i] + + TWO * (gupxy * fxy[i] + gupxz * fxz[i] + gupyz * fyz[i]) + - ((Gamx[i] + (gupxx * chix[i] + gupxy * chiy[i] + gupxz * chiz[i]) * HALF / chin1) * Kx[i] + + (Gamy[i] + (gupxy * chix[i] + gupyy * chiy[i] + gupyz * chiz[i]) * HALF / chin1) * Ky[i] + + (Gamz[i] + (gupxz * chix[i] + gupyz * chiy[i] + gupzz * chiz[i]) * HALF / chin1) * Kz[i]); + srhs = srhs * alpn1 + + (gupxx * Lapx[i] * Kx[i] + gupxy * Lapx[i] * Ky[i] + gupxz * Lapx[i] * Kz[i] + + gupxy * Lapy[i] * Kx[i] + gupyy * Lapy[i] * Ky[i] + gupyz * Lapy[i] * Kz[i] + + gupxz * Lapz[i] * Kx[i] + gupyz * Lapz[i] * Ky[i] + gupzz * Lapz[i] * Kz[i]); + Spi_rhs[i] = srhs * chin1 + alpn1 * (trK[i] * Spi[i] - dVdSphi); + + rho[i] = chin1 * ((gupxx * Kx[i] * Kx[i] + gupyy * Ky[i] * Ky[i] + gupzz * Kz[i] * Kz[i]) * HALF + + gupxy * Kx[i] * Ky[i] + gupxz * Kx[i] * Kz[i] + gupyz * Ky[i] * Kz[i]) + + HALF * Spi[i] * Spi[i] + V; + Sx[i] = -Spi[i] * Kx[i]; + Sy[i] = -Spi[i] * Ky[i]; + Sz[i] = -Spi[i] * Kz[i]; + const double f = (rho[i] - Spi[i] * Spi[i]) / chin1; + Sxx[i] = Kx[i] * Kx[i] - f * gxx; + Sxy[i] = Kx[i] * Ky[i] - f * lgxy; + Sxz[i] = Kx[i] * Kz[i] - f * lgxz; + Syy[i] = Ky[i] * Ky[i] - f * gyy; + Syz[i] = Ky[i] * Kz[i] - f * lgyz; + Szz[i] = Kz[i] * Kz[i] - f * gzz; + } +} + static bool resident_key_matches(const StepContext &ctx, int bank, double **host_key) { if (!host_key || bank < 0 || bank >= BSSN_RESIDENT_BANK_COUNT) @@ -6885,6 +6997,97 @@ int f_compute_rhs_bssn(int *ex, double &T, return 0; } +extern "C" +int bssn_cuda_compute_escalar_matter(void *block_tag, + int *ex, double *X, double *Y, double *Z, + double **state_host_in, + double *Sphi_host, + double *Spi_host, + double *Sphi_rhs_host, + double *Spi_rhs_host, + double a2, + int &Symmetry, + int &Lev, + double &eps, + int &co, + int &apply_enforce_ga) +{ + init_gpu_dispatch(); + CUDA_CHECK(cudaSetDevice(g_dispatch.my_device)); + + const size_t all = (size_t)ex[0] * ex[1] * ex[2]; + const size_t bytes = all * sizeof(double); + setup_grid_params(ex, X, Y, Z, Symmetry, eps, co); + StepContext &ctx = ensure_step_ctx(block_tag, all); + + const int input_bank = ensure_resident_bank(ctx, state_host_in, all, true); + mark_resident_current_bank(ctx, input_bank); + bind_state_input_slots(ctx.d_resident[input_bank]); + if (apply_enforce_ga) { + kern_enforce_ga_cuda<<>>(g_buf.slot[S_dxx], g_buf.slot[S_gxy], g_buf.slot[S_gxz], + g_buf.slot[S_dyy], g_buf.slot[S_gyz], g_buf.slot[S_dzz], + g_buf.slot[S_Axx], g_buf.slot[S_Axy], g_buf.slot[S_Axz], + g_buf.slot[S_Ayy], g_buf.slot[S_Ayz], g_buf.slot[S_Azz]); + set_resident_host_clean(ctx, input_bank, false); + } + CUDA_CHECK(cudaMemcpyAsync(g_buf.slot[S_S_arr], Sphi_host, bytes, cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpyAsync(g_buf.slot[S_f_arr], Spi_host, bytes, cudaMemcpyHostToDevice)); + + double *src_fields[3] = { + g_buf.slot[S_chi], g_buf.slot[S_Lap], g_buf.slot[S_S_arr] + }; + double *fx_fields[3] = { + g_buf.slot[S_chix], g_buf.slot[S_Lapx], g_buf.slot[S_Kx] + }; + double *fy_fields[3] = { + g_buf.slot[S_chiy], g_buf.slot[S_Lapy], g_buf.slot[S_Ky] + }; + double *fz_fields[3] = { + g_buf.slot[S_chiz], g_buf.slot[S_Lapz], g_buf.slot[S_Kz] + }; + int soa_signs[9] = { + 1, 1, 1, + 1, 1, 1, + 1, 1, 1 + }; + gpu_fderivs_batch(3, src_fields, fx_fields, fy_fields, fz_fields, + soa_signs, (int)all); + + double *fd_src[1] = { g_buf.slot[S_S_arr] }; + double *fxx_fields[1] = { g_buf.slot[S_fxx] }; + double *fxy_fields[1] = { g_buf.slot[S_fxy] }; + double *fxz_fields[1] = { g_buf.slot[S_fxz] }; + double *fyy_fields[1] = { g_buf.slot[S_fyy] }; + double *fyz_fields[1] = { g_buf.slot[S_fyz] }; + double *fzz_fields[1] = { g_buf.slot[S_fzz] }; + gpu_fdderivs_batch(1, fd_src, fxx_fields, fxy_fields, fxz_fields, + fyy_fields, fyz_fields, fzz_fields, soa_signs, (int)all); + + kern_escalar_matter_rhs<<>>( + g_buf.slot[S_chi], g_buf.slot[S_trK], + g_buf.slot[S_dxx], g_buf.slot[S_gxy], g_buf.slot[S_gxz], + g_buf.slot[S_dyy], g_buf.slot[S_gyz], g_buf.slot[S_dzz], + g_buf.slot[S_Gamx], g_buf.slot[S_Gamy], g_buf.slot[S_Gamz], + g_buf.slot[S_Lap], + g_buf.slot[S_S_arr], g_buf.slot[S_f_arr], + g_buf.slot[S_chix], g_buf.slot[S_chiy], g_buf.slot[S_chiz], + g_buf.slot[S_Lapx], g_buf.slot[S_Lapy], g_buf.slot[S_Lapz], + g_buf.slot[S_Kx], g_buf.slot[S_Ky], g_buf.slot[S_Kz], + g_buf.slot[S_fxx], g_buf.slot[S_fxy], g_buf.slot[S_fxz], + g_buf.slot[S_fyy], g_buf.slot[S_fyz], g_buf.slot[S_fzz], + g_buf.slot[S_Gamxa], g_buf.slot[S_Gamya], + ctx.d_matter[0], ctx.d_matter[1], ctx.d_matter[2], ctx.d_matter[3], + ctx.d_matter[4], ctx.d_matter[5], ctx.d_matter[6], + ctx.d_matter[7], ctx.d_matter[8], ctx.d_matter[9], + a2); + CUDA_CHECK(cudaMemcpyAsync(Sphi_rhs_host, g_buf.slot[S_Gamxa], bytes, cudaMemcpyDeviceToHost)); + CUDA_CHECK(cudaMemcpyAsync(Spi_rhs_host, g_buf.slot[S_Gamya], bytes, cudaMemcpyDeviceToHost)); + CUDA_CHECK(cudaDeviceSynchronize()); + ctx.matter_ready = true; + (void)Lev; + return 0; +} + extern "C" int bssn_cuda_rk4_substep(void *block_tag, int *ex, double *X, double *Y, double *Z, @@ -6968,6 +7171,8 @@ int bssn_cuda_rk4_substep(void *block_tag, if (RK4 == 0) { if (use_zero_matter) { if (!ctx.matter_ready) zero_matter_cache(ctx, all); + } else if (!matter_host) { + if (!ctx.matter_ready) return 1; } else { upload_matter_cache(ctx, matter_host, all); } @@ -6979,7 +7184,8 @@ int bssn_cuda_rk4_substep(void *block_tag, cudaMemcpyDeviceToDevice)); } else if (!ctx.matter_ready) { if (use_zero_matter) zero_matter_cache(ctx, all); - else upload_matter_cache(ctx, matter_host, all); + else if (matter_host) upload_matter_cache(ctx, matter_host, all); + else return 1; } bind_matter_slots(ctx); if (profile) { @@ -7107,6 +7313,20 @@ int bssn_cuda_download_resident_state_if_present(void *block_tag, return 0; } +extern "C" +int bssn_cuda_resident_state_matches(void *block_tag, + double **state_host_key) +{ + init_gpu_dispatch(); + CUDA_CHECK(cudaSetDevice(g_dispatch.my_device)); + auto it = g_step_ctx.find(block_tag); + if (it == g_step_ctx.end() || !resident_key_usable(state_host_key)) + return 0; + StepContext &ctx = it->second; + const int bank = find_resident_bank(ctx, state_host_key); + return (bank >= 0 && ctx.resident_valid[bank]) ? 1 : 0; +} + extern "C" int bssn_cuda_download_constraint_outputs(int *ex, double **constraint_host_out) diff --git a/AMSS_NCKU_source/bssn_rhs_cuda.h b/AMSS_NCKU_source/bssn_rhs_cuda.h index 41e92f8..63f5f59 100644 --- a/AMSS_NCKU_source/bssn_rhs_cuda.h +++ b/AMSS_NCKU_source/bssn_rhs_cuda.h @@ -55,6 +55,20 @@ int bssn_cuda_rk4_substep(void *block_tag, int &apply_enforce_ga, double &chitiny); +int bssn_cuda_compute_escalar_matter(void *block_tag, + int *ex, double *X, double *Y, double *Z, + double **state_host_in, + double *Sphi_host, + double *Spi_host, + double *Sphi_rhs_host, + double *Spi_rhs_host, + double a2, + int &Symmetry, + int &Lev, + double &eps, + int &co, + int &apply_enforce_ga); + int bssn_cuda_copy_state_region_to_host(void *block_tag, int state_index, double *host_state, @@ -77,6 +91,9 @@ int bssn_cuda_download_resident_state_if_present(void *block_tag, int *ex, double **state_host_out); +int bssn_cuda_resident_state_matches(void *block_tag, + double **state_host_key); + int bssn_cuda_download_constraint_outputs(int *ex, double **constraint_host_out); diff --git a/AMSS_NCKU_source/macrodef.fh b/AMSS_NCKU_source/macrodef.fh index 84fbf50..77036c3 100644 --- a/AMSS_NCKU_source/macrodef.fh +++ b/AMSS_NCKU_source/macrodef.fh @@ -13,7 +13,7 @@ #define ABV 0 -#define EScalar_CC 2 +#define EScalar_CC 2 #if 0 diff --git a/AMSS_NCKU_source/macrodef.h b/AMSS_NCKU_source/macrodef.h index 7160d1c..3b3614e 100644 --- a/AMSS_NCKU_source/macrodef.h +++ b/AMSS_NCKU_source/macrodef.h @@ -10,7 +10,7 @@ #define GaussInt -#define ABEtype 2 +#define ABEtype 1 //#define With_AHF #define Psi4type 0 diff --git a/makefile_and_run.py b/makefile_and_run.py index 7a86028..c3f5dd6 100755 --- a/makefile_and_run.py +++ b/makefile_and_run.py @@ -152,6 +152,9 @@ def _gpu_runtime_env(): "AMSS_CUDA_AMR_RESTRICT_BATCH": "0", "AMSS_CUDA_DEVICE_SEGMENT_BATCH": "0", } + if getattr(input_data, "Equation_Class", "") == "Z4C": + defaults["AMSS_CUDA_Z4C_KEEP_RESIDENT_AFTER_STEP"] = "0" + defaults["AMSS_CUDA_KEEP_ALL_LEVELS"] = "0" for key, value in defaults.items(): runtime_env.setdefault(key, value)