diff --git a/.github/workflows/generate-reference-data.yml b/.github/workflows/generate-reference-data.yml index 53f9e2f95..5eb8d1fac 100644 --- a/.github/workflows/generate-reference-data.yml +++ b/.github/workflows/generate-reference-data.yml @@ -30,7 +30,7 @@ jobs: - uses: ursg/gcc-problem-matcher@master - name: Compile vlasiator (Testpackage build) run: | - srun -M carrington --job-name CI_ref_tp_compile --interactive --nodes=1 -n 1 -c 16 --mem=40G -p short -t 0:10:0 bash -c 'module purge; module load GCC/11.2.0; module load OpenMPI/4.1.1-GCC-11.2.0 ; module load PMIx/4.1.0-GCCcore-11.2.0; module load PAPI/6.0.0.1-GCCcore-11.2.0; export VLASIATOR_ARCH=carrington_gcc_openmpi; make -j 16 testpackage; sleep 10s' + srun -M carrington --job-name CI_ref_tp_compile --interactive --nodes=1 -n 1 -c 16 --mem=40G -p short -t 0:10:0 bash -c 'module purge; module load GCC/13.2.0; module load OpenMPI/4.1.6-GCC-13.2.0 ; module load PMIx/4.2.6-GCCcore-13.2.0; module load PAPI/7.1.0-GCCcore-13.2.0; export VLASIATOR_ARCH=carrington_gcc_openmpi; make -j 16 testpackage; sleep 10s' - name: Make sure the output binary is visible in lustre uses: nick-fields/retry@v3 with: @@ -59,7 +59,7 @@ jobs: - name: Checkout source uses: actions/checkout@v4 with: - submodules: true + submodules: false - name: Download testpackage binary uses: actions/download-artifact@v4 with: diff --git a/.github/workflows/github-ci.yml b/.github/workflows/github-ci.yml index 4fb151dc9..ee3c16c9a 100644 --- a/.github/workflows/github-ci.yml +++ b/.github/workflows/github-ci.yml @@ -18,7 +18,7 @@ jobs: # Build libraries for the current version of the docker image # (to be used in any subsequent compilation and run jobs on said image) runs-on: ubuntu-latest - container: ursg/vlasiator_ci:20240131_1 + container: ursg/vlasiator_ci:20241101_1 steps: - name: Checkout source @@ -98,7 +98,7 @@ jobs: build_production: # Build Vlasiator with production flags runs-on: ubuntu-latest - container: ursg/vlasiator_ci:20240131_1 + container: ursg/vlasiator_ci:20241101_1 needs: build_libraries steps: @@ -124,6 +124,64 @@ jobs: path: vlasiator if-no-files-found: error + build_production_1st_order_time: + # Build Vlasiator with production flags + runs-on: ubuntu-latest + container: ursg/vlasiator_ci:20241101_1 + needs: build_libraries + + steps: + - name: Checkout source + uses: actions/checkout@v4 + with: + submodules: true + - name: Download libraries + uses: actions/download-artifact@v4 + with: + name: libraries + - name: Unpack libraries + run: tar --zstd -xvf libraries.tar.zstd + - uses: ursg/gcc-problem-matcher@master + - name: Compile vlasiator (Release build) + run: | + VLASIATOR_ARCH=github_actions make clean + VLASIATOR_ARCH=github_actions COMPFLAGS="-DFS_1ST_ORDER_TIME " make -j 3 + - name: Upload release binary + uses: actions/upload-artifact@v4 + with: + name: vlasiator-1st_order_time + path: vlasiator + if-no-files-found: error + + build_production_1st_order_space: + # Build Vlasiator with production flags + runs-on: ubuntu-latest + container: ursg/vlasiator_ci:20241101_1 + needs: build_libraries + + steps: + - name: Checkout source + uses: actions/checkout@v4 + with: + submodules: true + - name: Download libraries + uses: actions/download-artifact@v4 + with: + name: libraries + - name: Unpack libraries + run: tar --zstd -xvf libraries.tar.zstd + - uses: ursg/gcc-problem-matcher@master + - name: Compile vlasiator (Release build) + run: | + VLASIATOR_ARCH=github_actions make clean + VLASIATOR_ARCH=github_actions COMPFLAGS="-DFS_1ST_ORDER_SPACE " make -j 3 + - name: Upload release binary + uses: actions/upload-artifact@v4 + with: + name: vlasiator-1st_order_space + path: vlasiator + if-no-files-found: error + build_testpackage: # Build Vlasiator with testpackage flags, on the carrington cluster @@ -339,7 +397,7 @@ jobs: build_ionosphereTests: # Build IonosphereSolverTests miniApp runs-on: ubuntu-latest - container: ursg/vlasiator_ci:20240131_1 + container: ursg/vlasiator_ci:20241101_1 needs: build_libraries steps: - name: Checkout source @@ -375,7 +433,7 @@ jobs: run_ionosphere_LFMtest: runs-on: ubuntu-latest - container: ursg/vlasiator_ci:20240131_1 + container: ursg/vlasiator_ci:20241101_1 needs: build_ionosphereTests steps: - name: Checkout source @@ -420,7 +478,7 @@ jobs: check_cfg_files: runs-on: ubuntu-latest - container: ursg/vlasiator_ci:20240131_1 + container: ursg/vlasiator_ci:20241101_1 needs: [build_libraries, build_production] steps: - name: Checkout source diff --git a/build_libraries.sh b/build_libraries.sh index 50b655c47..efd35466a 100755 --- a/build_libraries.sh +++ b/build_libraries.sh @@ -57,7 +57,7 @@ if [[ $PLATFORM != "-arriesgado" && $PLATFORM != "-appleM1" ]]; then # This fai fi # Build jemalloc -wget https://github.com/jemalloc/jemalloc/releases/download/5.3.0/jemalloc-5.3.0.tar.bz2 +curl -O -L https://github.com/jemalloc/jemalloc/releases/download/5.3.0/jemalloc-5.3.0.tar.bz2 tar xjf jemalloc-5.3.0.tar.bz2 cd jemalloc-5.3.0 ./configure --prefix=$WORKSPACE/libraries${PLATFORM} --with-jemalloc-prefix=je_ && make -j 4 && make install diff --git a/common.h b/common.h index 7f35fa601..253f3df4c 100644 --- a/common.h +++ b/common.h @@ -165,15 +165,27 @@ namespace CellParams { P_11, /*!< Pressure P_xx component, computed by Vlasov propagator. */ P_22, /*!< Pressure P_yy component, computed by Vlasov propagator. */ P_33, /*!< Pressure P_zz component, computed by Vlasov propagator. */ + P_23, /*!< Pressure P_yz component, computed by Vlasov propagator. */ + P_13, /*!< Pressure P_xz component, computed by Vlasov propagator. */ + P_12, /*!< Pressure P_xy component, computed by Vlasov propagator. */ P_11_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_11_R and P_11_V. */ P_22_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_22_R and P_22_V. */ P_33_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_33_R and P_33_V. */ + P_23_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_23_R and P_23_V. */ + P_13_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_13_R and P_13_V. */ + P_12_DT2, /*!< Intermediate step value for RK2 time stepping in field solver. Computed from P_12_R and P_12_V. */ P_11_R, /*!< P_xx component after propagation in ordinary space */ P_22_R, /*!< P_yy component after propagation in ordinary space */ P_33_R, /*!< P_zz component after propagation in ordinary space */ + P_23_R, /*!< P_yz component after propagation in ordinary space */ + P_13_R, /*!< P_xz component after propagation in ordinary space */ + P_12_R, /*!< P_xy component after propagation in ordinary space */ P_11_V, /*!< P_xx component after propagation in velocity space */ P_22_V, /*!< P_yy component after propagation in velocity space */ P_33_V, /*!< P_zz component after propagation in velocity space */ + P_23_V, /*!< P_yz component after propagation in velocity space */ + P_12_V, /*!< P_xy component after propagation in velocity space */ + P_13_V, /*!< P_xz component after propagation in velocity space */ EXVOL, /*!< Volume electric field averaged over spatial cell, x-component.*/ EYVOL, /*!< Volume electric field averaged over spatial cell, y-component.*/ EZVOL, /*!< Volume electric field averaged over spatial cell, z-component.*/ @@ -209,11 +221,12 @@ namespace CellParams { AMR_DB, AMR_ALPHA1, AMR_ALPHA2, + P_ANISOTROPY, + AMR_VORTICITY, RECENTLY_REFINED, BULKV_FORCING_X, /*! Externally forced drift velocity (ex. from the ionosphere) */ BULKV_FORCING_Y, /*! Externally forced drift velocity (ex. from the ionosphere) */ BULKV_FORCING_Z, /*! Externally forced drift velocity (ex. from the ionosphere) */ - FORCING_CELL_NUM, /*! Number of boundary cells that have forced a bulkv here */ N_SPATIAL_CELL_PARAMS }; } @@ -236,6 +249,22 @@ namespace bvolderivatives { }; } +namespace vderivatives { + // Essentially a copy from dmoments + enum { + dVxdx, /*!< Derivative of volume-averaged Vx to x-direction. */ + dVxdy, /*!< Derivative of volume-averaged Vx to y-direction. */ + dVxdz, /*!< Derivative of volume-averaged Vx to z-direction. */ + dVydx, /*!< Derivative of volume-averaged Vy to x-direction. */ + dVydy, /*!< Derivative of volume-averaged Vy to y-direction. */ + dVydz, /*!< Derivative of volume-averaged Vy to z-direction. */ + dVzdx, /*!< Derivative of volume-averaged Vz to x-direction. */ + dVzdy, /*!< Derivative of volume-averaged Vz to y-direction. */ + dVzdz, /*!< Derivative of volume-averaged Vz to z-direction. */ + N_V_DERIVATIVES + }; +} + // FsGrid< std::array, FS_STENCIL_WIDTH> & perBGrid, // FsGrid< std::array, FS_STENCIL_WIDTH> & perBDt2Grid, // FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, diff --git a/datareduction/datareducer.cpp b/datareduction/datareducer.cpp index 309915ae9..dda5607d0 100644 --- a/datareduction/datareducer.cpp +++ b/datareduction/datareducer.cpp @@ -2962,7 +2962,21 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti } if(P::systemWriteAllDROs || lowercase == "vg_amr_alpha2") { outputReducer->addOperator(new DRO::DataReductionOperatorCellParams("vg_amr_alpha2",CellParams::AMR_ALPHA2,1)); - outputReducer->addMetadata(outputReducer->size()-1,"","","$\\alpha_2",""); + outputReducer->addMetadata(outputReducer->size()-1,"","","$\\alpha_2$",""); + if(!P::systemWriteAllDROs) { + continue; + } + } + if(P::systemWriteAllDROs || lowercase == "vg_pressure_anisotropy") { + outputReducer->addOperator(new DRO::DataReductionOperatorCellParams("vg_pressure_anisotropy",CellParams::P_ANISOTROPY,1)); + outputReducer->addMetadata(outputReducer->size()-1,"","","$P_\\perp / P_\\parallel$",""); + if(!P::systemWriteAllDROs) { + continue; + } + } + if(P::systemWriteAllDROs || lowercase == "vg_amr_vorticity") { + outputReducer->addOperator(new DRO::DataReductionOperatorCellParams("vg_amr_vorticity",CellParams::AMR_VORTICITY,1)); + outputReducer->addMetadata(outputReducer->size()-1,"","","Vorticity",""); if(!P::systemWriteAllDROs) { continue; } diff --git a/fieldsolver/derivatives.cpp b/fieldsolver/derivatives.cpp index b1aad577c..8112d08ac 100644 --- a/fieldsolver/derivatives.cpp +++ b/fieldsolver/derivatives.cpp @@ -25,6 +25,7 @@ #include "fs_common.h" #include "derivatives.hpp" #include "fs_limiters.h" +#include /*! \brief Low-level spatial derivatives calculation. * @@ -43,7 +44,7 @@ * \param dMomentsGrid fsGrid holding the derviatives of moments * \param technicalGrid fsGrid holding technical information (such as boundary types) * \param sysBoundaries System boundary conditions existing - * \param RKCase Element in the enum defining the Runge-Kutta method steps + * \param calculateMoments Bool telling whether the derivatives for moments need updating too. * * \sa calculateDerivativesSimple calculateBVOLDerivativesSimple calculateBVOLDerivatives */ @@ -57,7 +58,7 @@ void calculateDerivatives( FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, SysBoundary& sysBoundaries, - cint& RKCase + const bool calculateMoments ) { std::array * dPerB = dPerBGrid.get(i,j,k); std::array * dMoments = dMomentsGrid.get(i,j,k); @@ -119,31 +120,37 @@ void calculateDerivatives( #endif if(sysBoundaryLayer == 1 || (sysBoundaryLayer == 2 && sysBoundaryFlag == sysboundarytype::NOT_SYSBOUNDARY)) { - dMoments->at(fsgrids::dmoments::drhomdx) = (rghtMoments->at(fsgrids::moments::RHOM)-leftMoments->at(fsgrids::moments::RHOM))/2; - dMoments->at(fsgrids::dmoments::drhoqdx) = (rghtMoments->at(fsgrids::moments::RHOQ)-leftMoments->at(fsgrids::moments::RHOQ))/2; - dMoments->at(fsgrids::dmoments::dp11dx) = (rghtMoments->at(fsgrids::moments::P_11)-leftMoments->at(fsgrids::moments::P_11))/2; - dMoments->at(fsgrids::dmoments::dp22dx) = (rghtMoments->at(fsgrids::moments::P_22)-leftMoments->at(fsgrids::moments::P_22))/2; - dMoments->at(fsgrids::dmoments::dp33dx) = (rghtMoments->at(fsgrids::moments::P_33)-leftMoments->at(fsgrids::moments::P_33))/2; - dMoments->at(fsgrids::dmoments::dVxdx) = (rghtMoments->at(fsgrids::moments::VX)-leftMoments->at(fsgrids::moments::VX))/2; - dMoments->at(fsgrids::dmoments::dVydx) = (rghtMoments->at(fsgrids::moments::VY)-leftMoments->at(fsgrids::moments::VY))/2; - dMoments->at(fsgrids::dmoments::dVzdx) = (rghtMoments->at(fsgrids::moments::VZ)-leftMoments->at(fsgrids::moments::VZ))/2; + if(calculateMoments) { + dMoments->at(fsgrids::dmoments::drhomdx) = (rghtMoments->at(fsgrids::moments::RHOM)-leftMoments->at(fsgrids::moments::RHOM))/2; + dMoments->at(fsgrids::dmoments::drhoqdx) = (rghtMoments->at(fsgrids::moments::RHOQ)-leftMoments->at(fsgrids::moments::RHOQ))/2; + dMoments->at(fsgrids::dmoments::dp11dx) = (rghtMoments->at(fsgrids::moments::P_11)-leftMoments->at(fsgrids::moments::P_11))/2; + dMoments->at(fsgrids::dmoments::dp22dx) = (rghtMoments->at(fsgrids::moments::P_22)-leftMoments->at(fsgrids::moments::P_22))/2; + dMoments->at(fsgrids::dmoments::dp33dx) = (rghtMoments->at(fsgrids::moments::P_33)-leftMoments->at(fsgrids::moments::P_33))/2; + dMoments->at(fsgrids::dmoments::dVxdx) = (rghtMoments->at(fsgrids::moments::VX)-leftMoments->at(fsgrids::moments::VX))/2; + dMoments->at(fsgrids::dmoments::dVydx) = (rghtMoments->at(fsgrids::moments::VY)-leftMoments->at(fsgrids::moments::VY))/2; + dMoments->at(fsgrids::dmoments::dVzdx) = (rghtMoments->at(fsgrids::moments::VZ)-leftMoments->at(fsgrids::moments::VZ))/2; + } dPerB->at(fsgrids::dperb::dPERBydx) = (rghtPerB->at(fsgrids::bfield::PERBY)-leftPerB->at(fsgrids::bfield::PERBY))/2; dPerB->at(fsgrids::dperb::dPERBzdx) = (rghtPerB->at(fsgrids::bfield::PERBZ)-leftPerB->at(fsgrids::bfield::PERBZ))/2; } else { - dMoments->at(fsgrids::dmoments::drhomdx) = limiter(leftMoments->at(fsgrids::moments::RHOM),centMoments->at(fsgrids::moments::RHOM),rghtMoments->at(fsgrids::moments::RHOM)); - dMoments->at(fsgrids::dmoments::drhoqdx) = limiter(leftMoments->at(fsgrids::moments::RHOQ),centMoments->at(fsgrids::moments::RHOQ),rghtMoments->at(fsgrids::moments::RHOQ)); - dMoments->at(fsgrids::dmoments::dp11dx) = limiter(leftMoments->at(fsgrids::moments::P_11),centMoments->at(fsgrids::moments::P_11),rghtMoments->at(fsgrids::moments::P_11)); - dMoments->at(fsgrids::dmoments::dp22dx) = limiter(leftMoments->at(fsgrids::moments::P_22),centMoments->at(fsgrids::moments::P_22),rghtMoments->at(fsgrids::moments::P_22)); - dMoments->at(fsgrids::dmoments::dp33dx) = limiter(leftMoments->at(fsgrids::moments::P_33),centMoments->at(fsgrids::moments::P_33),rghtMoments->at(fsgrids::moments::P_33)); - dMoments->at(fsgrids::dmoments::dVxdx) = limiter(leftMoments->at(fsgrids::moments::VX), centMoments->at(fsgrids::moments::VX), rghtMoments->at(fsgrids::moments::VX)); - dMoments->at(fsgrids::dmoments::dVydx) = limiter(leftMoments->at(fsgrids::moments::VY), centMoments->at(fsgrids::moments::VY), rghtMoments->at(fsgrids::moments::VY)); - dMoments->at(fsgrids::dmoments::dVzdx) = limiter(leftMoments->at(fsgrids::moments::VZ), centMoments->at(fsgrids::moments::VZ), rghtMoments->at(fsgrids::moments::VZ)); + if(calculateMoments) { + dMoments->at(fsgrids::dmoments::drhomdx) = limiter(leftMoments->at(fsgrids::moments::RHOM),centMoments->at(fsgrids::moments::RHOM),rghtMoments->at(fsgrids::moments::RHOM)); + dMoments->at(fsgrids::dmoments::drhoqdx) = limiter(leftMoments->at(fsgrids::moments::RHOQ),centMoments->at(fsgrids::moments::RHOQ),rghtMoments->at(fsgrids::moments::RHOQ)); + dMoments->at(fsgrids::dmoments::dp11dx) = limiter(leftMoments->at(fsgrids::moments::P_11),centMoments->at(fsgrids::moments::P_11),rghtMoments->at(fsgrids::moments::P_11)); + dMoments->at(fsgrids::dmoments::dp22dx) = limiter(leftMoments->at(fsgrids::moments::P_22),centMoments->at(fsgrids::moments::P_22),rghtMoments->at(fsgrids::moments::P_22)); + dMoments->at(fsgrids::dmoments::dp33dx) = limiter(leftMoments->at(fsgrids::moments::P_33),centMoments->at(fsgrids::moments::P_33),rghtMoments->at(fsgrids::moments::P_33)); + dMoments->at(fsgrids::dmoments::dVxdx) = limiter(leftMoments->at(fsgrids::moments::VX), centMoments->at(fsgrids::moments::VX), rghtMoments->at(fsgrids::moments::VX)); + dMoments->at(fsgrids::dmoments::dVydx) = limiter(leftMoments->at(fsgrids::moments::VY), centMoments->at(fsgrids::moments::VY), rghtMoments->at(fsgrids::moments::VY)); + dMoments->at(fsgrids::dmoments::dVzdx) = limiter(leftMoments->at(fsgrids::moments::VZ), centMoments->at(fsgrids::moments::VZ), rghtMoments->at(fsgrids::moments::VZ)); + } dPerB->at(fsgrids::dperb::dPERBydx) = limiter(leftPerB->at(fsgrids::bfield::PERBY),centPerB->at(fsgrids::bfield::PERBY),rghtPerB->at(fsgrids::bfield::PERBY)); dPerB->at(fsgrids::dperb::dPERBzdx) = limiter(leftPerB->at(fsgrids::bfield::PERBZ),centPerB->at(fsgrids::bfield::PERBZ),rghtPerB->at(fsgrids::bfield::PERBZ)); } - // pres_e = const * np.power(rho_e, index) - dMoments->at(fsgrids::dmoments::dPedx) = Peconst * limiter(pow(leftMoments->at(fsgrids::moments::RHOQ)/physicalconstants::CHARGE,Parameters::electronPTindex),pow(centMoments->at(fsgrids::moments::RHOQ)/physicalconstants::CHARGE,Parameters::electronPTindex),pow(rghtMoments->at(fsgrids::moments::RHOQ)/physicalconstants::CHARGE,Parameters::electronPTindex)); + if(calculateMoments) { + // pres_e = const * np.power(rho_e, index) + dMoments->at(fsgrids::dmoments::dPedx) = Peconst * limiter(pow(leftMoments->at(fsgrids::moments::RHOQ)/physicalconstants::CHARGE,Parameters::electronPTindex),pow(centMoments->at(fsgrids::moments::RHOQ)/physicalconstants::CHARGE,Parameters::electronPTindex),pow(rghtMoments->at(fsgrids::moments::RHOQ)/physicalconstants::CHARGE,Parameters::electronPTindex)); + } if (Parameters::ohmHallTerm < 2 || sysBoundaryLayer == 1) { dPerB->at(fsgrids::dperb::dPERBydxx) = 0.0; @@ -168,25 +175,29 @@ void calculateDerivatives( } if(sysBoundaryLayer == 1 || (sysBoundaryLayer == 2 && sysBoundaryFlag == sysboundarytype::NOT_SYSBOUNDARY)) { - dMoments->at(fsgrids::dmoments::drhomdy) = (rghtMoments->at(fsgrids::moments::RHOM)-leftMoments->at(fsgrids::moments::RHOM))/2; - dMoments->at(fsgrids::dmoments::drhoqdy) = (rghtMoments->at(fsgrids::moments::RHOQ)-leftMoments->at(fsgrids::moments::RHOQ))/2; - dMoments->at(fsgrids::dmoments::dp11dy) = (rghtMoments->at(fsgrids::moments::P_11)-leftMoments->at(fsgrids::moments::P_11))/2; - dMoments->at(fsgrids::dmoments::dp22dy) = (rghtMoments->at(fsgrids::moments::P_22)-leftMoments->at(fsgrids::moments::P_22))/2; - dMoments->at(fsgrids::dmoments::dp33dy) = (rghtMoments->at(fsgrids::moments::P_33)-leftMoments->at(fsgrids::moments::P_33))/2; - dMoments->at(fsgrids::dmoments::dVxdy) = (rghtMoments->at(fsgrids::moments::VX)-leftMoments->at(fsgrids::moments::VX))/2; - dMoments->at(fsgrids::dmoments::dVydy) = (rghtMoments->at(fsgrids::moments::VY)-leftMoments->at(fsgrids::moments::VY))/2; - dMoments->at(fsgrids::dmoments::dVzdy) = (rghtMoments->at(fsgrids::moments::VZ)-leftMoments->at(fsgrids::moments::VZ))/2; + if(calculateMoments) { + dMoments->at(fsgrids::dmoments::drhomdy) = (rghtMoments->at(fsgrids::moments::RHOM)-leftMoments->at(fsgrids::moments::RHOM))/2; + dMoments->at(fsgrids::dmoments::drhoqdy) = (rghtMoments->at(fsgrids::moments::RHOQ)-leftMoments->at(fsgrids::moments::RHOQ))/2; + dMoments->at(fsgrids::dmoments::dp11dy) = (rghtMoments->at(fsgrids::moments::P_11)-leftMoments->at(fsgrids::moments::P_11))/2; + dMoments->at(fsgrids::dmoments::dp22dy) = (rghtMoments->at(fsgrids::moments::P_22)-leftMoments->at(fsgrids::moments::P_22))/2; + dMoments->at(fsgrids::dmoments::dp33dy) = (rghtMoments->at(fsgrids::moments::P_33)-leftMoments->at(fsgrids::moments::P_33))/2; + dMoments->at(fsgrids::dmoments::dVxdy) = (rghtMoments->at(fsgrids::moments::VX)-leftMoments->at(fsgrids::moments::VX))/2; + dMoments->at(fsgrids::dmoments::dVydy) = (rghtMoments->at(fsgrids::moments::VY)-leftMoments->at(fsgrids::moments::VY))/2; + dMoments->at(fsgrids::dmoments::dVzdy) = (rghtMoments->at(fsgrids::moments::VZ)-leftMoments->at(fsgrids::moments::VZ))/2; + } dPerB->at(fsgrids::dperb::dPERBxdy) = (rghtPerB->at(fsgrids::bfield::PERBX)-leftPerB->at(fsgrids::bfield::PERBX))/2; dPerB->at(fsgrids::dperb::dPERBzdy) = (rghtPerB->at(fsgrids::bfield::PERBZ)-leftPerB->at(fsgrids::bfield::PERBZ))/2; } else { - dMoments->at(fsgrids::dmoments::drhomdy) = limiter(leftMoments->at(fsgrids::moments::RHOM),centMoments->at(fsgrids::moments::RHOM),rghtMoments->at(fsgrids::moments::RHOM)); - dMoments->at(fsgrids::dmoments::drhoqdy) = limiter(leftMoments->at(fsgrids::moments::RHOQ),centMoments->at(fsgrids::moments::RHOQ),rghtMoments->at(fsgrids::moments::RHOQ)); - dMoments->at(fsgrids::dmoments::dp11dy) = limiter(leftMoments->at(fsgrids::moments::P_11),centMoments->at(fsgrids::moments::P_11),rghtMoments->at(fsgrids::moments::P_11)); - dMoments->at(fsgrids::dmoments::dp22dy) = limiter(leftMoments->at(fsgrids::moments::P_22),centMoments->at(fsgrids::moments::P_22),rghtMoments->at(fsgrids::moments::P_22)); - dMoments->at(fsgrids::dmoments::dp33dy) = limiter(leftMoments->at(fsgrids::moments::P_33),centMoments->at(fsgrids::moments::P_33),rghtMoments->at(fsgrids::moments::P_33)); - dMoments->at(fsgrids::dmoments::dVxdy) = limiter(leftMoments->at(fsgrids::moments::VX), centMoments->at(fsgrids::moments::VX), rghtMoments->at(fsgrids::moments::VX)); - dMoments->at(fsgrids::dmoments::dVydy) = limiter(leftMoments->at(fsgrids::moments::VY), centMoments->at(fsgrids::moments::VY), rghtMoments->at(fsgrids::moments::VY)); - dMoments->at(fsgrids::dmoments::dVzdy) = limiter(leftMoments->at(fsgrids::moments::VZ), centMoments->at(fsgrids::moments::VZ), rghtMoments->at(fsgrids::moments::VZ)); + if(calculateMoments) { + dMoments->at(fsgrids::dmoments::drhomdy) = limiter(leftMoments->at(fsgrids::moments::RHOM),centMoments->at(fsgrids::moments::RHOM),rghtMoments->at(fsgrids::moments::RHOM)); + dMoments->at(fsgrids::dmoments::drhoqdy) = limiter(leftMoments->at(fsgrids::moments::RHOQ),centMoments->at(fsgrids::moments::RHOQ),rghtMoments->at(fsgrids::moments::RHOQ)); + dMoments->at(fsgrids::dmoments::dp11dy) = limiter(leftMoments->at(fsgrids::moments::P_11),centMoments->at(fsgrids::moments::P_11),rghtMoments->at(fsgrids::moments::P_11)); + dMoments->at(fsgrids::dmoments::dp22dy) = limiter(leftMoments->at(fsgrids::moments::P_22),centMoments->at(fsgrids::moments::P_22),rghtMoments->at(fsgrids::moments::P_22)); + dMoments->at(fsgrids::dmoments::dp33dy) = limiter(leftMoments->at(fsgrids::moments::P_33),centMoments->at(fsgrids::moments::P_33),rghtMoments->at(fsgrids::moments::P_33)); + dMoments->at(fsgrids::dmoments::dVxdy) = limiter(leftMoments->at(fsgrids::moments::VX), centMoments->at(fsgrids::moments::VX), rghtMoments->at(fsgrids::moments::VX)); + dMoments->at(fsgrids::dmoments::dVydy) = limiter(leftMoments->at(fsgrids::moments::VY), centMoments->at(fsgrids::moments::VY), rghtMoments->at(fsgrids::moments::VY)); + dMoments->at(fsgrids::dmoments::dVzdy) = limiter(leftMoments->at(fsgrids::moments::VZ), centMoments->at(fsgrids::moments::VZ), rghtMoments->at(fsgrids::moments::VZ)); + } dPerB->at(fsgrids::dperb::dPERBxdy) = limiter(leftPerB->at(fsgrids::bfield::PERBX),centPerB->at(fsgrids::bfield::PERBX),rghtPerB->at(fsgrids::bfield::PERBX)); dPerB->at(fsgrids::dperb::dPERBzdy) = limiter(leftPerB->at(fsgrids::bfield::PERBZ),centPerB->at(fsgrids::bfield::PERBZ),rghtPerB->at(fsgrids::bfield::PERBZ)); } @@ -216,31 +227,38 @@ void calculateDerivatives( } if(sysBoundaryLayer == 1 || (sysBoundaryLayer == 2 && sysBoundaryFlag == sysboundarytype::NOT_SYSBOUNDARY)) { - dMoments->at(fsgrids::dmoments::drhomdz) = (rghtMoments->at(fsgrids::moments::RHOM)-leftMoments->at(fsgrids::moments::RHOM))/2; - dMoments->at(fsgrids::dmoments::drhoqdz) = (rghtMoments->at(fsgrids::moments::RHOQ)-leftMoments->at(fsgrids::moments::RHOQ))/2; - dMoments->at(fsgrids::dmoments::dp11dz) = (rghtMoments->at(fsgrids::moments::P_11)-leftMoments->at(fsgrids::moments::P_11))/2; - dMoments->at(fsgrids::dmoments::dp22dz) = (rghtMoments->at(fsgrids::moments::P_22)-leftMoments->at(fsgrids::moments::P_22))/2; - dMoments->at(fsgrids::dmoments::dp33dz) = (rghtMoments->at(fsgrids::moments::P_33)-leftMoments->at(fsgrids::moments::P_33))/2; - dMoments->at(fsgrids::dmoments::dVxdz) = (rghtMoments->at(fsgrids::moments::VX)-leftMoments->at(fsgrids::moments::VX))/2; - dMoments->at(fsgrids::dmoments::dVydz) = (rghtMoments->at(fsgrids::moments::VY)-leftMoments->at(fsgrids::moments::VY))/2; - dMoments->at(fsgrids::dmoments::dVzdz) = (rghtMoments->at(fsgrids::moments::VZ)-leftMoments->at(fsgrids::moments::VZ))/2; + if(calculateMoments) { + dMoments->at(fsgrids::dmoments::drhomdz) = (rghtMoments->at(fsgrids::moments::RHOM)-leftMoments->at(fsgrids::moments::RHOM))/2; + dMoments->at(fsgrids::dmoments::drhoqdz) = (rghtMoments->at(fsgrids::moments::RHOQ)-leftMoments->at(fsgrids::moments::RHOQ))/2; + dMoments->at(fsgrids::dmoments::dp11dz) = (rghtMoments->at(fsgrids::moments::P_11)-leftMoments->at(fsgrids::moments::P_11))/2; + dMoments->at(fsgrids::dmoments::dp22dz) = (rghtMoments->at(fsgrids::moments::P_22)-leftMoments->at(fsgrids::moments::P_22))/2; + dMoments->at(fsgrids::dmoments::dp33dz) = (rghtMoments->at(fsgrids::moments::P_33)-leftMoments->at(fsgrids::moments::P_33))/2; + dMoments->at(fsgrids::dmoments::dVxdz) = (rghtMoments->at(fsgrids::moments::VX)-leftMoments->at(fsgrids::moments::VX))/2; + dMoments->at(fsgrids::dmoments::dVydz) = (rghtMoments->at(fsgrids::moments::VY)-leftMoments->at(fsgrids::moments::VY))/2; + dMoments->at(fsgrids::dmoments::dVzdz) = (rghtMoments->at(fsgrids::moments::VZ)-leftMoments->at(fsgrids::moments::VZ))/2; + } dPerB->at(fsgrids::dperb::dPERBxdz) = (rghtPerB->at(fsgrids::bfield::PERBX)-leftPerB->at(fsgrids::bfield::PERBX))/2; dPerB->at(fsgrids::dperb::dPERBydz) = (rghtPerB->at(fsgrids::bfield::PERBY)-leftPerB->at(fsgrids::bfield::PERBY))/2; } else { - dMoments->at(fsgrids::dmoments::drhomdz) = limiter(leftMoments->at(fsgrids::moments::RHOM),centMoments->at(fsgrids::moments::RHOM),rghtMoments->at(fsgrids::moments::RHOM)); - dMoments->at(fsgrids::dmoments::drhoqdz) = limiter(leftMoments->at(fsgrids::moments::RHOQ),centMoments->at(fsgrids::moments::RHOQ),rghtMoments->at(fsgrids::moments::RHOQ)); - dMoments->at(fsgrids::dmoments::dp11dz) = limiter(leftMoments->at(fsgrids::moments::P_11),centMoments->at(fsgrids::moments::P_11),rghtMoments->at(fsgrids::moments::P_11)); - dMoments->at(fsgrids::dmoments::dp22dz) = limiter(leftMoments->at(fsgrids::moments::P_22),centMoments->at(fsgrids::moments::P_22),rghtMoments->at(fsgrids::moments::P_22)); - dMoments->at(fsgrids::dmoments::dp33dz) = limiter(leftMoments->at(fsgrids::moments::P_33),centMoments->at(fsgrids::moments::P_33),rghtMoments->at(fsgrids::moments::P_33)); - dMoments->at(fsgrids::dmoments::dVxdz) = limiter(leftMoments->at(fsgrids::moments::VX), centMoments->at(fsgrids::moments::VX), rghtMoments->at(fsgrids::moments::VX)); - dMoments->at(fsgrids::dmoments::dVydz) = limiter(leftMoments->at(fsgrids::moments::VY), centMoments->at(fsgrids::moments::VY), rghtMoments->at(fsgrids::moments::VY)); - dMoments->at(fsgrids::dmoments::dVzdz) = limiter(leftMoments->at(fsgrids::moments::VZ), centMoments->at(fsgrids::moments::VZ), rghtMoments->at(fsgrids::moments::VZ)); + if(calculateMoments) { + dMoments->at(fsgrids::dmoments::drhomdz) = limiter(leftMoments->at(fsgrids::moments::RHOM),centMoments->at(fsgrids::moments::RHOM),rghtMoments->at(fsgrids::moments::RHOM)); + dMoments->at(fsgrids::dmoments::drhoqdz) = limiter(leftMoments->at(fsgrids::moments::RHOQ),centMoments->at(fsgrids::moments::RHOQ),rghtMoments->at(fsgrids::moments::RHOQ)); + dMoments->at(fsgrids::dmoments::dp11dz) = limiter(leftMoments->at(fsgrids::moments::P_11),centMoments->at(fsgrids::moments::P_11),rghtMoments->at(fsgrids::moments::P_11)); + dMoments->at(fsgrids::dmoments::dp22dz) = limiter(leftMoments->at(fsgrids::moments::P_22),centMoments->at(fsgrids::moments::P_22),rghtMoments->at(fsgrids::moments::P_22)); + dMoments->at(fsgrids::dmoments::dp33dz) = limiter(leftMoments->at(fsgrids::moments::P_33),centMoments->at(fsgrids::moments::P_33),rghtMoments->at(fsgrids::moments::P_33)); + dMoments->at(fsgrids::dmoments::dVxdz) = limiter(leftMoments->at(fsgrids::moments::VX), centMoments->at(fsgrids::moments::VX), rghtMoments->at(fsgrids::moments::VX)); + dMoments->at(fsgrids::dmoments::dVydz) = limiter(leftMoments->at(fsgrids::moments::VY), centMoments->at(fsgrids::moments::VY), rghtMoments->at(fsgrids::moments::VY)); + dMoments->at(fsgrids::dmoments::dVzdz) = limiter(leftMoments->at(fsgrids::moments::VZ), centMoments->at(fsgrids::moments::VZ), rghtMoments->at(fsgrids::moments::VZ)); + } dPerB->at(fsgrids::dperb::dPERBxdz) = limiter(leftPerB->at(fsgrids::bfield::PERBX),centPerB->at(fsgrids::bfield::PERBX),rghtPerB->at(fsgrids::bfield::PERBX)); dPerB->at(fsgrids::dperb::dPERBydz) = limiter(leftPerB->at(fsgrids::bfield::PERBY),centPerB->at(fsgrids::bfield::PERBY),rghtPerB->at(fsgrids::bfield::PERBY)); } - // pres_e = const * np.power(rho_e, index) - dMoments->at(fsgrids::dmoments::dPedz) = Peconst * limiter(pow(leftMoments->at(fsgrids::moments::RHOQ)/physicalconstants::CHARGE,Parameters::electronPTindex),pow(centMoments->at(fsgrids::moments::RHOQ)/physicalconstants::CHARGE,Parameters::electronPTindex),pow(rghtMoments->at(fsgrids::moments::RHOQ)/physicalconstants::CHARGE,Parameters::electronPTindex)); + if(calculateMoments) { + // pres_e = const * np.power(rho_e, index) + dMoments->at(fsgrids::dmoments::dPedz) = Peconst * limiter(pow(leftMoments->at(fsgrids::moments::RHOQ)/physicalconstants::CHARGE,Parameters::electronPTindex),pow(centMoments->at(fsgrids::moments::RHOQ)/physicalconstants::CHARGE,Parameters::electronPTindex),pow(rghtMoments->at(fsgrids::moments::RHOQ)/physicalconstants::CHARGE,Parameters::electronPTindex)); + } + if (Parameters::ohmHallTerm < 2 || sysBoundaryLayer == 1) { dPerB->at(fsgrids::dperb::dPERBxdzz) = 0.0; dPerB->at(fsgrids::dperb::dPERBydzz) = 0.0; @@ -303,6 +321,7 @@ void calculateDerivatives( * \param momentsDt2Grid fsGrid holding the moment quantities at runge-kutta t=0.5 * \param dPerBGrid fsGrid holding the derivatives of perturbed B * \param dMomentsGrid fsGrid holding the derviatives of moments + * \param dMomentsGridDt2 fsGrid holding the derviatives of moments at runge-kutta t=0.5 * \param technicalGrid fsGrid holding technical information (such as boundary types) * \param sysBoundaries System boundary conditions existing * \param RKCase Element in the enum defining the Runge-Kutta method steps @@ -317,10 +336,11 @@ void calculateDerivativesSimple( FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsDt2Grid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, SysBoundary& sysBoundaries, cint& RKCase, - const bool communicateMoments) { + const bool doMoments) { //const std::array gridDims = technicalGrid.getLocalSize(); const FsGridTools::FsIndex_t* gridDims = &technicalGrid.getLocalSize()[0]; const size_t N_cells = gridDims[0]*gridDims[1]*gridDims[2]; @@ -335,7 +355,7 @@ void calculateDerivativesSimple( // The update of PERB[XYZ] is needed after the system // boundary update of propagateMagneticFieldSimple. perBGrid.updateGhostCells(); - if(communicateMoments) { + if(doMoments) { momentsGrid.updateGhostCells(); } break; @@ -344,7 +364,7 @@ void calculateDerivativesSimple( // update of PERB[XYZ]_DT2 is needed after the system // boundary update of propagateMagneticFieldSimple. perBDt2Grid.updateGhostCells(); - if(communicateMoments) { + if(doMoments) { momentsDt2Grid.updateGhostCells(); } break; @@ -353,7 +373,7 @@ void calculateDerivativesSimple( // is needed after the system boundary update of // propagateMagneticFieldSimple. perBGrid.updateGhostCells(); - if(communicateMoments) { + if(doMoments) { momentsGrid.updateGhostCells(); } break; @@ -372,9 +392,9 @@ void calculateDerivativesSimple( for (FsGridTools::FsIndex_t j=0; jparameters[CellParams::RHOM]}; Real myU {calculateU(cell)}; + Real myV {std::sqrt(std::pow(cell->parameters[CellParams::VX], 2) + std::pow(cell->parameters[CellParams::VY], 2) + std::pow(cell->parameters[CellParams::VZ], 2))}; + Real maxV {myV}; std::array myP = getMomentumDensity(cell); std::array myB = getBVol(cell); for (SpatialCell* neighbor : neighbors) { Real otherRho = neighbor->parameters[CellParams::RHOM]; Real otherU = calculateU(neighbor); + Real otherV {std::sqrt(std::pow(neighbor->parameters[CellParams::VX], 2) + std::pow(neighbor->parameters[CellParams::VY], 2) + std::pow(neighbor->parameters[CellParams::VZ], 2))}; std::array otherP = getMomentumDensity(neighbor); std::array otherB = getBVol(neighbor); Real deltaBsq = pow(myB[0] - otherB[0], 2) + pow(myB[1] - otherB[1], 2) + pow(myB[2] - otherB[2], 2); + if (myV < EPS) { + maxV = std::max(maxV,otherV); + } + Real maxRho = std::max(myRho, otherRho); if (maxRho > EPS) { dRho = std::max(fabs(myRho - otherRho) / maxRho, dRho); @@ -801,6 +828,38 @@ void calculateScaledDeltas( Bperp = std::sqrt(Bperp); } + // Now, rotation matrix to get parallel and perpendicular pressure + //Eigen::Quaterniond q {Quaterniond::FromTwoVectors(Eigen::vector3d{0, 0, 1}, Eigen::vector3d{myB[0], myB[1], myB[2]})}; + //Eigen::Matrix3d rot = q.toRotationMatrix(); + Eigen::Matrix3d rot = Eigen::Quaterniond::FromTwoVectors(Eigen::Vector3d{myB[0], myB[1], myB[2]}, Eigen::Vector3d{0, 0, 1}).normalized().toRotationMatrix(); + Eigen::Matrix3d P { + {cell->parameters[CellParams::P_11], cell->parameters[CellParams::P_12], cell->parameters[CellParams::P_13]}, + {cell->parameters[CellParams::P_12], cell->parameters[CellParams::P_22], cell->parameters[CellParams::P_23]}, + {cell->parameters[CellParams::P_13], cell->parameters[CellParams::P_23], cell->parameters[CellParams::P_33]}, + }; + + Eigen::Matrix3d transposerot = rot.transpose(); + Eigen::Matrix3d Pprime = rot * P * transposerot; + + Real Panisotropy {0.0}; + if (Pprime(2, 2) > EPS) { + Panisotropy = (Pprime(0, 0) + Pprime(1, 1)) / (2 * Pprime(2, 2)); + } + + // Vorticity + Real dVxdy {cell->derivativesV[vderivatives::dVxdy]}; + Real dVxdz {cell->derivativesV[vderivatives::dVxdz]}; + Real dVydx {cell->derivativesV[vderivatives::dVydx]}; + Real dVydz {cell->derivativesV[vderivatives::dVydz]}; + Real dVzdx {cell->derivativesV[vderivatives::dVzdx]}; + Real dVzdy {cell->derivativesV[vderivatives::dVzdy]}; + Real vorticity {std::sqrt(std::pow(dVxdy - dVydz, 2) + std::pow(dVxdz - dVzdx, 2 ) + std::pow(dVydx - dVxdy, 2))}; + //Real vA {std::sqrt(Bsq / (physicalconstants::MU_0 * myRho))}; + Real amr_vorticity {0.0}; + if (maxV > EPS) { + amr_vorticity = vorticity * cell->parameters[CellParams::DX] / maxV; + } + cell->parameters[CellParams::AMR_DRHO] = dRho; cell->parameters[CellParams::AMR_DU] = dU; cell->parameters[CellParams::AMR_DPSQ] = dPsq; @@ -808,6 +867,9 @@ void calculateScaledDeltas( cell->parameters[CellParams::AMR_DB] = dB; cell->parameters[CellParams::AMR_ALPHA1] = alpha; cell->parameters[CellParams::AMR_ALPHA2] = cell->parameters[CellParams::DX] * J / (Bperp + EPS); // Epsilon in denominator so we don't get infinities + cell->parameters[CellParams::P_ANISOTROPY] = Panisotropy; + // Experimental, current scaling is bulk velocity + cell->parameters[CellParams::AMR_VORTICITY] = amr_vorticity; } /*! \brief High-level scaled gradient calculation wrapper function. diff --git a/fieldsolver/derivatives.hpp b/fieldsolver/derivatives.hpp index a076f92a5..0f1000fb2 100644 --- a/fieldsolver/derivatives.hpp +++ b/fieldsolver/derivatives.hpp @@ -36,6 +36,7 @@ void calculateDerivativesSimple( FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsDt2Grid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, SysBoundary& sysBoundaries, cint& RKCase, diff --git a/fieldsolver/fs_common.h b/fieldsolver/fs_common.h index c2704086d..6e7242cd8 100644 --- a/fieldsolver/fs_common.h +++ b/fieldsolver/fs_common.h @@ -62,27 +62,6 @@ static creal EPS = 1.0e-30; using namespace std; -bool initializeFieldPropagator( - FsGrid< std::array, FS_STENCIL_WIDTH> & perBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & perBDt2Grid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EDt2Grid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EHallGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & volGrid, - FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, - SysBoundary& sysBoundaries -); - -bool initializeFieldPropagatorAfterRebalance(); - -bool finalizeFieldPropagator(); - bool propagateFields( FsGrid< std::array, FS_STENCIL_WIDTH> & perBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & perBDt2Grid, @@ -90,10 +69,12 @@ bool propagateFields( FsGrid< std::array, FS_STENCIL_WIDTH> & EDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & EHallGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & volGrid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, diff --git a/fieldsolver/gridGlue.cpp b/fieldsolver/gridGlue.cpp index 7501588ad..4574d1e67 100644 --- a/fieldsolver/gridGlue.cpp +++ b/fieldsolver/gridGlue.cpp @@ -43,9 +43,9 @@ void filterMoments(dccrg::Dccrg& mpiGrid, // Kernel Characteristics - const int kernelOffset = 2; // offset of 5 pointstencil 3D kernel => (floor(stencilWidth/2);) - const Real inverseKernelSum = 1.0 / 729.0; // the inverse of the total kernel's sum - const static Real kernel[5][5][5] ={ + constexpr int kernelOffset = 2; // offset of 5 pointstencil 3D kernel => (floor(stencilWidth/2);) + constexpr Real inverseKernelSum = 1.0 / 729.0; // the inverse of the total kernel's sum + constexpr static Real kernel[5][5][5] ={ {{ 1 * inverseKernelSum, 2 * inverseKernelSum, 3 * inverseKernelSum, 2 * inverseKernelSum, 1 * inverseKernelSum}, { 2 * inverseKernelSum, 4 * inverseKernelSum, 6 * inverseKernelSum, 4 * inverseKernelSum, 2 * inverseKernelSum}, { 3 * inverseKernelSum, 6 * inverseKernelSum, 9 * inverseKernelSum, 6 * inverseKernelSum, 3 * inverseKernelSum}, @@ -94,31 +94,27 @@ void filterMoments(dccrg::Dccrg& mpiGrid, for (FsGridTools::FsIndex_t j = 0; j < mntDims[1]; j++){ for (FsGridTools::FsIndex_t i = 0; i < mntDims[0]; i++){ - // Get refLevel level int refLevel = technicalGrid.get(i, j, k)->refLevel; + auto* swap {swapGrid.get(i, j, k)}; - // Skip pass + // Skip pass and copy value if (blurPass >= P::numPasses.at(refLevel) || technicalGrid.get(i, j, k)->sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY) { + *swap = *momentsGrid.get(i, j, k); continue; } - // Get pointers to our cells - std::array *cell; - std::array *swap; - - // Set Cell to zero before passing filter - swap = swapGrid.get(i,j,k); + // Set moments to zero before passing filter for (int e = 0; e < fsgrids::moments::N_MOMENTS; ++e) { - swap->at(e)=0.0; + swap->at(e) = 0.0; } // Perform the blur for (int c=-kernelOffset; c<=kernelOffset; c++){ for (int b=-kernelOffset; b<=kernelOffset; b++){ for (int a=-kernelOffset; a<=kernelOffset; a++){ - cell = momentsGrid.get(i+a,j+b,k+c); + const auto* cell {momentsGrid.get(i+a,j+b,k+c)}; for (int e = 0; e < fsgrids::moments::N_MOMENTS; ++e) { - swap->at(e)+=cell->at(e) *kernel[kernelOffset+a][kernelOffset+b][kernelOffset+c]; + swap->at(e) += cell->at(e) * kernel[kernelOffset+a][kernelOffset+b][kernelOffset+c]; } } } @@ -127,12 +123,12 @@ void filterMoments(dccrg::Dccrg& mpiGrid, } } //spatial loops - // Copy swapGrid back to momentsGrid - momentsGrid=swapGrid; - // Update Ghost Cells + // Allows argument dependent lookup (ADL) + // i.e. use specialized swap if it exists, fall back on std + using std::swap; + swap(momentsGrid, swapGrid); momentsGrid.updateGhostCells(); - - } + } } void feedMomentsIntoFsGrid(dccrg::Dccrg& mpiGrid, @@ -233,6 +229,7 @@ void getFieldsFromFsGrid( FsGrid< std::array, FS_STENCIL_WIDTH> & volumeFieldsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, dccrg::Dccrg& mpiGrid, const std::vector& cells @@ -315,7 +312,9 @@ void getFieldsFromFsGrid( std::array * volcell = volumeFieldsGrid.get(fsgridCell); std::array * bgcell = BgBGrid.get(fsgridCell); std::array * egradpecell = EGradPeGrid.get(fsgridCell); + std::array * dMomentscell = dMomentsGrid.get(fsgridCell); + // TODO consider pruning these and communicating only when required sendBuffer[ii].sums[FieldsToCommunicate::PERBXVOL] += volcell->at(fsgrids::volfields::PERBXVOL); sendBuffer[ii].sums[FieldsToCommunicate::PERBYVOL] += volcell->at(fsgrids::volfields::PERBYVOL); sendBuffer[ii].sums[FieldsToCommunicate::PERBZVOL] += volcell->at(fsgrids::volfields::PERBZVOL); @@ -328,6 +327,15 @@ void getFieldsFromFsGrid( sendBuffer[ii].sums[FieldsToCommunicate::dPERBZVOLdx] += volcell->at(fsgrids::volfields::dPERBZVOLdx) / technicalGrid.DX; sendBuffer[ii].sums[FieldsToCommunicate::dPERBZVOLdy] += volcell->at(fsgrids::volfields::dPERBZVOLdy) / technicalGrid.DY; sendBuffer[ii].sums[FieldsToCommunicate::dPERBZVOLdz] += volcell->at(fsgrids::volfields::dPERBZVOLdz) / technicalGrid.DZ; + sendBuffer[ii].sums[FieldsToCommunicate::dVxdx] += dMomentscell->at(fsgrids::dmoments::dVxdx) / technicalGrid.DX; + sendBuffer[ii].sums[FieldsToCommunicate::dVxdy] += dMomentscell->at(fsgrids::dmoments::dVxdy) / technicalGrid.DY; + sendBuffer[ii].sums[FieldsToCommunicate::dVxdz] += dMomentscell->at(fsgrids::dmoments::dVxdz) / technicalGrid.DZ; + sendBuffer[ii].sums[FieldsToCommunicate::dVydx] += dMomentscell->at(fsgrids::dmoments::dVydx) / technicalGrid.DX; + sendBuffer[ii].sums[FieldsToCommunicate::dVydy] += dMomentscell->at(fsgrids::dmoments::dVydy) / technicalGrid.DY; + sendBuffer[ii].sums[FieldsToCommunicate::dVydz] += dMomentscell->at(fsgrids::dmoments::dVydz) / technicalGrid.DZ; + sendBuffer[ii].sums[FieldsToCommunicate::dVzdx] += dMomentscell->at(fsgrids::dmoments::dVzdx) / technicalGrid.DX; + sendBuffer[ii].sums[FieldsToCommunicate::dVzdy] += dMomentscell->at(fsgrids::dmoments::dVzdy) / technicalGrid.DY; + sendBuffer[ii].sums[FieldsToCommunicate::dVzdz] += dMomentscell->at(fsgrids::dmoments::dVzdz) / technicalGrid.DZ; sendBuffer[ii].sums[FieldsToCommunicate::BGBXVOL] += bgcell->at(fsgrids::bgbfield::BGBXVOL); sendBuffer[ii].sums[FieldsToCommunicate::BGBYVOL] += bgcell->at(fsgrids::bgbfield::BGBYVOL); sendBuffer[ii].sums[FieldsToCommunicate::BGBZVOL] += bgcell->at(fsgrids::bgbfield::BGBZVOL); @@ -387,6 +395,15 @@ void getFieldsFromFsGrid( mpiGrid[cellAggregate.first]->derivativesBVOL[bvolderivatives::dPERBZVOLdx] = cellAggregate.second.sums[FieldsToCommunicate::dPERBZVOLdx] / cellAggregate.second.cells; mpiGrid[cellAggregate.first]->derivativesBVOL[bvolderivatives::dPERBZVOLdy] = cellAggregate.second.sums[FieldsToCommunicate::dPERBZVOLdy] / cellAggregate.second.cells; mpiGrid[cellAggregate.first]->derivativesBVOL[bvolderivatives::dPERBZVOLdz] = cellAggregate.second.sums[FieldsToCommunicate::dPERBZVOLdz] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVxdx] = cellAggregate.second.sums[FieldsToCommunicate::dVxdx] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVxdy] = cellAggregate.second.sums[FieldsToCommunicate::dVxdy] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVxdz] = cellAggregate.second.sums[FieldsToCommunicate::dVxdz] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVydx] = cellAggregate.second.sums[FieldsToCommunicate::dVydx] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVydy] = cellAggregate.second.sums[FieldsToCommunicate::dVydy] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVydz] = cellAggregate.second.sums[FieldsToCommunicate::dVydz] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVzdx] = cellAggregate.second.sums[FieldsToCommunicate::dVzdx] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVzdy] = cellAggregate.second.sums[FieldsToCommunicate::dVzdy] / cellAggregate.second.cells; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVzdz] = cellAggregate.second.sums[FieldsToCommunicate::dVzdz] / cellAggregate.second.cells; cellParams[CellParams::BGBXVOL] = cellAggregate.second.sums[FieldsToCommunicate::BGBXVOL] / cellAggregate.second.cells; cellParams[CellParams::BGBYVOL] = cellAggregate.second.sums[FieldsToCommunicate::BGBYVOL] / cellAggregate.second.cells; cellParams[CellParams::BGBZVOL] = cellAggregate.second.sums[FieldsToCommunicate::BGBZVOL] / cellAggregate.second.cells; @@ -414,6 +431,15 @@ void getFieldsFromFsGrid( mpiGrid[cellAggregate.first]->derivativesBVOL[bvolderivatives::dPERBZVOLdx] = 0; mpiGrid[cellAggregate.first]->derivativesBVOL[bvolderivatives::dPERBZVOLdy] = 0; mpiGrid[cellAggregate.first]->derivativesBVOL[bvolderivatives::dPERBZVOLdz] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVxdx] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVxdy] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVxdz] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVydx] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVydy] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVydz] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVzdx] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVzdy] = 0; + mpiGrid[cellAggregate.first]->derivativesV[vderivatives::dVzdz] = 0; cellParams[CellParams::BGBXVOL] = 0; cellParams[CellParams::BGBYVOL] = 0; cellParams[CellParams::BGBZVOL] = 0; diff --git a/fieldsolver/gridGlue.hpp b/fieldsolver/gridGlue.hpp index 43a03d417..f58d40dbc 100644 --- a/fieldsolver/gridGlue.hpp +++ b/fieldsolver/gridGlue.hpp @@ -36,6 +36,15 @@ enum FieldsToCommunicate { CURVATUREX, CURVATUREY, CURVATUREZ, + dVxdx, + dVxdy, + dVxdz, + dVydx, + dVydy, + dVydz, + dVzdx, + dVzdy, + dVzdz, N_FIELDSTOCOMMUNICATE }; @@ -64,12 +73,14 @@ void feedMomentsIntoFsGrid(dccrg::Dccrg& * * This function assumes that proper grid coupling has been set up. */ -void getFieldsFromFsGrid(FsGrid< std::array, FS_STENCIL_WIDTH> & volumeFieldsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, - FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, - dccrg::Dccrg& mpiGrid, - const std::vector& cells +void getFieldsFromFsGrid( + FsGrid< std::array, FS_STENCIL_WIDTH> & volumeFieldsGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, + FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, + dccrg::Dccrg& mpiGrid, + const std::vector& cells ); /*! Copy background B fields and store them into DCCRG diff --git a/fieldsolver/ldz_electric_field.cpp b/fieldsolver/ldz_electric_field.cpp index 7d9b52387..2453a50c2 100644 --- a/fieldsolver/ldz_electric_field.cpp +++ b/fieldsolver/ldz_electric_field.cpp @@ -1634,7 +1634,8 @@ void calculateElectricField( * \param momentsGrid fsGrid holding the moment quantities at runge-kutta t=0 * \param momentsDt2Grid fsGrid holding the moment quantities at runge-kutta t=0.5 * \param dPerBGrid fsGrid holding the derivatives of perturbed B - * \param dMomentsGrid fsGrid holding the derviatives of moments + * \param dMomentsGrid fsGrid holding the derivatives of moments + * \param dMomentsDt2Grid fsGrid holding the derivatives of moments at runge-kutta t=0.5 * \param BgBGrid fsGrid holding the background B quantities * \param technicalGrid fsGrid holding technical information (such as boundary types) * \param sysBoundaries System boundary conditions existing @@ -1649,14 +1650,17 @@ void calculateUpwindedElectricFieldSimple( FsGrid< std::array, FS_STENCIL_WIDTH> & EDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & EHallGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, SysBoundary& sysBoundaries, - cint& RKCase + cint& RKCase, + const bool communicateEGradPeOrMomentsDerivatives ) { //const std::array gridDims = technicalGrid.getLocalSize(); const FsGridTools::FsIndex_t* gridDims = &technicalGrid.getLocalSize()[0]; @@ -1669,14 +1673,22 @@ void calculateUpwindedElectricFieldSimple( if(P::ohmHallTerm > 0) { EHallGrid.updateGhostCells(); } - if(P::ohmGradPeTerm > 0) { - EGradPeGrid.updateGhostCells(); + if(P::ohmGradPeTerm > 0 && communicateEGradPeOrMomentsDerivatives) { + if (RKCase == RK_ORDER1 || RKCase == RK_ORDER2_STEP2) { + EGradPeGrid.updateGhostCells(); + } else { + EGradPeDt2Grid.updateGhostCells(); + } } if(P::ohmHallTerm == 0) { dPerBGrid.updateGhostCells(); } - if(P::ohmHallTerm == 0 && P::ohmGradPeTerm == 0) { - dMomentsGrid.updateGhostCells(); + if(P::ohmHallTerm == 0 && P::ohmGradPeTerm == 0 && communicateEGradPeOrMomentsDerivatives) { + if (RKCase == RK_ORDER1 || RKCase == RK_ORDER2_STEP2) { + dMomentsGrid.updateGhostCells(); + } else { + dMomentsDt2Grid.updateGhostCells(); + } } mpiTimer.stop(); @@ -1705,16 +1717,16 @@ void calculateUpwindedElectricFieldSimple( k, sysBoundaries, RKCase - ); + ); } else { // RKCase == RK_ORDER2_STEP1 calculateElectricField( perBDt2Grid, EDt2Grid, EHallGrid, - EGradPeGrid, + EGradPeDt2Grid, momentsDt2Grid, dPerBGrid, - dMomentsGrid, + dMomentsDt2Grid, BgBGrid, technicalGrid, i, @@ -1722,7 +1734,7 @@ void calculateUpwindedElectricFieldSimple( k, sysBoundaries, RKCase - ); + ); } } } diff --git a/fieldsolver/ldz_electric_field.hpp b/fieldsolver/ldz_electric_field.hpp index f286c5965..3cb4be205 100644 --- a/fieldsolver/ldz_electric_field.hpp +++ b/fieldsolver/ldz_electric_field.hpp @@ -29,14 +29,17 @@ void calculateUpwindedElectricFieldSimple( FsGrid< std::array, FS_STENCIL_WIDTH> & EDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & EHallGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, SysBoundary& sysBoundaries, - cint& RKCase + cint& RKCase, + const bool communicateEGradPeOrMomentsDerivatives ); #endif diff --git a/fieldsolver/ldz_gradpe.cpp b/fieldsolver/ldz_gradpe.cpp index 0f3d35f4a..8044de448 100644 --- a/fieldsolver/ldz_gradpe.cpp +++ b/fieldsolver/ldz_gradpe.cpp @@ -152,9 +152,11 @@ void calculateGradPeTerm( void calculateGradPeTermSimple( FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsDt2Grid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, SysBoundary& sysBoundaries, cint& RKCase @@ -166,7 +168,11 @@ void calculateGradPeTermSimple( int computeTimerId {phiprof::initializeTimer("EgradPe compute cells")}; phiprof::Timer mpiTimer {"EgradPe field update ghosts MPI", {"MPI"}}; - dMomentsGrid.updateGhostCells(); + if (RKCase == RK_ORDER1 || RKCase == RK_ORDER2_STEP2) { + dMomentsGrid.updateGhostCells(); + } else { + dMomentsDt2Grid.updateGhostCells(); + } mpiTimer.stop(); // Calculate GradPe term @@ -180,7 +186,7 @@ void calculateGradPeTermSimple( if (RKCase == RK_ORDER1 || RKCase == RK_ORDER2_STEP2) { calculateGradPeTerm(EGradPeGrid, momentsGrid, dMomentsGrid, technicalGrid, i, j, k, sysBoundaries); } else { - calculateGradPeTerm(EGradPeGrid, momentsDt2Grid, dMomentsGrid, technicalGrid, i, j, k, sysBoundaries); + calculateGradPeTerm(EGradPeDt2Grid, momentsDt2Grid, dMomentsDt2Grid, technicalGrid, i, j, k, sysBoundaries); } } } diff --git a/fieldsolver/ldz_gradpe.hpp b/fieldsolver/ldz_gradpe.hpp index 23caaf2f1..9e76b5cc7 100644 --- a/fieldsolver/ldz_gradpe.hpp +++ b/fieldsolver/ldz_gradpe.hpp @@ -24,9 +24,11 @@ void calculateGradPeTermSimple( FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsDt2Grid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, SysBoundary& sysBoundaries, cint& RKCase diff --git a/fieldsolver/ldz_hall.cpp b/fieldsolver/ldz_hall.cpp index b403c4a18..6b9357cbc 100644 --- a/fieldsolver/ldz_hall.cpp +++ b/fieldsolver/ldz_hall.cpp @@ -423,7 +423,7 @@ REAL JXBZ_110_111( * \param EHallGrid fsGrid holding the Hall contributions to the electric field * \param momentsGrid fsGrid holding the moment quantities * \param dPerBGrid fsGrid holding the derivatives of perturbed B - * \param dMomentsGrid fsGrid holding the derviatives of moments + * \param dMomentsGrid fsGrid holding the derivatives of moments * \param BgBGrid fsGrid holding the background B quantities * \param technicalGrid fsGrid holding technical information (such as boundary types) * \param perturbedCoefficients Reconstruction coefficients @@ -521,7 +521,7 @@ void calculateEdgeHallTermXComponents( * \param EHallGrid fsGrid holding the Hall contributions to the electric field * \param momentsGrid fsGrid holding the moment quantities * \param dPerBGrid fsGrid holding the derivatives of perturbed B - * \param dMomentsGrid fsGrid holding the derviatives of moments + * \param dMomentsGrid fsGrid holding the derivatives of moments * \param BgBGrid fsGrid holding the background B quantities * \param technicalGrid fsGrid holding technical information (such as boundary types) * \param perturbedCoefficients Reconstruction coefficients @@ -619,7 +619,7 @@ void calculateEdgeHallTermYComponents( * \param EHallGrid fsGrid holding the Hall contributions to the electric field * \param momentsGrid fsGrid holding the moment quantities * \param dPerBGrid fsGrid holding the derivatives of perturbed B - * \param dMomentsGrid fsGrid holding the derviatives of moments + * \param dMomentsGrid fsGrid holding the derivatives of moments * \param BgBGrid fsGrid holding the background B quantities * \param technicalGrid fsGrid holding technical information (such as boundary types) * \param perturbedCoefficients Reconstruction coefficients @@ -715,7 +715,7 @@ void calculateEdgeHallTermZComponents( * \param EHallGrid fsGrid holding the Hall contributions to the electric field * \param momentsGrid fsGrid holding the moment quantities * \param dPerBGrid fsGrid holding the derivatives of perturbed B - * \param dMomentsGrid fsGrid holding the derviatives of moments + * \param dMomentsGrid fsGrid holding the derivatives of moments * \param BgBGrid fsGrid holding the background B quantities * \param technicalGrid fsGrid holding technical information (such as boundary types) * \param sysBoundaries System boundary condition functions. @@ -784,12 +784,13 @@ void calculateHallTerm( * \param momentsGrid fsGrid holding the moment quantities * \param momentsDt2Grid fsGrid holding the moment quantities at runge-kutta half step * \param dPerBGrid fsGrid holding the derivatives of perturbed B - * \param dMomentsGrid fsGrid holding the derviatives of moments + * \param dMomentsGrid fsGrid holding the derivatives of moments + * \param dMomentsDt2Grid fsGrid holding the derivatives of moments at runge-kutta half step * \param BgBGrid fsGrid holding the background B quantities * \param technicalGrid fsGrid holding technical information (such as boundary types) * \param sysBoundaries System boundary condition functions. * \param RKCase Element in the enum defining the Runge-Kutta method steps - * \param communicateMomentsDerivatives whether to communicate derivatves with the neighbour CPUs + * \param communicateMomentsDerivatives whether to communicate derivatives with the neighbour CPUs * * \sa calculateHallTerm */ @@ -801,10 +802,12 @@ void calculateHallTermSimple( FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, SysBoundary& sysBoundaries, - cint& RKCase + cint& RKCase, + const bool communicateMomentsDerivatives ) { //const std::array gridDims = technicalGrid.getLocalSize(); const FsGridTools::FsIndex_t* gridDims = &technicalGrid.getLocalSize()[0]; @@ -814,8 +817,12 @@ void calculateHallTermSimple( phiprof::Timer mpiTimer {"EHall ghost updates MPI", {"MPI"}}; int computeTimerId {phiprof::initializeTimer("EHall compute cells")}; dPerBGrid.updateGhostCells(); - if(P::ohmGradPeTerm == 0) { - dMomentsGrid.updateGhostCells(); + if(P::ohmGradPeTerm == 0 && communicateMomentsDerivatives) { + if (RKCase == RK_ORDER1 || RKCase == RK_ORDER2_STEP2) { + dMomentsGrid.updateGhostCells(); + } else { + dMomentsDt2Grid.updateGhostCells(); + } } mpiTimer.stop(); @@ -829,7 +836,7 @@ void calculateHallTermSimple( if (RKCase == RK_ORDER1 || RKCase == RK_ORDER2_STEP2) { calculateHallTerm(perBGrid, EHallGrid, momentsGrid, dPerBGrid, dMomentsGrid, BgBGrid, technicalGrid,sysBoundaries, i, j, k); } else { - calculateHallTerm(perBDt2Grid, EHallGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, BgBGrid, technicalGrid,sysBoundaries, i, j, k); + calculateHallTerm(perBDt2Grid, EHallGrid, momentsDt2Grid, dPerBGrid, dMomentsDt2Grid, BgBGrid, technicalGrid,sysBoundaries, i, j, k); } } } diff --git a/fieldsolver/ldz_hall.hpp b/fieldsolver/ldz_hall.hpp index 2b5aa3545..b9d1e156b 100644 --- a/fieldsolver/ldz_hall.hpp +++ b/fieldsolver/ldz_hall.hpp @@ -30,10 +30,12 @@ void calculateHallTermSimple( FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, SysBoundary& sysBoundaries, - cint& RKCase + cint& RKCase, + const bool communicateMomentsDerivatives ); #endif diff --git a/fieldsolver/ldz_main.cpp b/fieldsolver/ldz_main.cpp index 2f1d769e4..a85897b41 100644 --- a/fieldsolver/ldz_main.cpp +++ b/fieldsolver/ldz_main.cpp @@ -54,22 +54,6 @@ #include "mpiconversion.h" #include "../fieldtracing/fieldtracing.h" -/*! Re-initialize field propagator after rebalance. E, BGB, RHO, RHO_V, - cell_dimensions, sysboundaryflag need to be up to date for the - extended neighborhood - */ -bool initializeFieldPropagatorAfterRebalance() { - // Assume static background field, they are not communicated here - // but are assumed to be ok after each load balance as that - // communicates all spatial data - - return true; -} - -bool finalizeFieldPropagator() { - return true; -} - /*! \brief Top-level field propagation function. * * Propagates the magnetic field, computes the derivatives and the upwinded @@ -91,10 +75,12 @@ bool propagateFields( FsGrid< std::array, FS_STENCIL_WIDTH> & EDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & EHallGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsDt2Grid, FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & volGrid, FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, @@ -123,9 +109,9 @@ bool propagateFields( if (subcycles == 1) { #ifdef FS_1ST_ORDER_TIME propagateMagneticFieldSimple(perBGrid, perBDt2Grid, BgBGrid, EGrid, EDt2Grid, technicalGrid, sysBoundaries, dt, RK_ORDER1); - calculateDerivativesSimple(perBGrid, perBDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, technicalGrid, sysBoundaries, RK_ORDER1, true); + calculateDerivativesSimple(perBGrid, perBDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, dMomentsDt2Grid, technicalGrid, sysBoundaries, RK_ORDER1, true/*doMoments*/); if(P::ohmGradPeTerm > 0){ - calculateGradPeTermSimple(EGradPeGrid, momentsGrid, momentsDt2Grid, dMomentsGrid, technicalGrid, sysBoundaries, RK_ORDER1); + calculateGradPeTermSimple(EGradPeGrid, EGradPeDt2Grid, momentsGrid, momentsDt2Grid, dMomentsGrid, dMomentsDt2Grid, technicalGrid, sysBoundaries, RK_ORDER1); } if(P::ohmHallTerm > 0) { calculateHallTermSimple( @@ -136,10 +122,12 @@ bool propagateFields( momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, BgBGrid, technicalGrid, sysBoundaries, - RK_ORDER1 + RK_ORDER1, + true // communicateMomentsDerivatives ); } calculateUpwindedElectricFieldSimple( @@ -149,20 +137,23 @@ bool propagateFields( EDt2Grid, EHallGrid, EGradPeGrid, + EGradPeDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, BgBGrid, technicalGrid, sysBoundaries, - RK_ORDER1 + RK_ORDER1, + true // communicateEGradPeOrMomentsDerivatives ); #else propagateMagneticFieldSimple(perBGrid, perBDt2Grid, BgBGrid, EGrid, EDt2Grid, technicalGrid, sysBoundaries, dt, RK_ORDER2_STEP1); - calculateDerivativesSimple(perBGrid, perBDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, technicalGrid, sysBoundaries, RK_ORDER2_STEP1, true); + calculateDerivativesSimple(perBGrid, perBDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, dMomentsDt2Grid, technicalGrid, sysBoundaries, RK_ORDER2_STEP1, true/*doMoments*/); if(P::ohmGradPeTerm > 0) { - calculateGradPeTermSimple(EGradPeGrid, momentsGrid, momentsDt2Grid, dMomentsGrid, technicalGrid, sysBoundaries, RK_ORDER2_STEP1); + calculateGradPeTermSimple(EGradPeGrid, EGradPeDt2Grid, momentsGrid, momentsDt2Grid, dMomentsGrid, dMomentsDt2Grid, technicalGrid, sysBoundaries, RK_ORDER2_STEP1); } if(P::ohmHallTerm > 0) { calculateHallTermSimple( @@ -173,10 +164,12 @@ bool propagateFields( momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, BgBGrid, technicalGrid, sysBoundaries, - RK_ORDER2_STEP1 + RK_ORDER2_STEP1, + true // communicateMomentsDerivatives ); } calculateUpwindedElectricFieldSimple( @@ -186,20 +179,23 @@ bool propagateFields( EDt2Grid, EHallGrid, EGradPeGrid, + EGradPeDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, BgBGrid, technicalGrid, sysBoundaries, - RK_ORDER2_STEP1 + RK_ORDER2_STEP1, + true // communicateEGradPeOrMomentsDerivatives ); propagateMagneticFieldSimple(perBGrid, perBDt2Grid, BgBGrid, EGrid, EDt2Grid, technicalGrid, sysBoundaries, dt, RK_ORDER2_STEP2); - calculateDerivativesSimple(perBGrid, perBDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, technicalGrid, sysBoundaries, RK_ORDER2_STEP2, true); + calculateDerivativesSimple(perBGrid, perBDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, dMomentsDt2Grid, technicalGrid, sysBoundaries, RK_ORDER2_STEP2, true/*doMoments*/); if(P::ohmGradPeTerm > 0) { - calculateGradPeTermSimple(EGradPeGrid, momentsGrid, momentsDt2Grid, dMomentsGrid, technicalGrid, sysBoundaries, RK_ORDER2_STEP2); + calculateGradPeTermSimple(EGradPeGrid, EGradPeDt2Grid, momentsGrid, momentsDt2Grid, dMomentsGrid, dMomentsDt2Grid, technicalGrid, sysBoundaries, RK_ORDER2_STEP2); } if(P::ohmHallTerm > 0) { calculateHallTermSimple( @@ -210,10 +206,12 @@ bool propagateFields( momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, BgBGrid, technicalGrid, sysBoundaries, - RK_ORDER2_STEP2 + RK_ORDER2_STEP2, + true // communicateMomentsDerivatives ); } calculateUpwindedElectricFieldSimple( @@ -223,14 +221,17 @@ bool propagateFields( EDt2Grid, EHallGrid, EGradPeGrid, + EGradPeDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, BgBGrid, technicalGrid, sysBoundaries, - RK_ORDER2_STEP2 + RK_ORDER2_STEP2, + true // communicateEGradPeOrMomentsDerivatives ); #endif } else { @@ -248,9 +249,9 @@ bool propagateFields( // We need to calculate derivatives of the moments at every substep, but the moments only // need to be communicated in the first one. - calculateDerivativesSimple(perBGrid, perBDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, technicalGrid, sysBoundaries, RK_ORDER2_STEP1, (subcycleCount==0)); + calculateDerivativesSimple(perBGrid, perBDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, dMomentsDt2Grid, technicalGrid, sysBoundaries, RK_ORDER2_STEP1, (subcycleCount==0)/*doMoments*/); if(P::ohmGradPeTerm > 0 && subcycleCount==0) { - calculateGradPeTermSimple(EGradPeGrid, momentsGrid, momentsDt2Grid, dMomentsGrid, technicalGrid, sysBoundaries, RK_ORDER2_STEP1); + calculateGradPeTermSimple(EGradPeGrid, EGradPeDt2Grid, momentsGrid, momentsDt2Grid, dMomentsGrid, dMomentsDt2Grid, technicalGrid, sysBoundaries, RK_ORDER2_STEP1); } if(P::ohmHallTerm > 0) { calculateHallTermSimple( @@ -261,10 +262,12 @@ bool propagateFields( momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, BgBGrid, technicalGrid, sysBoundaries, - RK_ORDER2_STEP1 + RK_ORDER2_STEP1, + subcycleCount==0 // communicateMomentsDerivatives ); } calculateUpwindedElectricFieldSimple( @@ -274,23 +277,26 @@ bool propagateFields( EDt2Grid, EHallGrid, EGradPeGrid, + EGradPeDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, BgBGrid, technicalGrid, sysBoundaries, - RK_ORDER2_STEP1 + RK_ORDER2_STEP1, + subcycleCount==0 // communicateEGradPeOrMomentsDerivatives ); propagateMagneticFieldSimple(perBGrid, perBDt2Grid, BgBGrid, EGrid, EDt2Grid, technicalGrid, sysBoundaries, subcycleDt, RK_ORDER2_STEP2); // We need to calculate derivatives of the moments at every substep, but the moments only // need to be communicated in the first one. - calculateDerivativesSimple(perBGrid, perBDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, technicalGrid, sysBoundaries, RK_ORDER2_STEP2, (subcycleCount==0)); + calculateDerivativesSimple(perBGrid, perBDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, dMomentsDt2Grid, technicalGrid, sysBoundaries, RK_ORDER2_STEP2, (subcycleCount==0)/*doMoments*/); if(P::ohmGradPeTerm > 0 && subcycleCount==0) { - calculateGradPeTermSimple(EGradPeGrid, momentsGrid, momentsDt2Grid, dMomentsGrid, technicalGrid, sysBoundaries, RK_ORDER2_STEP2); + calculateGradPeTermSimple(EGradPeGrid, EGradPeDt2Grid, momentsGrid, momentsDt2Grid, dMomentsGrid, dMomentsDt2Grid, technicalGrid, sysBoundaries, RK_ORDER2_STEP2); } if(P::ohmHallTerm > 0) { calculateHallTermSimple( @@ -301,10 +307,12 @@ bool propagateFields( momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, BgBGrid, technicalGrid, sysBoundaries, - RK_ORDER2_STEP2 + RK_ORDER2_STEP2, + subcycleCount==0 // communicateMomentsDerivatives ); } calculateUpwindedElectricFieldSimple( @@ -314,14 +322,17 @@ bool propagateFields( EDt2Grid, EHallGrid, EGradPeGrid, + EGradPeDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, BgBGrid, technicalGrid, sysBoundaries, - RK_ORDER2_STEP2 + RK_ORDER2_STEP2, + subcycleCount==0 // communicateEGradPeOrMomentsDerivatives ); phiprof::Timer subcyclingTimer {"FS subcycle stuff"}; diff --git a/fieldtracing/fieldtracing.cpp b/fieldtracing/fieldtracing.cpp index 78c70d349..b57c451e9 100644 --- a/fieldtracing/fieldtracing.cpp +++ b/fieldtracing/fieldtracing.cpp @@ -356,7 +356,7 @@ namespace FieldTracing { const std::array x_in = x; Real r_in = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]); // Take a step back and find the innerRadius crossing point - stepFieldLine(x,v, stepSize,fieldTracingParameters.min_tracer_dx_ionospere_coupling,fieldTracingParameters.max_tracer_dx_ionospere_coupling,fieldTracingParameters.tracingMethod,dipoleFieldOnly,false); + stepFieldLine(x,v, stepSize,fieldTracingParameters.min_tracer_dx_ionospere_coupling,fieldTracingParameters.max_tracer_dx_ionospere_coupling,fieldTracingParameters.tracingMethod,dipoleFieldOnly,true); Real r_out = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]); Real alpha = (SBC::Ionosphere::innerRadius - r_in)/(r_out - r_in); alpha = std::fmax(std::fmin(alpha,1.0),0.0); diff --git a/grid.cpp b/grid.cpp index 12b543256..eb4e2531c 100644 --- a/grid.cpp +++ b/grid.cpp @@ -95,6 +95,7 @@ void initializeGrids( FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & volGrid, @@ -356,7 +357,7 @@ void initializeGrids( volGrid.updateGhostCells(); fsGridGhostTimer.stop(); phiprof::Timer getFieldsTimer {"getFieldsFromFsGrid"}; - getFieldsFromFsGrid(volGrid, BgBGrid, EGradPeGrid, technicalGrid, mpiGrid, cells); + getFieldsFromFsGrid(volGrid, BgBGrid, EGradPeGrid, dMomentsGrid, technicalGrid, mpiGrid, cells); getFieldsTimer.stop(); setBTimer.stop(); @@ -682,16 +683,6 @@ void balanceLoad(dccrg::Dccrg& mpiGrid, S phiprof::Timer updateBoundariesTimer {"update sysboundaries"}; sysBoundaries.updateSysBoundariesAfterLoadBalance( mpiGrid ); updateBoundariesTimer.stop(); - - phiprof::Timer initSolversTimer {"Init solvers"}; - // Initialize field propagator (only if in use): - if (Parameters::propagateField == true) { - if (initializeFieldPropagatorAfterRebalance() == false) { - logFile << "(MAIN): Field propagator did not initialize correctly!" << endl << writeVerbose; - exit(1); - } - } - initSolversTimer.stop(); // Record ranks of face neighbors if(P::amrMaxSpatialRefLevel > 0) { diff --git a/grid.h b/grid.h index 0b4c6628a..7f389a819 100644 --- a/grid.h +++ b/grid.h @@ -41,6 +41,7 @@ void initializeGrids( FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, + FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, FsGrid< std::array, FS_STENCIL_WIDTH> & volGrid, diff --git a/ioread.cpp b/ioread.cpp index ece9cb2cd..cad300904 100644 --- a/ioread.cpp +++ b/ioread.cpp @@ -1296,7 +1296,6 @@ bool exec_readGrid(dccrg::Dccrg& mpiGrid, if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"max_r_dt",CellParams::MAXRDT,1,mpiGrid); } if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"max_fields_dt",CellParams::MAXFDT,1,mpiGrid); } if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"vg_drift",CellParams::BULKV_FORCING_X,3,mpiGrid); } - if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"vg_bulk_forcing_flag",CellParams::FORCING_CELL_NUM,1,mpiGrid); } if (P::refineOnRestart) { // Refinement indices alpha_1 and alpha_2 if(success) { success=readCellParamsVariable(file,fileCells,localCellStartOffset,localCells,"vg_amr_alpha1",CellParams::AMR_ALPHA1,1,mpiGrid); } diff --git a/iowrite.cpp b/iowrite.cpp index 219641907..d609149fd 100644 --- a/iowrite.cpp +++ b/iowrite.cpp @@ -1661,7 +1661,6 @@ bool writeRestart( restartReducer.addOperator(new DRO::DataReductionOperatorCellParams("max_r_dt",CellParams::MAXRDT,1)); restartReducer.addOperator(new DRO::DataReductionOperatorCellParams("max_fields_dt",CellParams::MAXFDT,1)); restartReducer.addOperator(new DRO::DataReductionOperatorCellParams("vg_drift",CellParams::BULKV_FORCING_X,3)); - restartReducer.addOperator(new DRO::DataReductionOperatorCellParams("vg_bulk_forcing_flag",CellParams::FORCING_CELL_NUM,1)); restartReducer.addOperator(new DRO::DataReductionOperatorCellParams("vg_amr_alpha1",CellParams::AMR_ALPHA1,1)); restartReducer.addOperator(new DRO::DataReductionOperatorCellParams("vg_amr_alpha2",CellParams::AMR_ALPHA2,1)); restartReducer.addOperator(new DRO::VariableBVol); diff --git a/parameters.cpp b/parameters.cpp index b6b5eb4bb..f360fde4a 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -167,6 +167,13 @@ Real P::alpha1CoarsenThreshold = -1.0; bool P::useAlpha2 = true; Real P::alpha2RefineThreshold = 0.5; Real P::alpha2CoarsenThreshold = -1.0; +bool P::useVorticity = false; +Real P::vorticityRefineThreshold = 0.5; +Real P::vorticityCoarsenThreshold = -1.0; +bool P::useAnisotropy = false; +Real P::anisotropyRefineThreshold = 0.5; +Real P::anisotropyCoarsenThreshold = -1.0; +int P::anisotropyMaxReflevel = 2; Real P::alphaDRhoWeight = 1.0; Real P::alphaDUWeight = 1.0; Real P::alphaDPSqWeight = 1.0; @@ -409,7 +416,7 @@ bool P::addParameters() { "ig_precipitation ig_deltaphi "+ "ig_inplanecurrent ig_b ig_e vg_drift vg_ionospherecoupling vg_connection vg_fluxrope fg_curvature "+ "vg_amr_drho vg_amr_du vg_amr_dpsq vg_amr_dbsq vg_amr_db vg_amr_alpha1 vg_amr_reflevel vg_amr_alpha2 "+ - "vg_amr_translate_comm vg_gridcoordinates fg_gridcoordinates "); + "vg_amr_translate_comm vg_gridcoordinates fg_gridcoordinates vg_pressure_anisotropy vg_amr_vorticity"); RP::addComposing( "variables_deprecated.output", @@ -484,6 +491,13 @@ bool P::addParameters() { RP::add("AMR.use_alpha2","Use J/B_perp as a refinement index", true); RP::add("AMR.alpha2_refine_threshold","Determines the minimum value of alpha_2 to refine cells", 0.5); RP::add("AMR.alpha2_coarsen_threshold","Determines the maximum value of alpha_2 to unrefine cells, default half of the refine threshold", -1.0); + RP::add("AMR.use_vorticity","Use vorticity as a refinement index", false); + RP::add("AMR.vorticity_refine_threshold","Determines the minimum value of vorticity to refine cells", 0.5); + RP::add("AMR.vorticity_coarsen_threshold","Determines the maximum value of vorticity to unrefine cells, default half of the refine threshold", -1.0); + RP::add("AMR.use_anisotropy","Use pressure anisotropy as a refinement index", false); + RP::add("AMR.anisotropy_refine_threshold","Determines the maximum value of pressure anisotropy to refine cells", 0.5); + RP::add("AMR.anisotropy_coarsen_threshold","Determines the minimum value of pressure anisotropy to unrefine cells, default twice the refine threshold", -1.0); + RP::add("AMR.anisotropy_max_reflevel","When anisotropy is below the refine threshold, defines the maximum level to refine to", 2); RP::add("AMR.refine_cadence","Refine every nth load balance", 5); RP::add("AMR.refine_after","Start refinement after this many simulation seconds", 0.0); RP::add("AMR.refine_radius","Maximum distance from origin to allow refinement within. Only induced refinement allowed outside this radius.", LARGE_REAL); @@ -788,6 +802,7 @@ void Parameters::getParameters() { } MPI_Abort(MPI_COMM_WORLD, 1); } + RP::get("AMR.use_alpha2",P::useAlpha2); RP::get("AMR.alpha2_refine_threshold",P::alpha2RefineThreshold); RP::get("AMR.alpha2_coarsen_threshold",P::alpha2CoarsenThreshold); @@ -801,6 +816,33 @@ void Parameters::getParameters() { MPI_Abort(MPI_COMM_WORLD, 1); } + RP::get("AMR.use_vorticity",P::useVorticity); + RP::get("AMR.vorticity_refine_threshold",P::vorticityRefineThreshold); + RP::get("AMR.vorticity_coarsen_threshold",P::vorticityCoarsenThreshold); + if (P::useVorticity && P::vorticityCoarsenThreshold < 0) { + P::vorticityCoarsenThreshold = P::vorticityRefineThreshold / 2.0; + } + if (P::useVorticity && P::vorticityRefineThreshold < 0) { + if (myRank == MASTER_RANK) { + cerr << "ERROR invalid vorticity refine threshold" << endl; + } + MPI_Abort(MPI_COMM_WORLD, 1); + } + + RP::get("AMR.use_anisotropy",P::useAnisotropy); + RP::get("AMR.anisotropy_refine_threshold", P::anisotropyRefineThreshold); + RP::get("AMR.anisotropy_coarsen_threshold", P::anisotropyCoarsenThreshold); + RP::get("AMR.anisotropy_max_reflevel", P::anisotropyMaxReflevel); + if (P::useAnisotropy && P::anisotropyCoarsenThreshold < 0) { + P::anisotropyCoarsenThreshold = P::anisotropyRefineThreshold * 2.0; + } + if (P::useAnisotropy && P::anisotropyRefineThreshold < 0) { + if (myRank == MASTER_RANK) { + cerr << "ERROR invalid anisotropy refine threshold" << endl; + } + MPI_Abort(MPI_COMM_WORLD, 1); + } + RP::get("AMR.refine_cadence",P::refineCadence); RP::get("AMR.refine_after",P::refineAfter); RP::get("AMR.refine_radius",P::refineRadius); diff --git a/parameters.h b/parameters.h index 87838ef9f..bd63bdbad 100644 --- a/parameters.h +++ b/parameters.h @@ -198,6 +198,14 @@ struct Parameters { static bool useAlpha2; static Real alpha2RefineThreshold; static Real alpha2CoarsenThreshold; + // TODO: consider renaming to alpha3 or something to that effect + static bool useVorticity; + static Real vorticityRefineThreshold; + static Real vorticityCoarsenThreshold; + static bool useAnisotropy; + static Real anisotropyRefineThreshold; + static Real anisotropyCoarsenThreshold; + static int anisotropyMaxReflevel; static uint refineCadence; static Real refineAfter; static Real refineRadius; diff --git a/projects/project.cpp b/projects/project.cpp index 11d21c052..29151d2db 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -563,13 +563,88 @@ namespace projects { return cell->sysBoundaryFlag == sysboundarytype::NOT_SYSBOUNDARY && (cell->sysBoundaryLayer == 0 || cell->sysBoundaryLayer > 2); } + bool Project::shouldRefineCell(dccrg::Dccrg& mpiGrid, CellID id, Real r_max2) const { + // Evaluate possible refinement for this cell + + // Cells too far from the ionosphere should not be refined but + // induced refinement still possible just beyond this r_max2 limit. + + std::array xyz {mpiGrid.get_center(id)}; + SpatialCell* cell {mpiGrid[id]}; + int refLevel {mpiGrid.get_refinement_level(id)}; + Real r2 {pow(xyz[0], 2) + pow(xyz[1], 2) + pow(xyz[2], 2)}; + + bool alpha1ShouldRefine = (P::useAlpha1 && cell->parameters[CellParams::AMR_ALPHA1] > P::alpha1RefineThreshold); + bool alpha2ShouldRefine = (P::useAlpha2 && cell->parameters[CellParams::AMR_ALPHA2] > P::alpha2RefineThreshold); + bool vorticityShouldRefine = (P::useVorticity && cell->parameters[CellParams::AMR_VORTICITY] > P::vorticityRefineThreshold); + bool anisotropyShouldRefine = (P::useAnisotropy && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && refLevel < P::anisotropyMaxReflevel); + + bool shouldRefine { + (r2 < r_max2) && ( + alpha1ShouldRefine || + alpha2ShouldRefine || + vorticityShouldRefine || + anisotropyShouldRefine + ) + }; + + if( + // If this cell is planned to be refined, but is outside the allowed refinement region, cancel that refinement. + // Induced refinement still possible just beyond that limit. + (xyz[0] < P::refinementMinX) || (xyz[0] > P::refinementMaxX) + || (xyz[1] < P::refinementMinY) || (xyz[1] > P::refinementMaxY) + || (xyz[2] < P::refinementMinZ) || (xyz[2] > P::refinementMaxZ)) { + shouldRefine = false; + } + + return shouldRefine; + } + + bool Project::shouldUnrefineCell(dccrg::Dccrg& mpiGrid, CellID id, Real r_max2) const { + // Evaluate possible unrefinement for this cell + + // Cells too far from the ionosphere should be unrefined but + // induced refinement still possible just beyond this r_max2 limit. + + std::array xyz {mpiGrid.get_center(id)}; + SpatialCell* cell {mpiGrid[id]}; + int refLevel {mpiGrid.get_refinement_level(id)}; + Real r2 {pow(xyz[0], 2) + pow(xyz[1], 2) + pow(xyz[2], 2)}; + + bool alpha1ShouldUnrefine = (!P::useAlpha1 || cell->parameters[CellParams::AMR_ALPHA1] < P::alpha1CoarsenThreshold); + bool alpha2ShouldUnrefine = (!P::useAlpha2 || cell->parameters[CellParams::AMR_ALPHA2] < P::alpha2CoarsenThreshold); + bool vorticityShouldUnrefine = (!P::useVorticity || cell->parameters[CellParams::AMR_VORTICITY] < P::vorticityCoarsenThreshold); + bool anisotropyShouldUnrefine = (!P::useAnisotropy || cell->parameters[CellParams::P_ANISOTROPY] > P::anisotropyCoarsenThreshold || refLevel > P::anisotropyMaxReflevel); + + bool shouldUnrefine { + (r2 > r_max2) || ( + alpha1ShouldUnrefine && + alpha2ShouldUnrefine && + vorticityShouldUnrefine && + anisotropyShouldUnrefine + ) + }; + + if( + // If this cell is planned to remain at the current refinement level, but is outside the allowed refinement region, + // attempt to unrefine it instead. (If it is already at the lowest refinement level, DCCRG should not go belly-up.) + // Induced refinement still possible just beyond that limit. + (xyz[0] < P::refinementMinX) || (xyz[0] > P::refinementMaxX) + || (xyz[1] < P::refinementMinY) || (xyz[1] > P::refinementMaxY) + || (xyz[2] < P::refinementMinZ) || (xyz[2] > P::refinementMaxZ)) { + shouldUnrefine = true; + } + + return shouldUnrefine; + } + uint64_t Project::adaptRefinement( dccrg::Dccrg& mpiGrid ) const { phiprof::Timer refinesTimer {"Set refines"}; int myRank; MPI_Comm_rank(MPI_COMM_WORLD,&myRank); uint64_t refines {0}; - if (!P::useAlpha1 && !P::useAlpha2) { + if (!P::useAlpha1 && !P::useAlpha2 && !P::useAnisotropy && !P::useVorticity) { if (myRank == MASTER_RANK) { std::cout << "WARNING All refinement indices disabled" << std::endl; } @@ -581,10 +656,7 @@ namespace projects { //#pragma omp parallel for for (CellID id : cells) { - std::array xyz {mpiGrid.get_center(id)}; - SpatialCell* cell {mpiGrid[id]}; int refLevel {mpiGrid.get_refinement_level(id)}; - Real r2 {pow(xyz[0], 2) + pow(xyz[1], 2) + pow(xyz[2], 2)}; if (!canRefine(mpiGrid[id])) { // Skip refining, touching boundaries during runtime breaks everything @@ -593,58 +665,18 @@ namespace projects { } else { // Evaluate possible refinement or unrefinement for this cell - // Cells too far from the ionosphere should be unrefined but - // induced refinement still possible just beyond this r_max2 limit. - bool shouldRefine {(r2 < r_max2) && ((P::useAlpha1 ? cell->parameters[CellParams::AMR_ALPHA1] > P::alpha1RefineThreshold : false) || (P::useAlpha2 ? cell->parameters[CellParams::AMR_ALPHA2] > P::alpha2RefineThreshold : false))}; - bool shouldUnrefine {(r2 > r_max2) || ((P::useAlpha1 ? cell->parameters[CellParams::AMR_ALPHA1] < P::alpha1CoarsenThreshold : true) && (P::useAlpha2 ? cell->parameters[CellParams::AMR_ALPHA2] < P::alpha2CoarsenThreshold : true))}; - - if(shouldRefine - // If this cell is planned to be refined, but is outside the allowed refinement region, cancel that refinement. - // Induced refinement still possible just beyond that limit. - && ((xyz[0] < P::refinementMinX) || (xyz[0] > P::refinementMaxX) - || (xyz[1] < P::refinementMinY) || (xyz[1] > P::refinementMaxY) - || (xyz[2] < P::refinementMinZ) || (xyz[2] > P::refinementMaxZ))) { - shouldRefine = false; - } - if(!shouldUnrefine - // If this cell is planned to remain at the current refinement level, but is outside the allowed refinement region, - // attempt to unrefine it instead. (If it is already at the lowest refinement level, DCCRG should not go belly-up.) - // Induced refinement still possible just beyond that limit. - && ((xyz[0] < P::refinementMinX) || (xyz[0] > P::refinementMaxX) - || (xyz[1] < P::refinementMinY) || (xyz[1] > P::refinementMaxY) - || (xyz[2] < P::refinementMinZ) || (xyz[2] > P::refinementMaxZ))) { - shouldUnrefine = true; - } + bool shouldRefine = shouldRefineCell(mpiGrid, id, r_max2); + bool shouldUnrefine = shouldUnrefineCell(mpiGrid, id, r_max2); // Finally, check neighbors int refined_neighbors {0}; int coarser_neighbors {0}; for (const auto& [neighbor, dir] : mpiGrid.get_face_neighbors_of(id)) { // Evaluate all face neighbors of the current cell - std::array neighborXyz {mpiGrid.get_center(neighbor)}; - Real neighborR2 {pow(neighborXyz[0], 2) + pow(neighborXyz[1], 2) + pow(neighborXyz[2], 2)}; - const int neighborRef = mpiGrid.get_refinement_level(neighbor); - // Induced refinement still possible just beyond the r_max2 limit. - bool shouldRefineNeighbor {(neighborR2 < r_max2) && ((P::useAlpha1 ? mpiGrid[neighbor]->parameters[CellParams::AMR_ALPHA1] > P::alpha1RefineThreshold : false) || (P::useAlpha2 ? mpiGrid[neighbor]->parameters[CellParams::AMR_ALPHA2] > P::alpha2RefineThreshold : false))}; - if(shouldRefineNeighbor && - // If the neighbor is planned to be refined, but is outside the allowed refinement region, cancel that refinement. - // Induced refinement still possible just beyond that limit. - ((neighborXyz[0] < P::refinementMinX) || (neighborXyz[0] > P::refinementMaxX) - || (neighborXyz[1] < P::refinementMinY) || (neighborXyz[1] > P::refinementMaxY) - || (neighborXyz[2] < P::refinementMinZ) || (neighborXyz[2] > P::refinementMaxZ))) { - shouldRefineNeighbor = false; - } - // Induced refinement still possible just beyond the r_max2 limit. - bool shouldUnrefineNeighbor {(neighborR2 > r_max2) || ((P::useAlpha1 ? mpiGrid[neighbor]->parameters[CellParams::AMR_ALPHA1] < P::alpha1CoarsenThreshold : true) && (P::useAlpha2 ? mpiGrid[neighbor]->parameters[CellParams::AMR_ALPHA2] < P::alpha2CoarsenThreshold : true))}; - if(!shouldUnrefineNeighbor && - // If the neighbor is planned to remain at the current refinement level, but is outside the allowed refinement region, - // consider it as unrefining instead for purposes of evaluating the neighbors of this cell. - // Induced refinement still possible just beyond that limit. - ((neighborXyz[0] < P::refinementMinX) || (neighborXyz[0] > P::refinementMaxX) - || (neighborXyz[1] < P::refinementMinY) || (neighborXyz[1] > P::refinementMaxY) - || (neighborXyz[2] < P::refinementMinZ) || (neighborXyz[2] > P::refinementMaxZ))) { - shouldUnrefineNeighbor = true; - } + bool shouldRefineNeighbor = shouldRefineCell(mpiGrid, neighbor, r_max2); + bool shouldUnrefineNeighbor = shouldUnrefineCell(mpiGrid, neighbor, r_max2); + int neighborRef {mpiGrid.get_refinement_level(neighbor)}; + if (neighborRef > refLevel && !shouldUnrefineNeighbor) { ++refined_neighbors; } else if (neighborRef < refLevel && !shouldRefineNeighbor) { diff --git a/projects/project.h b/projects/project.h index 2980b8c3e..2d766915b 100644 --- a/projects/project.h +++ b/projects/project.h @@ -85,6 +85,10 @@ namespace projects { virtual bool canRefine(spatial_cell::SpatialCell* cell) const; + virtual bool shouldRefineCell(dccrg::Dccrg& mpiGrid, CellID id, Real r_max2) const; + + virtual bool shouldUnrefineCell(dccrg::Dccrg& mpiGrid, CellID id, Real r_max2) const; + virtual bool refineSpatialCells( dccrg::Dccrg& mpiGrid ) const; /*!\brief Adapts refinement by one level according to the project. Returns true if any cells were refined, false if not. diff --git a/spatial_cell_cpu.cpp b/spatial_cell_cpu.cpp index b7f0f2c0d..3fcb17d5d 100644 --- a/spatial_cell_cpu.cpp +++ b/spatial_cell_cpu.cpp @@ -743,7 +743,7 @@ namespace spatial_cell { // Refinement parameters if ((SpatialCell::mpi_transfer_type & Transfer::REFINEMENT_PARAMETERS)){ displacements.push_back(reinterpret_cast(this->parameters.data() + CellParams::AMR_ALPHA1) - reinterpret_cast(this)); - block_lengths.push_back(sizeof(Real) * (CellParams::AMR_ALPHA2 - CellParams::AMR_ALPHA1 + 1)); // This is just 2, but let's be explicit + block_lengths.push_back(sizeof(Real) * (CellParams::AMR_VORTICITY - CellParams::AMR_ALPHA1 + 1)); // This is just 4, but let's be explicit. } // Copy random number generator state variables diff --git a/spatial_cell_cpu.hpp b/spatial_cell_cpu.hpp index 55281d2fc..d9ec4f5ee 100644 --- a/spatial_cell_cpu.hpp +++ b/spatial_cell_cpu.hpp @@ -316,6 +316,7 @@ namespace spatial_cell { //random_data* get_rng_data_buffer(); // Member variables // + std::array derivativesV; // Derivatives of V for AMR std::array derivativesBVOL; /**< Derivatives of BVOL needed by the acceleration. * Separate array because it does not need to be communicated.*/ //Real parameters[CellParams::N_SPATIAL_CELL_PARAMS]; /**< Bulk variables in this spatial cell.*/ diff --git a/submodules/fsgrid b/submodules/fsgrid index 793a60e68..bf8fe7539 160000 --- a/submodules/fsgrid +++ b/submodules/fsgrid @@ -1 +1 @@ -Subproject commit 793a60e68d648c315e724a0d24cc885262be106c +Subproject commit bf8fe7539e8e589b131fe8f6e35dcfb8802b4dd7 diff --git a/sysboundary/ionosphere.cpp b/sysboundary/ionosphere.cpp index 9a161c29b..9525002fe 100644 --- a/sysboundary/ionosphere.cpp +++ b/sysboundary/ionosphere.cpp @@ -2212,7 +2212,7 @@ namespace SBC { Readparameters::add("ionosphere.atmosphericModelFile", "Filename to read the MSIS atmosphere data from (default: NRLMSIS.dat)", std::string("NRLMSIS.dat")); Readparameters::add("ionosphere.recombAlpha", "Ionospheric recombination parameter (m^3/s)", 2.4e-13); // Default value from Schunck & Nagy, Table 8.5 Readparameters::add("ionosphere.ionizationModel", "Ionospheric electron production rate model. Options are: Rees1963, Rees1989, SergienkoIvanov (default).", std::string("SergienkoIvanov")); - Readparameters::add("ionosphere.innerBoundaryVDFmode", "Inner boundary VDF construction method. Options ar: FixedMoments, AverageMoments, AverageAllMoments, CopyAndLosscone, ForceL2EXB.", std::string("FixedMoments")); + Readparameters::add("ionosphere.innerBoundaryVDFmode", "Inner boundary VDF construction method. Options ar: FixedMoments, AverageMoments, AverageAllMoments, CopyAndLosscone.", std::string("FixedMoments")); Readparameters::add("ionosphere.F10_7", "Solar 10.7 cm radio flux (sfu = 10^{-22} W/m^2)", 100); Readparameters::add("ionosphere.backgroundIonisation", "Background ionoisation due to cosmic rays (mho)", 0.5); Readparameters::add("ionosphere.solverMaxIterations", "Maximum number of iterations for the conjugate gradient solver", 2000); @@ -2271,14 +2271,12 @@ namespace SBC { Readparameters::get("ionosphere.innerBoundaryVDFmode", VDFmodeString); if(VDFmodeString == "FixedMoments") { boundaryVDFmode = FixedMoments; - } else if(VDFmodeString == "AveragMoments") { + } else if(VDFmodeString == "AverageMoments") { boundaryVDFmode = AverageMoments; - } else if(VDFmodeString == "AveragAllMoments") { + } else if(VDFmodeString == "AverageAllMoments") { boundaryVDFmode = AverageAllMoments; } else if(VDFmodeString == "CopyAndLosscone") { boundaryVDFmode = CopyAndLosscone; - } else if(VDFmodeString == "ForceL2EXB") { - boundaryVDFmode = ForceL2EXB; } else { cerr << "(IONOSPHERE) Unknown inner boundary VDF mode \"" << VDFmodeString << "\". Aborting." << endl; abort(); @@ -3098,7 +3096,6 @@ namespace SBC { cellParams[CellParams::BULKV_FORCING_X] = (E[1] * B[2] - E[2] * B[1])/Bsqr; cellParams[CellParams::BULKV_FORCING_Y] = (E[2] * B[0] - E[0] * B[2])/Bsqr; cellParams[CellParams::BULKV_FORCING_Z] = (E[0] * B[1] - E[1] * B[0])/Bsqr; - cellParams[CellParams::FORCING_CELL_NUM]=1; } void Ionosphere::vlasovBoundaryCondition( @@ -3127,20 +3124,6 @@ namespace SBC { switch(boundaryVDFmode) { #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wimplicit-fallthrough" - case ForceL2EXB: - { - // EXB forcing is assigned to the L2 Neighbour cells here, so they can update their VDFs in acceleration - const vector& closeCells = getAllCloseNonsysboundaryCells(cellID); - for (CellID celli : closeCells) { - #pragma omp critical(L2) - { - if(mpiGrid[celli]->parameters[CellParams::FORCING_CELL_NUM] == 0) { - mapCellPotentialAndGetEXBDrift(mpiGrid[celli]->parameters); - } - } - } - // Fall through, to handle L1 in the same way as fixed moments - } case FixedMoments: density = speciesParams[popID].rho; temperature = speciesParams[popID].T; @@ -3184,8 +3167,7 @@ namespace SBC { switch(boundaryVDFmode) { case FixedMoments: case AverageAllMoments: - case AverageMoments: - case ForceL2EXB: + case AverageMoments: { // Fill velocity space with new maxwellian data SpatialCell& cell = *mpiGrid[cellID]; diff --git a/sysboundary/ionosphere.h b/sysboundary/ionosphere.h index 867a52c0f..73f8b6b2d 100644 --- a/sysboundary/ionosphere.h +++ b/sysboundary/ionosphere.h @@ -58,11 +58,10 @@ namespace SBC { }; enum IonosphereBoundaryVDFmode { // How are inner boundary VDFs constructed from the ionosphere - FixedMoments, // Predefine temperature, density and V = 0 on the inner boundary. - AverageMoments, // Copy averaged density and temperature from nearest cells, V = 0 - AverageAllMoments, // Same as above, but also copy V - CopyAndLosscone, - ForceL2EXB + FixedMoments, // Predefine temperature, density and V = EXB drift on the inner boundary. + AverageMoments, // Copy averaged density and temperature from nearest cells, V = EXB drift + AverageAllMoments, // Same as above, but also copy V + add EXB drift to it + CopyAndLosscone }; extern IonosphereBoundaryVDFmode boundaryVDFmode; @@ -384,7 +383,7 @@ namespace SBC { cint k, cuint component ); - // Compute and store the EXB drift into the cell's BULKV_FORCING_X/Y/Z fields and set counter to 1 + // Compute and store the EXB drift into the cell's BULKV_FORCING_X/Y/Z fields virtual void mapCellPotentialAndGetEXBDrift( std::array& cellParams ); diff --git a/sysboundary/sysboundary.cpp b/sysboundary/sysboundary.cpp index 008ba8417..6b355f3dd 100644 --- a/sysboundary/sysboundary.cpp +++ b/sysboundary/sysboundary.cpp @@ -354,11 +354,10 @@ void SysBoundary::classifyCells(dccrg::Dccrg& cells = getLocalCells(); auto localSize = technicalGrid.getLocalSize().data(); - /*set all cells to default value, not_sysboundary and no forcing of the bulkv */ + /*set all cells to default value, not_sysboundary */ #pragma omp parallel for for (uint i = 0; i < cells.size(); i++) { mpiGrid[cells[i]]->sysBoundaryFlag = sysboundarytype::NOT_SYSBOUNDARY; - mpiGrid[cells[i]]->parameters[CellParams::FORCING_CELL_NUM] = -1; } #pragma omp parallel for collapse(2) for (FsGridTools::FsIndex_t z = 0; z < localSize[2]; ++z) { @@ -483,11 +482,6 @@ void SysBoundary::classifyCells(dccrg::Dccrgparameters[CellParams::YCRD]+ mpiGrid[cells[i]]->parameters[CellParams::DY]; creal z = mpiGrid[cells[i]]->parameters[CellParams::ZCRD]+ mpiGrid[cells[i]]->parameters[CellParams::DZ]; creal radius2 = x*x + y*y + z*z; - if ((mpiGrid[cells[i]]->sysBoundaryLayer == 2 || mpiGrid[cells[i]]->sysBoundaryLayer == 1) && - radius2 < ionosphereDownmapRadius*ionosphereDownmapRadius - ) { - mpiGrid[cells[i]]->parameters[CellParams::FORCING_CELL_NUM] = 0; - } } } diff --git a/tools/fileHandling/chunkDataToTape.sh b/tools/fileHandling/chunkDataToTape.sh index 87419e19d..5f8f5ec3b 100755 --- a/tools/fileHandling/chunkDataToTape.sh +++ b/tools/fileHandling/chunkDataToTape.sh @@ -19,7 +19,7 @@ path=$2 start=$3 end=$4 -chunkSize=$((50 * 1024*1024*1024)) +chunkSize=$((5 * 1024*1024*1024)) fileSize=$( /bin/ls -la ${file} | gawk '{ print $5 }' ) chunkNumber=$(( fileSize / chunkSize )) diff --git a/vlasiator.cpp b/vlasiator.cpp index 97dcae141..1e6edc9a6 100644 --- a/vlasiator.cpp +++ b/vlasiator.cpp @@ -289,60 +289,14 @@ void recalculateLocalCellsCache() { Parameters::localCells = mpiGrid.get_cells(); } -int main(int argn,char* args[]) { +int simulate(int argn,char* args[]) { int myRank, doBailout=0; const creal DT_EPSILON=1e-12; typedef Parameters P; Real newDt; bool dtIsChanged {false}; - // Before MPI_Init we hardwire some settings, if we are in OpenMPI - int required=MPI_THREAD_FUNNELED; - int provided, resultlen; - char mpiversion[MPI_MAX_LIBRARY_VERSION_STRING]; - bool overrideMCAompio = false; - - MPI_Get_library_version(mpiversion, &resultlen); - string versionstr = string(mpiversion); - stringstream mpiioMessage; - - if(versionstr.find("Open MPI") != std::string::npos) { - #ifdef VLASIATOR_ALLOW_MCA_OMPIO - mpiioMessage << "We detected OpenMPI but the compilation flag VLASIATOR_ALLOW_MCA_OMPIO was set so we do not override the default MCA io flag." << endl; - #else - overrideMCAompio = true; - int index, count; - char io_value[64]; - MPI_T_cvar_handle io_handle; - - MPI_T_init_thread(required, &provided); - MPI_T_cvar_get_index("io", &index); - MPI_T_cvar_handle_alloc(index, NULL, &io_handle, &count); - MPI_T_cvar_write(io_handle, "^ompio"); - MPI_T_cvar_read(io_handle, io_value); - MPI_T_cvar_handle_free(&io_handle); - mpiioMessage << "We detected OpenMPI so we set the cvars value to disable ompio, MCA io: " << io_value << endl; - #endif - } - - // After the MPI_T settings we can init MPI all right. - MPI_Init_thread(&argn,&args,required,&provided); MPI_Comm_rank(MPI_COMM_WORLD,&myRank); - if (required > provided){ - if(myRank==MASTER_RANK) { - cerr << "(MAIN): MPI_Init_thread failed! Got " << provided << ", need "<, FS_STENCIL_WIDTH> EDt2Grid(fsGridDimensions, MPI_COMM_WORLD, periodicity,P::manualFsGridDecomposition); FsGrid< std::array, FS_STENCIL_WIDTH> EHallGrid(fsGridDimensions, MPI_COMM_WORLD, periodicity,P::manualFsGridDecomposition); FsGrid< std::array, FS_STENCIL_WIDTH> EGradPeGrid(fsGridDimensions, MPI_COMM_WORLD, periodicity,P::manualFsGridDecomposition); + FsGrid< std::array, FS_STENCIL_WIDTH> EGradPeDt2Grid(fsGridDimensions, MPI_COMM_WORLD, periodicity,P::manualFsGridDecomposition); FsGrid< std::array, FS_STENCIL_WIDTH> momentsGrid(fsGridDimensions, MPI_COMM_WORLD, periodicity,P::manualFsGridDecomposition); FsGrid< std::array, FS_STENCIL_WIDTH> momentsDt2Grid(fsGridDimensions, MPI_COMM_WORLD, periodicity,P::manualFsGridDecomposition); FsGrid< std::array, FS_STENCIL_WIDTH> dPerBGrid(fsGridDimensions, MPI_COMM_WORLD, periodicity,P::manualFsGridDecomposition); FsGrid< std::array, FS_STENCIL_WIDTH> dMomentsGrid(fsGridDimensions, MPI_COMM_WORLD, periodicity,P::manualFsGridDecomposition); + FsGrid< std::array, FS_STENCIL_WIDTH> dMomentsDt2Grid(fsGridDimensions, MPI_COMM_WORLD, periodicity,P::manualFsGridDecomposition); FsGrid< std::array, FS_STENCIL_WIDTH> BgBGrid(fsGridDimensions, MPI_COMM_WORLD, periodicity,P::manualFsGridDecomposition); FsGrid< std::array, FS_STENCIL_WIDTH> volGrid(fsGridDimensions, MPI_COMM_WORLD, periodicity,P::manualFsGridDecomposition); FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> technicalGrid(fsGridDimensions, MPI_COMM_WORLD, periodicity,P::manualFsGridDecomposition); @@ -481,19 +437,19 @@ int main(int argn,char* args[]) { // Set DX, DY and DZ // TODO: This is currently just taking the values from cell 1, and assuming them to be // constant throughout the simulation. - perBGrid.DX = perBDt2Grid.DX = EGrid.DX = EDt2Grid.DX = EHallGrid.DX = EGradPeGrid.DX = momentsGrid.DX - = momentsDt2Grid.DX = dPerBGrid.DX = dMomentsGrid.DX = BgBGrid.DX = volGrid.DX = technicalGrid.DX + perBGrid.DX = perBDt2Grid.DX = EGrid.DX = EDt2Grid.DX = EHallGrid.DX = EGradPeGrid.DX = EGradPeDt2Grid.DX = momentsGrid.DX + = momentsDt2Grid.DX = dPerBGrid.DX = dMomentsGrid.DX = dMomentsDt2Grid.DX = BgBGrid.DX = volGrid.DX = technicalGrid.DX = P::dx_ini / pow(2, P::amrMaxSpatialRefLevel); - perBGrid.DY = perBDt2Grid.DY = EGrid.DY = EDt2Grid.DY = EHallGrid.DY = EGradPeGrid.DY = momentsGrid.DY - = momentsDt2Grid.DY = dPerBGrid.DY = dMomentsGrid.DY = BgBGrid.DY = volGrid.DY = technicalGrid.DY + perBGrid.DY = perBDt2Grid.DY = EGrid.DY = EDt2Grid.DY = EHallGrid.DY = EGradPeGrid.DY = EGradPeDt2Grid.DY = momentsGrid.DY + = momentsDt2Grid.DY = dPerBGrid.DY = dMomentsGrid.DY = dMomentsDt2Grid.DY = BgBGrid.DY = volGrid.DY = technicalGrid.DY = P::dy_ini / pow(2, P::amrMaxSpatialRefLevel); - perBGrid.DZ = perBDt2Grid.DZ = EGrid.DZ = EDt2Grid.DZ = EHallGrid.DZ = EGradPeGrid.DZ = momentsGrid.DZ - = momentsDt2Grid.DZ = dPerBGrid.DZ = dMomentsGrid.DZ = BgBGrid.DZ = volGrid.DZ = technicalGrid.DZ + perBGrid.DZ = perBDt2Grid.DZ = EGrid.DZ = EDt2Grid.DZ = EHallGrid.DZ = EGradPeGrid.DZ = EGradPeDt2Grid.DZ = momentsGrid.DZ + = momentsDt2Grid.DZ = dPerBGrid.DZ = dMomentsGrid.DZ = dMomentsDt2Grid.DZ = BgBGrid.DZ = volGrid.DZ = technicalGrid.DZ = P::dz_ini / pow(2, P::amrMaxSpatialRefLevel); // Set the physical start (lower left corner) X, Y, Z perBGrid.physicalGlobalStart = perBDt2Grid.physicalGlobalStart = EGrid.physicalGlobalStart = EDt2Grid.physicalGlobalStart - = EHallGrid.physicalGlobalStart = EGradPeGrid.physicalGlobalStart = momentsGrid.physicalGlobalStart - = momentsDt2Grid.physicalGlobalStart = dPerBGrid.physicalGlobalStart = dMomentsGrid.physicalGlobalStart + = EHallGrid.physicalGlobalStart = EGradPeGrid.physicalGlobalStart = EGradPeDt2Grid.physicalGlobalStart = momentsGrid.physicalGlobalStart + = momentsDt2Grid.physicalGlobalStart = dPerBGrid.physicalGlobalStart = dMomentsGrid.physicalGlobalStart = dMomentsDt2Grid.physicalGlobalStart = BgBGrid.physicalGlobalStart = volGrid.physicalGlobalStart = technicalGrid.physicalGlobalStart = {P::xmin, P::ymin, P::zmin}; @@ -525,6 +481,7 @@ int main(int argn,char* args[]) { BgBGrid, momentsGrid, momentsDt2Grid, + dMomentsGrid, EGrid, EGradPeGrid, volGrid, @@ -586,7 +543,7 @@ int main(int argn,char* args[]) { EHallGrid, EGradPeGrid, momentsGrid, - dPerBGrid, + dPerBGrid, dMomentsGrid, BgBGrid, volGrid, @@ -616,10 +573,12 @@ int main(int argn,char* args[]) { EDt2Grid.finalize(); EHallGrid.finalize(); EGradPeGrid.finalize(); + EGradPeDt2Grid.finalize(); momentsGrid.finalize(); momentsDt2Grid.finalize(); dPerBGrid.finalize(); dMomentsGrid.finalize(); + dMomentsDt2Grid.finalize(); BgBGrid.finalize(); volGrid.finalize(); technicalGrid.finalize(); @@ -639,10 +598,12 @@ int main(int argn,char* args[]) { EDt2Grid, EHallGrid, EGradPeGrid, + EGradPeDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, BgBGrid, volGrid, technicalGrid, @@ -652,7 +613,7 @@ int main(int argn,char* args[]) { phiprof::Timer getFieldsTimer {"getFieldsFromFsGrid"}; volGrid.updateGhostCells(); - getFieldsFromFsGrid(volGrid, BgBGrid, EGradPeGrid, technicalGrid, mpiGrid, cells); + getFieldsFromFsGrid(volGrid, BgBGrid, EGradPeGrid, dMomentsGrid, technicalGrid, mpiGrid, cells); getFieldsTimer.stop(); // Build communicator for ionosphere solving @@ -666,6 +627,7 @@ int main(int argn,char* args[]) { momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, technicalGrid, sysBoundaryContainer, RK_ORDER1, // Update and compute on non-dt2 grids. @@ -790,7 +752,10 @@ int main(int argn,char* args[]) { CellParams::RHOQ, CellParams::P_11, CellParams::P_22, - CellParams::P_33 + CellParams::P_33, + CellParams::P_23, + CellParams::P_13, + CellParams::P_12 ); computeMomentsTimer.stop(); } @@ -1205,7 +1170,10 @@ int main(int argn,char* args[]) { CellParams::RHOQ_DT2, CellParams::P_11_DT2, CellParams::P_22_DT2, - CellParams::P_33_DT2 + CellParams::P_33_DT2, + CellParams::P_23_DT2, + CellParams::P_13_DT2, + CellParams::P_12_DT2 ); momentsTimer.stop(); @@ -1228,10 +1196,12 @@ int main(int argn,char* args[]) { EDt2Grid, EHallGrid, EGradPeGrid, + EGradPeDt2Grid, momentsGrid, momentsDt2Grid, dPerBGrid, dMomentsGrid, + dMomentsDt2Grid, BgBGrid, volGrid, technicalGrid, @@ -1244,7 +1214,7 @@ int main(int argn,char* args[]) { // Copy results back from fsgrid. volGrid.updateGhostCells(); technicalGrid.updateGhostCells(); - getFieldsFromFsGrid(volGrid, BgBGrid, EGradPeGrid, technicalGrid, mpiGrid, cells); + getFieldsFromFsGrid(volGrid, BgBGrid, EGradPeGrid, dMomentsGrid, technicalGrid, mpiGrid, cells); getFieldsTimer.stop(); propagateTimer.stop(cells.size(),"SpatialCells"); addTimedBarrier("barrier-after-field-solver"); @@ -1281,13 +1251,6 @@ int main(int argn,char* args[]) { << endl; SBC::Ionosphere::solveCount++; globalflags::ionosphereJustSolved = true; - // Reset flag in all cells - #pragma omp parallel for - for(size_t i=0; iparameters[CellParams::FORCING_CELL_NUM] == 1) { - mpiGrid[cells[i]]->parameters[CellParams::FORCING_CELL_NUM] = 0; - } - } } phiprof::Timer vspaceTimer {"Velocity-space"}; @@ -1320,7 +1283,10 @@ int main(int argn,char* args[]) { CellParams::RHOQ, CellParams::P_11, CellParams::P_22, - CellParams::P_33 + CellParams::P_33, + CellParams::P_23, + CellParams::P_13, + CellParams::P_12 ); momentsTimer.stop(); @@ -1348,9 +1314,6 @@ int main(int argn,char* args[]) { simulationTimer.stop(); phiprof::Timer finalizationTimer {"Finalization"}; - if (P::propagateField ) { - finalizeFieldPropagator(); - } if (myRank == MASTER_RANK) { if (doBailout > 0) { logFile << "(BAILOUT): Bailing out, see error log for details." << endl; @@ -1377,27 +1340,73 @@ int main(int argn,char* args[]) { phiprof::print(MPI_COMM_WORLD,"phiprof"); - if (myRank == MASTER_RANK) logFile << "(MAIN): Exiting." << endl << writeVerbose; + if (myRank == MASTER_RANK) { + logFile << "(MAIN): Exiting." << endl << writeVerbose; + } logFile.close(); - if (P::diagnosticInterval != 0) diagnostic.close(); + if (P::diagnosticInterval != 0) { + diagnostic.close(); + } + + return 0; +} + +int main(int argn, char* args[]) { + // Before MPI_Init we hardwire some settings, if we are in OpenMPI + int myRank; + int required=MPI_THREAD_FUNNELED; + int provided, resultlen; + char mpiversion[MPI_MAX_LIBRARY_VERSION_STRING]; + bool overrideMCAompio = false; + + MPI_Get_library_version(mpiversion, &resultlen); + string versionstr = string(mpiversion); + stringstream mpiioMessage; + + if(versionstr.find("Open MPI") != std::string::npos) { + #ifdef VLASIATOR_ALLOW_MCA_OMPIO + mpiioMessage << "We detected OpenMPI but the compilation flag VLASIATOR_ALLOW_MCA_OMPIO was set so we do not override the default MCA io flag." << endl; + #else + overrideMCAompio = true; + int index, count; + char io_value[64]; + MPI_T_cvar_handle io_handle; + + MPI_T_init_thread(required, &provided); + MPI_T_cvar_get_index("io", &index); + MPI_T_cvar_handle_alloc(index, NULL, &io_handle, &count); + MPI_T_cvar_write(io_handle, "^ompio"); + MPI_T_cvar_read(io_handle, io_value); + MPI_T_cvar_handle_free(&io_handle); + mpiioMessage << "We detected OpenMPI so we set the cvars value to disable ompio, MCA io: " << io_value << endl; + #endif + } - perBGrid.finalize(); - perBDt2Grid.finalize(); - EGrid.finalize(); - EDt2Grid.finalize(); - EHallGrid.finalize(); - EGradPeGrid.finalize(); - momentsGrid.finalize(); - momentsDt2Grid.finalize(); - dPerBGrid.finalize(); - dMomentsGrid.finalize(); - BgBGrid.finalize(); - volGrid.finalize(); - technicalGrid.finalize(); + // After the MPI_T settings we can init MPI all right. + MPI_Init_thread(&argn,&args,required,&provided); + MPI_Comm_rank(MPI_COMM_WORLD,&myRank); + if (required > provided){ + if(myRank==MASTER_RANK) { + cerr << "(MAIN): MPI_Init_thread failed! Got " << provided << ", need "< compute_acceleration_transformation( } } - // If a bulk velocity is being forced here, perform that last, after things were gyrated in the Hall frame - // If a cell is a remote L2 and was not caught in the loop over neighbours of L1 cells, compute its forcing here - if(globalflags::ionosphereJustSolved - && spatial_cell->parameters[CellParams::FORCING_CELL_NUM] == 0 - && SBC::boundaryVDFmode == SBC::ForceL2EXB - ) { - getObjectWrapper().sysBoundaryContainer.getSysBoundary(sysboundarytype::IONOSPHERE)->mapCellPotentialAndGetEXBDrift(spatial_cell->parameters); // This sets the FORCING_CELL_NUM to 1 - } - if(spatial_cell->parameters[CellParams::FORCING_CELL_NUM] > 0) { - Eigen::Matrix forced_bulkv(spatial_cell->parameters[CellParams::BULKV_FORCING_X], - spatial_cell->parameters[CellParams::BULKV_FORCING_Y], - spatial_cell->parameters[CellParams::BULKV_FORCING_Z]); - - Eigen::Matrix bulkDeltaV = forced_bulkv - bulk_velocity; - total_transform=Translation(bulkDeltaV) * total_transform; - } - return total_transform; } diff --git a/vlasovsolver/cpu_moments.cpp b/vlasovsolver/cpu_moments.cpp index 995dc12f5..7e032674a 100644 --- a/vlasovsolver/cpu_moments.cpp +++ b/vlasovsolver/cpu_moments.cpp @@ -56,6 +56,9 @@ void calculateCellMoments(spatial_cell::SpatialCell* cell, cell->parameters[CellParams::P_11] = 0.0; cell->parameters[CellParams::P_22] = 0.0; cell->parameters[CellParams::P_33] = 0.0; + cell->parameters[CellParams::P_23] = 0.0; + cell->parameters[CellParams::P_13] = 0.0; + cell->parameters[CellParams::P_12] = 0.0; } // Loop over all particle species @@ -125,10 +128,7 @@ void calculateCellMoments(spatial_cell::SpatialCell* cell, const Real mass = getObjectWrapper().particleSpecies[popID].mass; // Temporary array for storing moments - Real array[3]; - for (int i=0; i<3; ++i) { - array[i] = 0.0; - } + std::vector array(6, 0.0); // Calculate species' contribution to second velocity moments Population & pop = cell->get_population(popID); @@ -147,9 +147,12 @@ void calculateCellMoments(spatial_cell::SpatialCell* cell, pop.P[2] = mass*array[2]; if (!computePopulationMomentsOnly) { - cell->parameters[CellParams::P_11] += pop.P[0]; - cell->parameters[CellParams::P_22] += pop.P[1]; - cell->parameters[CellParams::P_33] += pop.P[2]; + cell->parameters[CellParams::P_11] += mass * array[0]; + cell->parameters[CellParams::P_22] += mass * array[1]; + cell->parameters[CellParams::P_33] += mass * array[2]; + cell->parameters[CellParams::P_23] += mass * array[3]; + cell->parameters[CellParams::P_13] += mass * array[4]; + cell->parameters[CellParams::P_12] += mass * array[5]; } } // for-loop over particle species } @@ -187,6 +190,9 @@ void calculateMoments_R( cell->parameters[CellParams::P_11_R] = 0.0; cell->parameters[CellParams::P_22_R] = 0.0; cell->parameters[CellParams::P_33_R] = 0.0; + cell->parameters[CellParams::P_23_R] = 0.0; + cell->parameters[CellParams::P_13_R] = 0.0; + cell->parameters[CellParams::P_12_R] = 0.0; } vmesh::VelocityBlockContainer& blockContainer = cell->get_velocity_blocks(popID); @@ -277,10 +283,7 @@ void calculateMoments_R( const Real mass = getObjectWrapper().particleSpecies[popID].mass; // Temporary array where species' contribution to 2nd moments is accumulated - Real array[3]; - for (int i=0; i<3; ++i) { - array[i] = 0.0; - } + std::vector array(6, 0.0); // Calculate species' contribution to second velocity moments Population & pop = cell->get_population(popID); @@ -298,9 +301,12 @@ void calculateMoments_R( pop.P_R[1] = mass*array[1]; pop.P_R[2] = mass*array[2]; - cell->parameters[CellParams::P_11_R] += pop.P_R[0]; - cell->parameters[CellParams::P_22_R] += pop.P_R[1]; - cell->parameters[CellParams::P_33_R] += pop.P_R[2]; + cell->parameters[CellParams::P_11_R] += mass * array[0]; + cell->parameters[CellParams::P_22_R] += mass * array[1]; + cell->parameters[CellParams::P_33_R] += mass * array[2]; + cell->parameters[CellParams::P_23_R] += mass * array[3]; + cell->parameters[CellParams::P_13_R] += mass * array[4]; + cell->parameters[CellParams::P_12_R] += mass * array[5]; } // for-loop over spatial cells } // for-loop over particle species @@ -341,6 +347,9 @@ void calculateMoments_V( cell->parameters[CellParams::P_11_V] = 0.0; cell->parameters[CellParams::P_22_V] = 0.0; cell->parameters[CellParams::P_33_V] = 0.0; + cell->parameters[CellParams::P_23_V] = 0.0; + cell->parameters[CellParams::P_13_V] = 0.0; + cell->parameters[CellParams::P_12_V] = 0.0; } vmesh::VelocityBlockContainer& blockContainer = cell->get_velocity_blocks(popID); @@ -417,10 +426,7 @@ void calculateMoments_V( const Real mass = getObjectWrapper().particleSpecies[popID].mass; // Temporary array where moments are stored - Real array[3]; - for (int i=0; i<3; ++i) { - array[i] = 0.0; - } + std::vector array(6, 0.0); // Calculate species' contribution to second velocity moments Population & pop = cell->get_population(popID); @@ -439,9 +445,12 @@ void calculateMoments_V( pop.P_V[1] = mass*array[1]; pop.P_V[2] = mass*array[2]; - cell->parameters[CellParams::P_11_V] += pop.P_V[0]; - cell->parameters[CellParams::P_22_V] += pop.P_V[1]; - cell->parameters[CellParams::P_33_V] += pop.P_V[2]; + cell->parameters[CellParams::P_11_V] += mass * array[0]; + cell->parameters[CellParams::P_22_V] += mass * array[1]; + cell->parameters[CellParams::P_33_V] += mass * array[2]; + cell->parameters[CellParams::P_23_V] += mass * array[3]; + cell->parameters[CellParams::P_13_V] += mass * array[4]; + cell->parameters[CellParams::P_12_V] += mass * array[5]; } // for-loop over spatial cells } // for-loop over particle species diff --git a/vlasovsolver/cpu_moments.h b/vlasovsolver/cpu_moments.h index ce98931cb..962920c41 100644 --- a/vlasovsolver/cpu_moments.h +++ b/vlasovsolver/cpu_moments.h @@ -43,7 +43,7 @@ void blockVelocityFirstMoments(const Realf* avgs,const Real* blockParams, template void blockVelocitySecondMoments(const Realf* avgs,const Real* blockParams, const REAL v[3], - REAL* array); + std::vector& array); void calculateMoments_R(dccrg::Dccrg& mpiGrid, const std::vector& cells, @@ -114,13 +114,16 @@ void blockVelocitySecondMoments( const REAL averageVX, const REAL averageVY, const REAL averageVZ, - REAL* array) { + std::vector& array) { const Real HALF = 0.5; Real nvx2_sum = 0.0; Real nvy2_sum = 0.0; Real nvz2_sum = 0.0; + Real nvyvz_sum = 0.0; + Real nvxvz_sum = 0.0; + Real nvxvy_sum = 0.0; for (uint k=0; k& cells = getLocalCells(); @@ -537,6 +540,9 @@ void calculateInterpolatedVelocityMoments( SC->parameters[cp_p11] = 0.5* ( SC->parameters[CellParams::P_11_R] + SC->parameters[CellParams::P_11_V] ); SC->parameters[cp_p22] = 0.5* ( SC->parameters[CellParams::P_22_R] + SC->parameters[CellParams::P_22_V] ); SC->parameters[cp_p33] = 0.5* ( SC->parameters[CellParams::P_33_R] + SC->parameters[CellParams::P_33_V] ); + SC->parameters[cp_p23] = 0.5* ( SC->parameters[CellParams::P_23_R] + SC->parameters[CellParams::P_23_V] ); + SC->parameters[cp_p13] = 0.5* ( SC->parameters[CellParams::P_13_R] + SC->parameters[CellParams::P_13_V] ); + SC->parameters[cp_p12] = 0.5* ( SC->parameters[CellParams::P_12_R] + SC->parameters[CellParams::P_12_V] ); for (uint popID=0; popIDget_population(popID);