From f132ec6a02c1f36c21cd4dde6204ee926cf27bc2 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Tue, 12 Nov 2024 13:55:54 +0200 Subject: [PATCH] Revert "Revert "Multipop pressure anisotropy"" This reverts commit 10e5207b749f2aacdf62852c9783eae46bf5046a. --- 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] ); } } }