From fb20eec9c85e028f964ec913892f016db584a043 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Wed, 28 Aug 2024 15:27:12 +0300 Subject: [PATCH 01/30] Off-diagonal pressure calculation --- common.h | 9 ++++++++ vlasovsolver/cpu_moments.cpp | 42 ++++++++++++++++++------------------ vlasovsolver/cpu_moments.h | 14 ++++++++++-- 3 files changed, 42 insertions(+), 23 deletions(-) diff --git a/common.h b/common.h index 7f35fa601..86472fac9 100644 --- a/common.h +++ b/common.h @@ -165,15 +165,24 @@ 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_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.*/ diff --git a/vlasovsolver/cpu_moments.cpp b/vlasovsolver/cpu_moments.cpp index 995dc12f5..8d00f87db 100644 --- a/vlasovsolver/cpu_moments.cpp +++ b/vlasovsolver/cpu_moments.cpp @@ -125,10 +125,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 +144,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 } @@ -277,10 +277,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 +295,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 @@ -417,10 +417,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 +436,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[4]; + cell->parameters[CellParams::P_13_V] += mass * array[5]; + cell->parameters[CellParams::P_12_V] += mass * array[6]; } // 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 Date: Fri, 30 Aug 2024 07:12:32 +0300 Subject: [PATCH 02/30] Pressure anisotropy calculation --- common.h | 1 + fieldsolver/derivatives.cpp | 17 +++++++++++++++++ 2 files changed, 18 insertions(+) diff --git a/common.h b/common.h index 86472fac9..722d34d08 100644 --- a/common.h +++ b/common.h @@ -219,6 +219,7 @@ namespace CellParams { AMR_ALPHA1, AMR_ALPHA2, RECENTLY_REFINED, + P_ANISOTROPY, 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) */ diff --git a/fieldsolver/derivatives.cpp b/fieldsolver/derivatives.cpp index b1aad577c..0597e369f 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. * @@ -801,6 +802,22 @@ 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{0, 0, 1}, Eigen::Vector3d{myB[0], myB[1], myB[2]}).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 Pprime = rot * P * rot.transpose(); + cell->parameters[CellParams::P_ANISOTROPY] = (Pprime(0, 0) + Pprime(1, 1)) / (2 * Pprime(2, 2)); + + // Vorticity: dmoments + // Pressure: devise a rotation matrix for B, see zulip and good old wikipedia + cell->parameters[CellParams::AMR_DRHO] = dRho; cell->parameters[CellParams::AMR_DU] = dU; cell->parameters[CellParams::AMR_DPSQ] = dPsq; From b5356f5738e224c60a3d1b948df71329871a0dd3 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Fri, 30 Aug 2024 07:12:42 +0300 Subject: [PATCH 03/30] Implement pressure anisotropy in refinement --- projects/project.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/projects/project.cpp b/projects/project.cpp index a4b174fec..99f3d8d20 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -598,6 +598,13 @@ namespace projects { 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))}; + // Pressure anisotropy forces refinement + // TODO make this a config thingie + if (cell->parameters[CellParams::P_ANISOTROPY] < 0.5) { + shouldRefine = true; + shouldUnrefine = false; + } + 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. From d65e2cc9e896d1b73864f5b9af990736dcd71c50 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Fri, 30 Aug 2024 11:56:26 +0300 Subject: [PATCH 04/30] Off by one --- vlasovsolver/cpu_moments.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/vlasovsolver/cpu_moments.cpp b/vlasovsolver/cpu_moments.cpp index 8d00f87db..184f50253 100644 --- a/vlasovsolver/cpu_moments.cpp +++ b/vlasovsolver/cpu_moments.cpp @@ -439,9 +439,9 @@ void calculateMoments_V( 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[4]; - cell->parameters[CellParams::P_13_V] += mass * array[5]; - cell->parameters[CellParams::P_12_V] += mass * array[6]; + 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 From ff3baa159d58f86791947d3d99db7e6a568fe0f1 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Fri, 30 Aug 2024 12:09:15 +0300 Subject: [PATCH 05/30] Pressure anisotropy threshold and DROs --- datareduction/datareducer.cpp | 9 ++++++++- parameters.cpp | 5 ++++- parameters.h | 1 + projects/project.cpp | 2 +- 4 files changed, 14 insertions(+), 3 deletions(-) diff --git a/datareduction/datareducer.cpp b/datareduction/datareducer.cpp index 309915ae9..f41f31d31 100644 --- a/datareduction/datareducer.cpp +++ b/datareduction/datareducer.cpp @@ -2962,7 +2962,14 @@ 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; } diff --git a/parameters.cpp b/parameters.cpp index 6fa269c92..035b0fa18 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -167,6 +167,7 @@ Real P::alpha1CoarsenThreshold = -1.0; bool P::useAlpha2 = true; Real P::alpha2RefineThreshold = 0.5; Real P::alpha2CoarsenThreshold = -1.0; +Real P::anisotropyThreshold = -1; Real P::alphaDRhoWeight = 1.0; Real P::alphaDUWeight = 1.0; Real P::alphaDPSqWeight = 1.0; @@ -409,7 +410,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"); RP::addComposing( "variables_deprecated.output", @@ -484,6 +485,7 @@ 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.anisotropy_threshold","Determines the maximum value of pressure anisotropy to refine cells", -1.0); 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); @@ -800,6 +802,7 @@ void Parameters::getParameters() { } MPI_Abort(MPI_COMM_WORLD, 1); } + RP::get("AMR.anisotropy_threshold", P::anisotropyThreshold); RP::get("AMR.refine_cadence",P::refineCadence); RP::get("AMR.refine_after",P::refineAfter); diff --git a/parameters.h b/parameters.h index 87838ef9f..9616559c4 100644 --- a/parameters.h +++ b/parameters.h @@ -198,6 +198,7 @@ struct Parameters { static bool useAlpha2; static Real alpha2RefineThreshold; static Real alpha2CoarsenThreshold; + static Real anisotropyThreshold; static uint refineCadence; static Real refineAfter; static Real refineRadius; diff --git a/projects/project.cpp b/projects/project.cpp index 99f3d8d20..ef65ef88d 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -600,7 +600,7 @@ namespace projects { // Pressure anisotropy forces refinement // TODO make this a config thingie - if (cell->parameters[CellParams::P_ANISOTROPY] < 0.5) { + if (P::anisotropyThreshold > 0 && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyThreshold) { shouldRefine = true; shouldUnrefine = false; } From dfdcaa12c0f078b0448c36012f44f07ae79a97b6 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Fri, 30 Aug 2024 12:25:03 +0300 Subject: [PATCH 06/30] Interpolate off-diagonal pressure --- common.h | 3 +++ vlasiator.cpp | 15 ++++++++++++--- vlasovmover.h | 5 ++++- vlasovsolver/vlasovmover.cpp | 8 +++++++- 4 files changed, 26 insertions(+), 5 deletions(-) diff --git a/common.h b/common.h index 722d34d08..79e28a729 100644 --- a/common.h +++ b/common.h @@ -171,6 +171,9 @@ namespace CellParams { 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 */ diff --git a/vlasiator.cpp b/vlasiator.cpp index c8f1df0ac..07233ea3b 100644 --- a/vlasiator.cpp +++ b/vlasiator.cpp @@ -787,7 +787,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(); } @@ -1202,7 +1205,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(); @@ -1317,7 +1323,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(); diff --git a/vlasovmover.h b/vlasovmover.h index 9289d8396..85dd72401 100644 --- a/vlasovmover.h +++ b/vlasovmover.h @@ -64,7 +64,10 @@ void calculateInterpolatedVelocityMoments( const int cp_rhoq, const int cp_p11, const int cp_p22, - const int cp_p33 + const int cp_p33, + const int cp_p23, + const int cp_p13, + const int cp_p12 ); /*! diff --git a/vlasovsolver/vlasovmover.cpp b/vlasovsolver/vlasovmover.cpp index 0332df94f..62e290a9c 100644 --- a/vlasovsolver/vlasovmover.cpp +++ b/vlasovsolver/vlasovmover.cpp @@ -520,7 +520,10 @@ void calculateInterpolatedVelocityMoments( const int cp_rhoq, const int cp_p11, const int cp_p22, - const int cp_p33 + const int cp_p33, + const int cp_p23, + const int cp_p13, + const int cp_p12 ) { const vector& 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); From 16e6f23a449598f2acbdea4bb6fbb786ee5b6909 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Mon, 2 Sep 2024 16:32:14 +0300 Subject: [PATCH 07/30] Vorticity calculations --- common.h | 19 ++++++++++++++++++- fieldsolver/derivatives.cpp | 15 ++++++++++++--- fieldsolver/gridGlue.cpp | 29 +++++++++++++++++++++++++++++ fieldsolver/gridGlue.hpp | 25 ++++++++++++++++++------- grid.cpp | 3 ++- grid.h | 1 + spatial_cell_cpu.hpp | 1 + vlasiator.cpp | 5 +++-- warnings.txt | 34 ++++++++++++++++++++++++++++++++++ 9 files changed, 118 insertions(+), 14 deletions(-) create mode 100644 warnings.txt diff --git a/common.h b/common.h index 79e28a729..c734d4513 100644 --- a/common.h +++ b/common.h @@ -221,8 +221,9 @@ namespace CellParams { AMR_DB, AMR_ALPHA1, AMR_ALPHA2, - RECENTLY_REFINED, 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) */ @@ -249,6 +250,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/fieldsolver/derivatives.cpp b/fieldsolver/derivatives.cpp index 0597e369f..ca97fc986 100644 --- a/fieldsolver/derivatives.cpp +++ b/fieldsolver/derivatives.cpp @@ -813,10 +813,16 @@ void calculateScaledDeltas( }; Eigen::Matrix3d Pprime = rot * P * rot.transpose(); - cell->parameters[CellParams::P_ANISOTROPY] = (Pprime(0, 0) + Pprime(1, 1)) / (2 * Pprime(2, 2)); - // Vorticity: dmoments - // Pressure: devise a rotation matrix for B, see zulip and good old wikipedia + // 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 velocity {std::sqrt(std::pow(myP[0], 2) + std::pow(myP[1], 2) + std::pow(myP[2], 2)) / myRho}; cell->parameters[CellParams::AMR_DRHO] = dRho; cell->parameters[CellParams::AMR_DU] = dU; @@ -825,6 +831,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] = (Pprime(0, 0) + Pprime(1, 1)) / (2 * Pprime(2, 2)); + // Experimental, consider other values for scaling (term here is V = myP / myRho) + cell->parameters[CellParams::AMR_VORTICITY] = vorticity * cell->parameters[CellParams::DX] / (velocity + EPS); } /*! \brief High-level scaled gradient calculation wrapper function. diff --git a/fieldsolver/gridGlue.cpp b/fieldsolver/gridGlue.cpp index 9a0f9a7e2..4979530d5 100644 --- a/fieldsolver/gridGlue.cpp +++ b/fieldsolver/gridGlue.cpp @@ -305,6 +305,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 @@ -396,6 +397,7 @@ 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); sendBuffer[ii].sums[FieldsToCommunicate::PERBXVOL] += volcell->at(fsgrids::volfields::PERBXVOL); sendBuffer[ii].sums[FieldsToCommunicate::PERBYVOL] += volcell->at(fsgrids::volfields::PERBYVOL); @@ -409,6 +411,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); @@ -468,6 +479,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; @@ -495,6 +515,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 c9cbe8fce..1d1d07b85 100644 --- a/fieldsolver/gridGlue.hpp +++ b/fieldsolver/gridGlue.hpp @@ -29,6 +29,15 @@ enum FieldsToCommunicate { CURVATUREX, CURVATUREY, CURVATUREZ, + dVxdx, + dVxdy, + dVxdz, + dVydx, + dVydy, + dVydz, + dVzdx, + dVzdy, + dVzdz, N_FIELDSTOCOMMUNICATE }; @@ -56,13 +65,15 @@ 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 * \param mpiGrid The DCCRG grid carrying fields. diff --git a/grid.cpp b/grid.cpp index ec1bf3c20..4e92bdd54 100644 --- a/grid.cpp +++ b/grid.cpp @@ -94,6 +94,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, @@ -353,7 +354,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(); diff --git a/grid.h b/grid.h index f124cdf99..1cd1c759b 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/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/vlasiator.cpp b/vlasiator.cpp index 07233ea3b..7331d9c22 100644 --- a/vlasiator.cpp +++ b/vlasiator.cpp @@ -522,6 +522,7 @@ int main(int argn,char* args[]) { BgBGrid, momentsGrid, momentsDt2Grid, + dMomentsGrid, EGrid, EGradPeGrid, volGrid, @@ -649,7 +650,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 @@ -1247,7 +1248,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"); diff --git a/warnings.txt b/warnings.txt new file mode 100644 index 000000000..6642ed85f --- /dev/null +++ b/warnings.txt @@ -0,0 +1,34 @@ +./generate_version.sh: line 78: xxd: command not found +./generate_version.sh: line 141: xxd: command not found +make: *** [Makefile:235: datareducer.o] Interrupt +make: *** [Makefile:235: datareductionoperator.o] Interrupt +make: *** [Makefile:225: vamr_refinement_criteria.o] Interrupt +make: *** [Makefile:240: donotcompute.o] Interrupt +make: *** [Makefile:240: ionosphere.o] Interrupt +make: *** [Makefile:240: copysphere.o] Interrupt +make: *** [Makefile:240: outflow.o] Interrupt +make: *** [Makefile:240: inflow.o] Interrupt +make: *** [Makefile:240: setmaxwellian.o] Interrupt +make: *** [Makefile:245: fieldtracing.o] Interrupt +make: *** [Makefile:240: sysboundary.o] Interrupt +make: *** [Makefile:240: sysboundarycondition.o] Interrupt +make: *** [Makefile:250: project.o] Interrupt +make: *** [Makefile:250: projectTriAxisSearch.o] Interrupt +make: *** [Makefile:256: Alfven.o] Interrupt +make: *** [Makefile:256: Diffusion.o] Interrupt +make: *** [Makefile:256: Dispersion.o] Interrupt +make: *** [Makefile:256: Distributions.o] Interrupt +make: *** [Makefile:256: Firehose.o] Interrupt +make: *** [Makefile:256: Flowthrough.o] Interrupt +make: *** [Makefile:256: Fluctuations.o] Interrupt +make: *** [Makefile:256: Harris.o] Interrupt +make: *** [Makefile:256: KHB.o] Interrupt +make: *** [Makefile:256: Larmor.o] Interrupt +make: *** [Makefile:256: Magnetosphere.o] Interrupt +make: *** [Makefile:256: MultiPeak.o] Interrupt +make: *** [Makefile:256: Riemann1.o] Interrupt +make: *** [Makefile:256: Shock.o] Interrupt +make: *** [Makefile:256: Template.o] Interrupt +make: *** [Makefile:256: test_fp.o] Interrupt +make: *** [Makefile:256: testHall.o] Interrupt +make: *** [Makefile:256: test_trans.o] Interrupt From 6474f3ef0018cc68c1b26b718c8c384a125924e6 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Mon, 2 Sep 2024 16:43:57 +0300 Subject: [PATCH 08/30] Integrate vorticity with AMR --- datareduction/datareducer.cpp | 7 +++++++ parameters.cpp | 23 +++++++++++++++++++++- parameters.h | 4 ++++ projects/project.cpp | 36 +++++++++++++++++++++++++++++------ 4 files changed, 63 insertions(+), 7 deletions(-) diff --git a/datareduction/datareducer.cpp b/datareduction/datareducer.cpp index f41f31d31..dda5607d0 100644 --- a/datareduction/datareducer.cpp +++ b/datareduction/datareducer.cpp @@ -2974,6 +2974,13 @@ void initializeDataReducers(DataReducer * outputReducer, DataReducer * diagnosti 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; + } + } if(P::systemWriteAllDROs || lowercase == "ig_latitude") { outputReducer->addOperator(new DRO::DataReductionOperatorIonosphereNode("ig_latitude", []( SBC::SphericalTriGrid& grid)->std::vector { diff --git a/parameters.cpp b/parameters.cpp index 035b0fa18..7c2bddf07 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -167,6 +167,9 @@ Real P::alpha1CoarsenThreshold = -1.0; bool P::useAlpha2 = true; Real P::alpha2RefineThreshold = 0.5; Real P::alpha2CoarsenThreshold = -1.0; +bool P::useVorticity = true; +Real P::vorticityRefineThreshold = 0.5; +Real P::vorticityCoarsenThreshold = -1.0; Real P::anisotropyThreshold = -1; Real P::alphaDRhoWeight = 1.0; Real P::alphaDUWeight = 1.0; @@ -410,7 +413,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_pressure_anisotropy"); + "vg_amr_translate_comm vg_gridcoordinates fg_gridcoordinates vg_pressure_anisotropy vg_amr_vorticity"); RP::addComposing( "variables_deprecated.output", @@ -485,6 +488,9 @@ 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_alpha2","Use vorticity as a refinement index", true); + RP::add("AMR.alpha2_refine_threshold","Determines the minimum value of vorticity to refine cells", 0.5); + RP::add("AMR.alpha2_coarsen_threshold","Determines the maximum value of vorticity to unrefine cells, default half of the refine threshold", -1.0); RP::add("AMR.anisotropy_threshold","Determines the maximum value of pressure anisotropy to refine cells", -1.0); 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); @@ -790,6 +796,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); @@ -802,6 +809,20 @@ 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.anisotropy_threshold", P::anisotropyThreshold); RP::get("AMR.refine_cadence",P::refineCadence); diff --git a/parameters.h b/parameters.h index 9616559c4..ee4e6b20d 100644 --- a/parameters.h +++ b/parameters.h @@ -198,6 +198,10 @@ 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 Real anisotropyThreshold; static uint refineCadence; static Real refineAfter; diff --git a/projects/project.cpp b/projects/project.cpp index ef65ef88d..1c8d8ccfb 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -595,8 +595,21 @@ namespace projects { // 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))}; + // TODO starting to look truly cursed. Might need refactoring + 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) || + (P::useVorticity ? cell->parameters[CellParams::AMR_VORTICITY] > P::vorticityRefineThreshold : 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) && + (P::useVorticity ? cell->parameters[CellParams::AMR_VORTICITY] < P::vorticityCoarsenThreshold : true) + ) + }; // Pressure anisotropy forces refinement // TODO make this a config thingie @@ -632,8 +645,13 @@ namespace projects { 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 && + 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) || + (P::useVorticity ? mpiGrid[neighbor]->parameters[CellParams::AMR_VORTICITY] > P::vorticityRefineThreshold : 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) @@ -642,8 +660,14 @@ namespace projects { 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 && + 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) && + (P::useVorticity ? mpiGrid[neighbor]->parameters[CellParams::AMR_VORTICITY] < P::vorticityCoarsenThreshold : 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. From 4449ff5b8072104fe98e6d2f853272a3c5cc7475 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Mon, 2 Sep 2024 16:46:21 +0300 Subject: [PATCH 09/30] Rename anisotropyThreshold --- parameters.cpp | 6 +++--- parameters.h | 2 +- projects/project.cpp | 5 ++--- 3 files changed, 6 insertions(+), 7 deletions(-) diff --git a/parameters.cpp b/parameters.cpp index 7c2bddf07..ee6667adc 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -170,7 +170,7 @@ Real P::alpha2CoarsenThreshold = -1.0; bool P::useVorticity = true; Real P::vorticityRefineThreshold = 0.5; Real P::vorticityCoarsenThreshold = -1.0; -Real P::anisotropyThreshold = -1; +Real P::anisotropyRefineThreshold = -1; Real P::alphaDRhoWeight = 1.0; Real P::alphaDUWeight = 1.0; Real P::alphaDPSqWeight = 1.0; @@ -491,7 +491,7 @@ bool P::addParameters() { RP::add("AMR.use_alpha2","Use vorticity as a refinement index", true); RP::add("AMR.alpha2_refine_threshold","Determines the minimum value of vorticity to refine cells", 0.5); RP::add("AMR.alpha2_coarsen_threshold","Determines the maximum value of vorticity to unrefine cells, default half of the refine threshold", -1.0); - RP::add("AMR.anisotropy_threshold","Determines the maximum value of pressure anisotropy to refine cells", -1.0); + RP::add("AMR.anisotropyRefineThreshold","Determines the maximum value of pressure anisotropy to refine cells", -1.0); 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); @@ -823,7 +823,7 @@ void Parameters::getParameters() { MPI_Abort(MPI_COMM_WORLD, 1); } - RP::get("AMR.anisotropy_threshold", P::anisotropyThreshold); + RP::get("AMR.anisotropy_refine_threshold", P::anisotropyRefineThreshold); RP::get("AMR.refine_cadence",P::refineCadence); RP::get("AMR.refine_after",P::refineAfter); diff --git a/parameters.h b/parameters.h index ee4e6b20d..3c21f87d1 100644 --- a/parameters.h +++ b/parameters.h @@ -202,7 +202,7 @@ struct Parameters { static bool useVorticity; static Real vorticityRefineThreshold; static Real vorticityCoarsenThreshold; - static Real anisotropyThreshold; + static Real anisotropyRefineThreshold; static uint refineCadence; static Real refineAfter; static Real refineRadius; diff --git a/projects/project.cpp b/projects/project.cpp index 1c8d8ccfb..b67005f07 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -611,9 +611,8 @@ namespace projects { ) }; - // Pressure anisotropy forces refinement - // TODO make this a config thingie - if (P::anisotropyThreshold > 0 && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyThreshold) { + // Pressure anisotropy forces refinement when under threshold + if (P::anisotropyRefineThreshold > 0 && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold) { shouldRefine = true; shouldUnrefine = false; } From 02bcb200d4b2a8c505b76504fa4fe3ec6b4cf23e Mon Sep 17 00:00:00 2001 From: lkotipal Date: Mon, 2 Sep 2024 16:55:30 +0300 Subject: [PATCH 10/30] oops --- parameters.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/parameters.cpp b/parameters.cpp index ee6667adc..a0209f316 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -488,9 +488,9 @@ 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_alpha2","Use vorticity as a refinement index", true); - RP::add("AMR.alpha2_refine_threshold","Determines the minimum value of vorticity to refine cells", 0.5); - RP::add("AMR.alpha2_coarsen_threshold","Determines the maximum value of vorticity to unrefine cells, default half of the refine threshold", -1.0); + RP::add("AMR.use_vorticity","Use vorticity as a refinement index", true); + 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.anisotropyRefineThreshold","Determines the maximum value of pressure anisotropy to refine cells", -1.0); 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); From d2ac51e19c6e82edf2daba1698d33e7931045617 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Tue, 3 Sep 2024 11:01:14 +0300 Subject: [PATCH 11/30] Oops x2 --- parameters.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parameters.cpp b/parameters.cpp index a0209f316..cdba63a41 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -491,7 +491,7 @@ bool P::addParameters() { RP::add("AMR.use_vorticity","Use vorticity as a refinement index", true); 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.anisotropyRefineThreshold","Determines the maximum value of pressure anisotropy to refine cells", -1.0); + RP::add("AMR.anisotropy_refine_threshold","Determines the maximum value of pressure anisotropy to refine cells", -1.0); 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); From decc207f9d1dd496d359deeead6f36d21cbbe827 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Wed, 4 Sep 2024 14:14:50 +0300 Subject: [PATCH 12/30] Scale vorticity with Alfven velocity --- fieldsolver/derivatives.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/fieldsolver/derivatives.cpp b/fieldsolver/derivatives.cpp index ca97fc986..d5f8d36f2 100644 --- a/fieldsolver/derivatives.cpp +++ b/fieldsolver/derivatives.cpp @@ -822,7 +822,7 @@ void calculateScaledDeltas( 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 velocity {std::sqrt(std::pow(myP[0], 2) + std::pow(myP[1], 2) + std::pow(myP[2], 2)) / myRho}; + Real vA {std::sqrt(Bsq / (physicalconstants::MU_0 * myRho))}; cell->parameters[CellParams::AMR_DRHO] = dRho; cell->parameters[CellParams::AMR_DU] = dU; @@ -832,8 +832,8 @@ void calculateScaledDeltas( 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] = (Pprime(0, 0) + Pprime(1, 1)) / (2 * Pprime(2, 2)); - // Experimental, consider other values for scaling (term here is V = myP / myRho) - cell->parameters[CellParams::AMR_VORTICITY] = vorticity * cell->parameters[CellParams::DX] / (velocity + EPS); + // Experimental, current scaling is Alfven velocity + cell->parameters[CellParams::AMR_VORTICITY] = vorticity * cell->parameters[CellParams::DX] / (vA + EPS); } /*! \brief High-level scaled gradient calculation wrapper function. From f1e6125306adeecad7b677c6c536cbde4e9fcee2 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Thu, 12 Sep 2024 13:55:39 +0300 Subject: [PATCH 13/30] Maximum reflevel for anisotropy based refinement --- parameters.cpp | 3 +++ parameters.h | 1 + projects/project.cpp | 2 +- 3 files changed, 5 insertions(+), 1 deletion(-) diff --git a/parameters.cpp b/parameters.cpp index cdba63a41..88ee037d2 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -171,6 +171,7 @@ bool P::useVorticity = true; Real P::vorticityRefineThreshold = 0.5; Real P::vorticityCoarsenThreshold = -1.0; Real P::anisotropyRefineThreshold = -1; +int P::anisotropyMaxReflevel = 2; Real P::alphaDRhoWeight = 1.0; Real P::alphaDUWeight = 1.0; Real P::alphaDPSqWeight = 1.0; @@ -492,6 +493,7 @@ bool P::addParameters() { 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.anisotropy_refine_threshold","Determines the maximum value of pressure anisotropy to refine cells", -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); @@ -824,6 +826,7 @@ void Parameters::getParameters() { } RP::get("AMR.anisotropy_refine_threshold", P::anisotropyRefineThreshold); + RP::get("AMR.anisotropy_max_reflevel", P::anisotropyMaxReflevel); RP::get("AMR.refine_cadence",P::refineCadence); RP::get("AMR.refine_after",P::refineAfter); diff --git a/parameters.h b/parameters.h index 3c21f87d1..eb1534863 100644 --- a/parameters.h +++ b/parameters.h @@ -203,6 +203,7 @@ struct Parameters { static Real vorticityRefineThreshold; static Real vorticityCoarsenThreshold; static Real anisotropyRefineThreshold; + static int anisotropyMaxReflevel; static uint refineCadence; static Real refineAfter; static Real refineRadius; diff --git a/projects/project.cpp b/projects/project.cpp index b67005f07..e93ce8537 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -612,7 +612,7 @@ namespace projects { }; // Pressure anisotropy forces refinement when under threshold - if (P::anisotropyRefineThreshold > 0 && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold) { + if (P::anisotropyRefineThreshold > 0 && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && refLevel < P::anisotropyMaxReflevel) { shouldRefine = true; shouldUnrefine = false; } From 2a602f1d0f3d13f84715f7bc7520d1f7d4e70f6b Mon Sep 17 00:00:00 2001 From: lkotipal Date: Thu, 12 Sep 2024 13:57:24 +0300 Subject: [PATCH 14/30] Bulk velocity scaling for vorticity --- fieldsolver/derivatives.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/fieldsolver/derivatives.cpp b/fieldsolver/derivatives.cpp index d5f8d36f2..affea1140 100644 --- a/fieldsolver/derivatives.cpp +++ b/fieldsolver/derivatives.cpp @@ -822,7 +822,8 @@ void calculateScaledDeltas( 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 vA {std::sqrt(Bsq / (physicalconstants::MU_0 * myRho))}; + Real myV {std::sqrt(std::pow(myP[0], 2) + std::pow(myP[1], 2) + std::pow(myP[2], 2)) / myRho}; cell->parameters[CellParams::AMR_DRHO] = dRho; cell->parameters[CellParams::AMR_DU] = dU; @@ -832,8 +833,8 @@ void calculateScaledDeltas( 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] = (Pprime(0, 0) + Pprime(1, 1)) / (2 * Pprime(2, 2)); - // Experimental, current scaling is Alfven velocity - cell->parameters[CellParams::AMR_VORTICITY] = vorticity * cell->parameters[CellParams::DX] / (vA + EPS); + // Experimental, current scaling is bulk velocity + cell->parameters[CellParams::AMR_VORTICITY] = vorticity * cell->parameters[CellParams::DX] / (myV + EPS); } /*! \brief High-level scaled gradient calculation wrapper function. From 3693c1b72623f4a34c9ee311801cc2f76e1c833f Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Wed, 18 Sep 2024 10:44:48 +0300 Subject: [PATCH 15/30] Make anisotropy AMR obey rmax and communicate vorticity and anisotropy --- projects/project.cpp | 5 +++-- spatial_cell_cpu.cpp | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/projects/project.cpp b/projects/project.cpp index e93ce8537..d06321c75 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -612,7 +612,7 @@ namespace projects { }; // Pressure anisotropy forces refinement when under threshold - if (P::anisotropyRefineThreshold > 0 && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && refLevel < P::anisotropyMaxReflevel) { + if (r2 < r_max2 && P::anisotropyRefineThreshold > 0 && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && refLevel < P::anisotropyMaxReflevel) { shouldRefine = true; shouldUnrefine = false; } @@ -648,7 +648,8 @@ namespace projects { (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) || - (P::useVorticity ? mpiGrid[neighbor]->parameters[CellParams::AMR_VORTICITY] > P::vorticityRefineThreshold : false) + (P::useVorticity ? mpiGrid[neighbor]->parameters[CellParams::AMR_VORTICITY] > P::vorticityRefineThreshold : false) || + (P::anisotropyRefineThreshold > 0 && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && refLevel < P::anisotropyMaxReflevel) )}; if (shouldRefineNeighbor && // If the neighbor is planned to be refined, but is outside the allowed refinement region, cancel that refinement. diff --git a/spatial_cell_cpu.cpp b/spatial_cell_cpu.cpp index b7f0f2c0d..1714d7f1a 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 2, but let's be explicit } // Copy random number generator state variables From ac5ddcceda8be2b733cb9b4489e886d554c7e6ea Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Wed, 18 Sep 2024 10:44:48 +0300 Subject: [PATCH 16/30] Make anisotropy AMR obey rmax and communicate vorticity and anisotropy From 01113dd15aadbea239d630f0d27cb3b3b0b85eb6 Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Wed, 18 Sep 2024 11:07:20 +0300 Subject: [PATCH 17/30] Fixed bugs and an inaccurate comment --- projects/project.cpp | 2 +- spatial_cell_cpu.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/projects/project.cpp b/projects/project.cpp index d06321c75..fc3655826 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -649,7 +649,7 @@ namespace projects { (P::useAlpha1 ? mpiGrid[neighbor]->parameters[CellParams::AMR_ALPHA1] > P::alpha1RefineThreshold : false) || (P::useAlpha2 ? mpiGrid[neighbor]->parameters[CellParams::AMR_ALPHA2] > P::alpha2RefineThreshold : false) || (P::useVorticity ? mpiGrid[neighbor]->parameters[CellParams::AMR_VORTICITY] > P::vorticityRefineThreshold : false) || - (P::anisotropyRefineThreshold > 0 && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && refLevel < P::anisotropyMaxReflevel) + (P::anisotropyRefineThreshold > 0 && mpiGrid[neighbor]->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && neighborRef < P::anisotropyMaxReflevel) )}; if (shouldRefineNeighbor && // If the neighbor is planned to be refined, but is outside the allowed refinement region, cancel that refinement. diff --git a/spatial_cell_cpu.cpp b/spatial_cell_cpu.cpp index 1714d7f1a..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_VORTICITY - 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 From ffde733e3d47534542656d9de5e88959afa77c19 Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Tue, 24 Sep 2024 11:02:21 +0300 Subject: [PATCH 18/30] Fix divide by zero in vorticity and anisotropy DROs --- fieldsolver/derivatives.cpp | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/fieldsolver/derivatives.cpp b/fieldsolver/derivatives.cpp index affea1140..cdf19055e 100644 --- a/fieldsolver/derivatives.cpp +++ b/fieldsolver/derivatives.cpp @@ -740,15 +740,22 @@ void calculateScaledDeltas( Real myRho {cell->parameters[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); @@ -814,6 +821,11 @@ void calculateScaledDeltas( Eigen::Matrix3d Pprime = rot * P * rot.transpose(); + 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]}; @@ -823,7 +835,10 @@ void calculateScaledDeltas( 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 myV {std::sqrt(std::pow(myP[0], 2) + std::pow(myP[1], 2) + std::pow(myP[2], 2)) / 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; @@ -832,9 +847,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] = (Pprime(0, 0) + Pprime(1, 1)) / (2 * Pprime(2, 2)); + cell->parameters[CellParams::P_ANISOTROPY] = Panisotropy; // Experimental, current scaling is bulk velocity - cell->parameters[CellParams::AMR_VORTICITY] = vorticity * cell->parameters[CellParams::DX] / (myV + EPS); + cell->parameters[CellParams::AMR_VORTICITY] = amr_vorticity; } /*! \brief High-level scaled gradient calculation wrapper function. From 73717f48d58e2b7f3fd9ca3ad1a3cc21f28d400a Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Tue, 24 Sep 2024 11:10:13 +0300 Subject: [PATCH 19/30] Bug fix: Change eps to EPS --- fieldsolver/derivatives.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fieldsolver/derivatives.cpp b/fieldsolver/derivatives.cpp index cdf19055e..be69314b7 100644 --- a/fieldsolver/derivatives.cpp +++ b/fieldsolver/derivatives.cpp @@ -836,7 +836,7 @@ void calculateScaledDeltas( 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) { + if (maxV > EPS) { amr_vorticity = vorticity * cell->parameters[CellParams::DX] / maxV; } From 3396cc6f6d4a737b3911a7a46bf394ad4c573928 Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Tue, 24 Sep 2024 16:54:01 +0300 Subject: [PATCH 20/30] Added use_anisotropy config parameter and fixed default values of vorticity and anisotropy params --- parameters.cpp | 17 +++++++++++++---- parameters.h | 1 + projects/project.cpp | 4 ++-- 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/parameters.cpp b/parameters.cpp index 88ee037d2..c16713987 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -167,10 +167,11 @@ Real P::alpha1CoarsenThreshold = -1.0; bool P::useAlpha2 = true; Real P::alpha2RefineThreshold = 0.5; Real P::alpha2CoarsenThreshold = -1.0; -bool P::useVorticity = true; +bool P::useVorticity = false; Real P::vorticityRefineThreshold = 0.5; Real P::vorticityCoarsenThreshold = -1.0; -Real P::anisotropyRefineThreshold = -1; +bool P::useAnisotropy = false; +Real P::anisotropyRefineThreshold = 0.5; int P::anisotropyMaxReflevel = 2; Real P::alphaDRhoWeight = 1.0; Real P::alphaDUWeight = 1.0; @@ -489,10 +490,11 @@ 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", true); + 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.anisotropy_refine_threshold","Determines the maximum value of pressure anisotropy to refine cells", -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_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); @@ -825,8 +827,15 @@ void Parameters::getParameters() { 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_max_reflevel", P::anisotropyMaxReflevel); + 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); diff --git a/parameters.h b/parameters.h index eb1534863..1fc379770 100644 --- a/parameters.h +++ b/parameters.h @@ -202,6 +202,7 @@ struct Parameters { static bool useVorticity; static Real vorticityRefineThreshold; static Real vorticityCoarsenThreshold; + static bool useAnisotropy; static Real anisotropyRefineThreshold; static int anisotropyMaxReflevel; static uint refineCadence; diff --git a/projects/project.cpp b/projects/project.cpp index fc3655826..4c14f50a8 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -612,7 +612,7 @@ namespace projects { }; // Pressure anisotropy forces refinement when under threshold - if (r2 < r_max2 && P::anisotropyRefineThreshold > 0 && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && refLevel < P::anisotropyMaxReflevel) { + if (r2 < r_max2 && P::useAnisotropy && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && refLevel < P::anisotropyMaxReflevel) { shouldRefine = true; shouldUnrefine = false; } @@ -649,7 +649,7 @@ namespace projects { (P::useAlpha1 ? mpiGrid[neighbor]->parameters[CellParams::AMR_ALPHA1] > P::alpha1RefineThreshold : false) || (P::useAlpha2 ? mpiGrid[neighbor]->parameters[CellParams::AMR_ALPHA2] > P::alpha2RefineThreshold : false) || (P::useVorticity ? mpiGrid[neighbor]->parameters[CellParams::AMR_VORTICITY] > P::vorticityRefineThreshold : false) || - (P::anisotropyRefineThreshold > 0 && mpiGrid[neighbor]->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && neighborRef < P::anisotropyMaxReflevel) + (P::useAnisotropy ? mpiGrid[neighbor]->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && neighborRef < P::anisotropyMaxReflevel : false) )}; if (shouldRefineNeighbor && // If the neighbor is planned to be refined, but is outside the allowed refinement region, cancel that refinement. From 1750e819f6a94cc3ec47740ff8c463a772288223 Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Mon, 30 Sep 2024 09:20:33 +0300 Subject: [PATCH 21/30] Make rotation correct and safer and clear offdiagonal cell params properly --- fieldsolver/derivatives.cpp | 5 +++-- vlasovsolver/cpu_moments.cpp | 9 +++++++++ 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/fieldsolver/derivatives.cpp b/fieldsolver/derivatives.cpp index be69314b7..5d0d96aeb 100644 --- a/fieldsolver/derivatives.cpp +++ b/fieldsolver/derivatives.cpp @@ -812,14 +812,15 @@ void calculateScaledDeltas( // 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{0, 0, 1}, Eigen::Vector3d{myB[0], myB[1], myB[2]}).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 Pprime = rot * P * rot.transpose(); + Eigen::Matrix3d transposerot = rot.transpose(); + Eigen::Matrix3d Pprime = rot * P * transposerot; Real Panisotropy {0.0}; if (Pprime(2, 2) > EPS) { diff --git a/vlasovsolver/cpu_moments.cpp b/vlasovsolver/cpu_moments.cpp index 184f50253..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 @@ -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); @@ -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); From 079d34ea50e7c17178820a70b692f35e8a1fe244 Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Wed, 9 Oct 2024 10:48:42 +0300 Subject: [PATCH 22/30] Fixed and refactored refinement and unrefinement logic --- parameters.cpp | 6 ++ parameters.h | 1 + projects/project.cpp | 157 ++++++++++++++++++++++--------------------- 3 files changed, 88 insertions(+), 76 deletions(-) diff --git a/parameters.cpp b/parameters.cpp index c16713987..9c91ffb2e 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -172,6 +172,7 @@ 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; @@ -495,6 +496,7 @@ bool P::addParameters() { 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); @@ -829,7 +831,11 @@ void Parameters::getParameters() { 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; diff --git a/parameters.h b/parameters.h index 1fc379770..bd63bdbad 100644 --- a/parameters.h +++ b/parameters.h @@ -204,6 +204,7 @@ struct Parameters { static Real vorticityCoarsenThreshold; static bool useAnisotropy; static Real anisotropyRefineThreshold; + static Real anisotropyCoarsenThreshold; static int anisotropyMaxReflevel; static uint refineCadence; static Real refineAfter; diff --git a/projects/project.cpp b/projects/project.cpp index 4c14f50a8..56654a31a 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -563,6 +563,81 @@ 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; + } + int Project::adaptRefinement( dccrg::Dccrg& mpiGrid ) const { phiprof::Timer refinesTimer {"Set refines"}; int myRank; @@ -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 @@ -596,86 +668,19 @@ namespace projects { // Cells too far from the ionosphere should be unrefined but // induced refinement still possible just beyond this r_max2 limit. // TODO starting to look truly cursed. Might need refactoring - 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) || - (P::useVorticity ? cell->parameters[CellParams::AMR_VORTICITY] > P::vorticityRefineThreshold : 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) && - (P::useVorticity ? cell->parameters[CellParams::AMR_VORTICITY] < P::vorticityCoarsenThreshold : true) - ) - }; - - // Pressure anisotropy forces refinement when under threshold - if (r2 < r_max2 && P::useAnisotropy && cell->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && refLevel < P::anisotropyMaxReflevel) { - shouldRefine = true; - shouldUnrefine = false; - } - 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) || - (P::useVorticity ? mpiGrid[neighbor]->parameters[CellParams::AMR_VORTICITY] > P::vorticityRefineThreshold : false) || - (P::useAnisotropy ? mpiGrid[neighbor]->parameters[CellParams::P_ANISOTROPY] < P::anisotropyRefineThreshold && neighborRef < P::anisotropyMaxReflevel : 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) && - (P::useVorticity ? mpiGrid[neighbor]->parameters[CellParams::AMR_VORTICITY] < P::vorticityCoarsenThreshold : 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; - } + shouldRefineNeighbor = shouldRefineCell(mpiGrid, neighbor, r_max2); + shouldUnrefineNeighbor = shouldUnrefineCell(mpiGrid, neighbor, r_max2); + int neighborRef {mpiGrid.get_refinement_level(neighbor)}; + if (neighborRef > refLevel && !shouldUnrefineNeighbor) { ++refined_neighbors; } else if (neighborRef < refLevel && !shouldRefineNeighbor) { From c17086460962e10aa89793e80fd4eea61767aa4b Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Wed, 9 Oct 2024 10:55:47 +0300 Subject: [PATCH 23/30] Add new functions to header --- projects/project.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/projects/project.h b/projects/project.h index 0e7f44dc3..ecf1fa38f 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. From ec57989145bd738962fcf18a78f432fb07ccbc65 Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Wed, 9 Oct 2024 11:01:48 +0300 Subject: [PATCH 24/30] Forgot to declare variables --- projects/project.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/projects/project.cpp b/projects/project.cpp index 56654a31a..20cc05cfb 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -677,8 +677,8 @@ namespace projects { int coarser_neighbors {0}; for (const auto& [neighbor, dir] : mpiGrid.get_face_neighbors_of(id)) { // Evaluate all face neighbors of the current cell - shouldRefineNeighbor = shouldRefineCell(mpiGrid, neighbor, r_max2); - shouldUnrefineNeighbor = shouldUnrefineCell(mpiGrid, neighbor, r_max2); + 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) { From 3d06814a03fe7794f75a4d66e149c7ed38a26af4 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Mon, 4 Nov 2024 16:48:35 +0200 Subject: [PATCH 25/30] Remove accidental commit --- warnings.txt | 34 ---------------------------------- 1 file changed, 34 deletions(-) delete mode 100644 warnings.txt diff --git a/warnings.txt b/warnings.txt deleted file mode 100644 index 6642ed85f..000000000 --- a/warnings.txt +++ /dev/null @@ -1,34 +0,0 @@ -./generate_version.sh: line 78: xxd: command not found -./generate_version.sh: line 141: xxd: command not found -make: *** [Makefile:235: datareducer.o] Interrupt -make: *** [Makefile:235: datareductionoperator.o] Interrupt -make: *** [Makefile:225: vamr_refinement_criteria.o] Interrupt -make: *** [Makefile:240: donotcompute.o] Interrupt -make: *** [Makefile:240: ionosphere.o] Interrupt -make: *** [Makefile:240: copysphere.o] Interrupt -make: *** [Makefile:240: outflow.o] Interrupt -make: *** [Makefile:240: inflow.o] Interrupt -make: *** [Makefile:240: setmaxwellian.o] Interrupt -make: *** [Makefile:245: fieldtracing.o] Interrupt -make: *** [Makefile:240: sysboundary.o] Interrupt -make: *** [Makefile:240: sysboundarycondition.o] Interrupt -make: *** [Makefile:250: project.o] Interrupt -make: *** [Makefile:250: projectTriAxisSearch.o] Interrupt -make: *** [Makefile:256: Alfven.o] Interrupt -make: *** [Makefile:256: Diffusion.o] Interrupt -make: *** [Makefile:256: Dispersion.o] Interrupt -make: *** [Makefile:256: Distributions.o] Interrupt -make: *** [Makefile:256: Firehose.o] Interrupt -make: *** [Makefile:256: Flowthrough.o] Interrupt -make: *** [Makefile:256: Fluctuations.o] Interrupt -make: *** [Makefile:256: Harris.o] Interrupt -make: *** [Makefile:256: KHB.o] Interrupt -make: *** [Makefile:256: Larmor.o] Interrupt -make: *** [Makefile:256: Magnetosphere.o] Interrupt -make: *** [Makefile:256: MultiPeak.o] Interrupt -make: *** [Makefile:256: Riemann1.o] Interrupt -make: *** [Makefile:256: Shock.o] Interrupt -make: *** [Makefile:256: Template.o] Interrupt -make: *** [Makefile:256: test_fp.o] Interrupt -make: *** [Makefile:256: testHall.o] Interrupt -make: *** [Makefile:256: test_trans.o] Interrupt From 3c4499982cee187070ccdbe558418983326f3118 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Mon, 4 Nov 2024 17:47:53 +0200 Subject: [PATCH 26/30] Multipop pressure anisotropy --- fieldsolver/derivatives.cpp | 54 +++++++++++++++++--------- spatial_cell_cpu.hpp | 17 +++++++-- vlasovsolver/cpu_moments.cpp | 73 +++++++++++++++++++----------------- vlasovsolver/vlasovmover.cpp | 6 ++- 4 files changed, 92 insertions(+), 58 deletions(-) diff --git a/fieldsolver/derivatives.cpp b/fieldsolver/derivatives.cpp index 520d7b3c8..101c51d0c 100644 --- a/fieldsolver/derivatives.cpp +++ b/fieldsolver/derivatives.cpp @@ -740,6 +740,32 @@ static Real calculateU(SpatialCell* cell) (rho > EPS ? (pow(p[0], 2) + pow(p[1], 2) + pow(p[2], 2)) / (2.0 * cell->parameters[CellParams::RHOM]) : 0.0); // Kinetic energy } +/*! \brief Calculates pressure anistotropy from B and P + * \param P elements of pressure order in order: P_11, P_22, P_33, P_23, P_13, P_12 + */ +static Real calculateAnistropy(const std::array& B, const std::array& P) +{ + // 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{B[0], B[1], B[2]}, Eigen::Vector3d{0, 0, 1}).normalized().toRotationMatrix(); + Eigen::Matrix3d Ptensor { + {P[0], P[5], P[4]}, + {P[5], P[1], P[3]}, + {P[4], P[3], P[2]}, + }; + + Eigen::Matrix3d transposerot = rot.transpose(); + Eigen::Matrix3d Pprime = rot * Ptensor * transposerot; + + Real Panisotropy {0.0}; + if (Pprime(2, 2) > EPS) { + Panisotropy = (Pprime(0, 0) + Pprime(1, 1)) / (2 * Pprime(2, 2)); + } + + return Panisotropy; +} + /*! \brief Low-level scaled gradients calculation * * For the SpatialCell* cell and its neighbors, calculate scaled gradients and their maximum alpha @@ -828,24 +854,6 @@ 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]}; @@ -860,6 +868,16 @@ void calculateScaledDeltas( amr_vorticity = vorticity * cell->parameters[CellParams::DX] / maxV; } + std::array myPressure {cell->parameters[CellParams::P_11], cell->parameters[CellParams::P_22], cell->parameters[CellParams::P_33], cell->parameters[CellParams::P_23], cell->parameters[CellParams::P_13], cell->parameters[CellParams::P_12]}; + Real Panisotropy {calculateAnistropy(myB, myPressure)}; + for (const auto& pop : cell->get_populations()) { + // TODO I hate this. Change all this crap to std::vectors? + std::array popP {pop.P[0], pop.P[1], pop.P[2], pop.P[3], pop.P[4], pop.P[5]}; + Real popPanisotropy {calculateAnistropy(myB, popP)}; + // low value refines + Panisotropy = std::min(Panisotropy, popPanisotropy); + } + cell->parameters[CellParams::AMR_DRHO] = dRho; cell->parameters[CellParams::AMR_DU] = dU; cell->parameters[CellParams::AMR_DPSQ] = dPsq; diff --git a/spatial_cell_cpu.hpp b/spatial_cell_cpu.hpp index d9ec4f5ee..4b917fb7a 100644 --- a/spatial_cell_cpu.hpp +++ b/spatial_cell_cpu.hpp @@ -156,9 +156,9 @@ namespace spatial_cell { Real V_R[3]; Real RHO_V; Real V_V[3]; - Real P[3]; - Real P_R[3]; - Real P_V[3]; + Real P[6]; + Real P_R[6]; + Real P_V[6]; Real RHOLOSSADJUST = 0.0; /*!< Counter for particle number loss from the destroying blocks in blockadjustment*/ Real max_dt[2]; /**< Element[0] is max_r_dt, element[1] max_v_dt.*/ Real velocityBlockMinValue; @@ -206,6 +206,9 @@ namespace spatial_cell { const Population & get_population(const uint popID) const; void set_population(const Population& pop, cuint popID); + std::vector& get_populations(); + const std::vector& get_populations() const; + uint8_t get_maximum_refinement_level(const uint popID); const Real& get_max_r_dt(const uint popID) const; const Real& get_max_v_dt(const uint popID) const; @@ -929,6 +932,14 @@ namespace spatial_cell { this->populations[popID] = pop; } + inline std::vector& SpatialCell::get_populations() { + return populations; + } + + inline const std::vector& SpatialCell::get_populations() const { + return populations; + } + inline const vmesh::LocalID* SpatialCell::get_velocity_grid_length(const uint popID,const uint8_t& refLevel) { return populations[popID].vmesh.getGridLength(refLevel); } diff --git a/vlasovsolver/cpu_moments.cpp b/vlasovsolver/cpu_moments.cpp index 7e032674a..b49b7c418 100644 --- a/vlasovsolver/cpu_moments.cpp +++ b/vlasovsolver/cpu_moments.cpp @@ -142,17 +142,17 @@ void calculateCellMoments(spatial_cell::SpatialCell* cell, } // Store species' contribution to bulk velocity moments - pop.P[0] = mass*array[0]; - pop.P[1] = mass*array[1]; - pop.P[2] = mass*array[2]; + for (size_t i = 0; i < array.size(); ++i) { + pop.P[i] = mass * array[i]; + } if (!computePopulationMomentsOnly) { - 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]; + 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_23] += pop.P[3]; + cell->parameters[CellParams::P_13] += pop.P[4]; + cell->parameters[CellParams::P_12] += pop.P[5]; } } // for-loop over particle species } @@ -199,9 +199,11 @@ void calculateMoments_R( Population & pop = cell->get_population(popID); if (blockContainer.size() == 0) { pop.RHO_R = 0; - for (int i=0; i<3; ++i) { - pop.V_R[i]=0; - pop.P_R[i]=0; + for (int i = 0; i < 3; ++i) { + pop.V_R[i] = 0; + } + for (int i = 0; i < 6; ++i) { + pop.P_R[i] = 0; } continue; } @@ -297,16 +299,16 @@ void calculateMoments_R( } // for-loop over velocity blocks // Store species' contribution to 2nd bulk velocity moments - pop.P_R[0] = mass*array[0]; - pop.P_R[1] = mass*array[1]; - pop.P_R[2] = mass*array[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 (size_t i = 0; i < array.size(); ++i) { + pop.P_R[i] = mass * array[i]; + } + + 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_23_R] += pop.P_R[3]; + cell->parameters[CellParams::P_13_R] += pop.P_R[4]; + cell->parameters[CellParams::P_12_R] += pop.P_R[5]; } // for-loop over spatial cells } // for-loop over particle species @@ -356,9 +358,11 @@ void calculateMoments_V( Population & pop = cell->get_population(popID); if (blockContainer.size() == 0) { pop.RHO_V = 0; - for (int i=0; i<3; ++i) { - pop.V_V[i]=0; - pop.P_V[i]=0; + for (int i = 0; i < 3; ++i) { + pop.V_V [i] = 0; + } + for (int i = 0; i < 6; ++i) { + pop.P_V [i] = 0; } continue; } @@ -441,17 +445,16 @@ void calculateMoments_V( } // for-loop over velocity blocks // Store species' contribution to 2nd bulk velocity moments - pop.P_V[0] = mass*array[0]; - pop.P_V[1] = mass*array[1]; - pop.P_V[2] = mass*array[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 (size_t i = 0; i < array.size(); ++i) { + pop.P_V[i] = mass * array[i]; + } + 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_23_V] += pop.P_V[3]; + cell->parameters[CellParams::P_13_V] += pop.P_V[4]; + cell->parameters[CellParams::P_12_V] += pop.P_V[5]; } // for-loop over spatial cells } // for-loop over particle species } diff --git a/vlasovsolver/vlasovmover.cpp b/vlasovsolver/vlasovmover.cpp index 62e290a9c..95db0a496 100644 --- a/vlasovsolver/vlasovmover.cpp +++ b/vlasovsolver/vlasovmover.cpp @@ -547,9 +547,11 @@ void calculateInterpolatedVelocityMoments( for (uint popID=0; popIDget_population(popID); pop.RHO = 0.5 * ( pop.RHO_R + pop.RHO_V ); - for(int i=0; i<3; i++) { + for (int i = 0; i < 3; i++) { pop.V[i] = 0.5 * ( pop.V_R[i] + pop.V_V[i] ); - pop.P[i] = 0.5 * ( pop.P_R[i] + pop.P_V[i] ); + } + for (int i = 0; i < 6; i++) { + pop.P[i] = 0.5 * ( pop.P_R[i] + pop.P_V[i] ); } } } From 10e5207b749f2aacdf62852c9783eae46bf5046a Mon Sep 17 00:00:00 2001 From: lkotipal Date: Mon, 4 Nov 2024 17:58:13 +0200 Subject: [PATCH 27/30] Revert "Multipop pressure anisotropy" This reverts commit 3c4499982cee187070ccdbe558418983326f3118. --- fieldsolver/derivatives.cpp | 54 +++++++++----------------- spatial_cell_cpu.hpp | 17 ++------- vlasovsolver/cpu_moments.cpp | 73 +++++++++++++++++------------------- vlasovsolver/vlasovmover.cpp | 6 +-- 4 files changed, 58 insertions(+), 92 deletions(-) diff --git a/fieldsolver/derivatives.cpp b/fieldsolver/derivatives.cpp index 101c51d0c..520d7b3c8 100644 --- a/fieldsolver/derivatives.cpp +++ b/fieldsolver/derivatives.cpp @@ -740,32 +740,6 @@ static Real calculateU(SpatialCell* cell) (rho > EPS ? (pow(p[0], 2) + pow(p[1], 2) + pow(p[2], 2)) / (2.0 * cell->parameters[CellParams::RHOM]) : 0.0); // Kinetic energy } -/*! \brief Calculates pressure anistotropy from B and P - * \param P elements of pressure order in order: P_11, P_22, P_33, P_23, P_13, P_12 - */ -static Real calculateAnistropy(const std::array& B, const std::array& P) -{ - // 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{B[0], B[1], B[2]}, Eigen::Vector3d{0, 0, 1}).normalized().toRotationMatrix(); - Eigen::Matrix3d Ptensor { - {P[0], P[5], P[4]}, - {P[5], P[1], P[3]}, - {P[4], P[3], P[2]}, - }; - - Eigen::Matrix3d transposerot = rot.transpose(); - Eigen::Matrix3d Pprime = rot * Ptensor * transposerot; - - Real Panisotropy {0.0}; - if (Pprime(2, 2) > EPS) { - Panisotropy = (Pprime(0, 0) + Pprime(1, 1)) / (2 * Pprime(2, 2)); - } - - return Panisotropy; -} - /*! \brief Low-level scaled gradients calculation * * For the SpatialCell* cell and its neighbors, calculate scaled gradients and their maximum alpha @@ -854,6 +828,24 @@ 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]}; @@ -868,16 +860,6 @@ void calculateScaledDeltas( amr_vorticity = vorticity * cell->parameters[CellParams::DX] / maxV; } - std::array myPressure {cell->parameters[CellParams::P_11], cell->parameters[CellParams::P_22], cell->parameters[CellParams::P_33], cell->parameters[CellParams::P_23], cell->parameters[CellParams::P_13], cell->parameters[CellParams::P_12]}; - Real Panisotropy {calculateAnistropy(myB, myPressure)}; - for (const auto& pop : cell->get_populations()) { - // TODO I hate this. Change all this crap to std::vectors? - std::array popP {pop.P[0], pop.P[1], pop.P[2], pop.P[3], pop.P[4], pop.P[5]}; - Real popPanisotropy {calculateAnistropy(myB, popP)}; - // low value refines - Panisotropy = std::min(Panisotropy, popPanisotropy); - } - cell->parameters[CellParams::AMR_DRHO] = dRho; cell->parameters[CellParams::AMR_DU] = dU; cell->parameters[CellParams::AMR_DPSQ] = dPsq; diff --git a/spatial_cell_cpu.hpp b/spatial_cell_cpu.hpp index 4b917fb7a..d9ec4f5ee 100644 --- a/spatial_cell_cpu.hpp +++ b/spatial_cell_cpu.hpp @@ -156,9 +156,9 @@ namespace spatial_cell { Real V_R[3]; Real RHO_V; Real V_V[3]; - Real P[6]; - Real P_R[6]; - Real P_V[6]; + Real P[3]; + Real P_R[3]; + Real P_V[3]; Real RHOLOSSADJUST = 0.0; /*!< Counter for particle number loss from the destroying blocks in blockadjustment*/ Real max_dt[2]; /**< Element[0] is max_r_dt, element[1] max_v_dt.*/ Real velocityBlockMinValue; @@ -206,9 +206,6 @@ namespace spatial_cell { const Population & get_population(const uint popID) const; void set_population(const Population& pop, cuint popID); - std::vector& get_populations(); - const std::vector& get_populations() const; - uint8_t get_maximum_refinement_level(const uint popID); const Real& get_max_r_dt(const uint popID) const; const Real& get_max_v_dt(const uint popID) const; @@ -932,14 +929,6 @@ namespace spatial_cell { this->populations[popID] = pop; } - inline std::vector& SpatialCell::get_populations() { - return populations; - } - - inline const std::vector& SpatialCell::get_populations() const { - return populations; - } - inline const vmesh::LocalID* SpatialCell::get_velocity_grid_length(const uint popID,const uint8_t& refLevel) { return populations[popID].vmesh.getGridLength(refLevel); } diff --git a/vlasovsolver/cpu_moments.cpp b/vlasovsolver/cpu_moments.cpp index b49b7c418..7e032674a 100644 --- a/vlasovsolver/cpu_moments.cpp +++ b/vlasovsolver/cpu_moments.cpp @@ -142,17 +142,17 @@ void calculateCellMoments(spatial_cell::SpatialCell* cell, } // Store species' contribution to bulk velocity moments - for (size_t i = 0; i < array.size(); ++i) { - pop.P[i] = mass * array[i]; - } + pop.P[0] = mass*array[0]; + pop.P[1] = mass*array[1]; + 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_23] += pop.P[3]; - cell->parameters[CellParams::P_13] += pop.P[4]; - cell->parameters[CellParams::P_12] += pop.P[5]; + 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 } @@ -199,11 +199,9 @@ void calculateMoments_R( Population & pop = cell->get_population(popID); if (blockContainer.size() == 0) { pop.RHO_R = 0; - for (int i = 0; i < 3; ++i) { - pop.V_R[i] = 0; - } - for (int i = 0; i < 6; ++i) { - pop.P_R[i] = 0; + for (int i=0; i<3; ++i) { + pop.V_R[i]=0; + pop.P_R[i]=0; } continue; } @@ -299,16 +297,16 @@ void calculateMoments_R( } // for-loop over velocity blocks // Store species' contribution to 2nd bulk velocity moments - for (size_t i = 0; i < array.size(); ++i) { - pop.P_R[i] = mass * array[i]; - } - - 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_23_R] += pop.P_R[3]; - cell->parameters[CellParams::P_13_R] += pop.P_R[4]; - cell->parameters[CellParams::P_12_R] += pop.P_R[5]; + pop.P_R[0] = mass*array[0]; + pop.P_R[1] = mass*array[1]; + pop.P_R[2] = mass*array[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 @@ -358,11 +356,9 @@ void calculateMoments_V( Population & pop = cell->get_population(popID); if (blockContainer.size() == 0) { pop.RHO_V = 0; - for (int i = 0; i < 3; ++i) { - pop.V_V [i] = 0; - } - for (int i = 0; i < 6; ++i) { - pop.P_V [i] = 0; + for (int i=0; i<3; ++i) { + pop.V_V[i]=0; + pop.P_V[i]=0; } continue; } @@ -445,16 +441,17 @@ void calculateMoments_V( } // for-loop over velocity blocks // Store species' contribution to 2nd bulk velocity moments - for (size_t i = 0; i < array.size(); ++i) { - pop.P_V[i] = mass * array[i]; - } + pop.P_V[0] = mass*array[0]; + pop.P_V[1] = mass*array[1]; + pop.P_V[2] = mass*array[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]; - 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_23_V] += pop.P_V[3]; - cell->parameters[CellParams::P_13_V] += pop.P_V[4]; - cell->parameters[CellParams::P_12_V] += pop.P_V[5]; } // for-loop over spatial cells } // for-loop over particle species } diff --git a/vlasovsolver/vlasovmover.cpp b/vlasovsolver/vlasovmover.cpp index 95db0a496..62e290a9c 100644 --- a/vlasovsolver/vlasovmover.cpp +++ b/vlasovsolver/vlasovmover.cpp @@ -547,11 +547,9 @@ void calculateInterpolatedVelocityMoments( for (uint popID=0; popIDget_population(popID); pop.RHO = 0.5 * ( pop.RHO_R + pop.RHO_V ); - for (int i = 0; i < 3; i++) { + for(int i=0; i<3; i++) { pop.V[i] = 0.5 * ( pop.V_R[i] + pop.V_V[i] ); - } - for (int i = 0; i < 6; i++) { - pop.P[i] = 0.5 * ( pop.P_R[i] + pop.P_V[i] ); + pop.P[i] = 0.5 * ( pop.P_R[i] + pop.P_V[i] ); } } } From a526c3c301deac9e69586e97e99f86db39ccd337 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Tue, 12 Nov 2024 10:29:08 +0200 Subject: [PATCH 28/30] Comment --- fieldsolver/gridGlue.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/fieldsolver/gridGlue.cpp b/fieldsolver/gridGlue.cpp index 624336420..4574d1e67 100644 --- a/fieldsolver/gridGlue.cpp +++ b/fieldsolver/gridGlue.cpp @@ -314,6 +314,7 @@ void getFieldsFromFsGrid( 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); From b36f471ea78287263d187f50447607cf7cb24200 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Tue, 12 Nov 2024 10:35:31 +0200 Subject: [PATCH 29/30] Early return logic --- projects/project.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/projects/project.cpp b/projects/project.cpp index 20cc05cfb..3f70ed52a 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -644,7 +644,7 @@ namespace projects { MPI_Comm_rank(MPI_COMM_WORLD,&myRank); int 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; } From fda23d91532a5c7f18ca87fed55abe6213067789 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Tue, 12 Nov 2024 10:38:28 +0200 Subject: [PATCH 30/30] Outdated comment --- projects/project.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/projects/project.cpp b/projects/project.cpp index 3f70ed52a..801709ea6 100644 --- a/projects/project.cpp +++ b/projects/project.cpp @@ -665,10 +665,6 @@ 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. - // TODO starting to look truly cursed. Might need refactoring - bool shouldRefine = shouldRefineCell(mpiGrid, id, r_max2); bool shouldUnrefine = shouldUnrefineCell(mpiGrid, id, r_max2);