From d9be0e80f3d7293b2adb56290ff7f699412c6264 Mon Sep 17 00:00:00 2001 From: ykempf Date: Thu, 21 Sep 2023 16:53:22 +0300 Subject: [PATCH 1/3] Dummy commit --- common.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/common.h b/common.h index 99581af0e..8d2bfe48d 100644 --- a/common.h +++ b/common.h @@ -1,4 +1,4 @@ -/* +/* * This file is part of Vlasiator. * Copyright 2010-2016 Finnish Meteorological Institute * From 5f7118c4f790dab48c4c2158fe5b33b9963da214 Mon Sep 17 00:00:00 2001 From: ykempf Date: Tue, 5 Dec 2023 09:45:38 +0200 Subject: [PATCH 2/3] Suggested changes to clang-format rules. --- .clang-format | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.clang-format b/.clang-format index 743a8ef78..f9800f199 100644 --- a/.clang-format +++ b/.clang-format @@ -10,5 +10,9 @@ AllowShortIfStatementsOnASingleLine: false IndentCaseLabels: false ColumnLimit: 120 AccessModifierOffset: -3 +BinPackArguments: false +BinPackParameters: false +ExperimentalAutoDetectBinPacking: false +AllowAllParametersOfDeclarationOnNextLine: false ... From b75ed38e6bd4cd84eece92ef0d9f0907e03252b1 Mon Sep 17 00:00:00 2001 From: ykempf Date: Tue, 5 Dec 2023 10:31:34 +0200 Subject: [PATCH 3/3] clang-format on fieldsolver ldz_electric_field --- fieldsolver/ldz_electric_field.cpp | 2508 +++++++++++++++------------- fieldsolver/ldz_electric_field.hpp | 29 +- 2 files changed, 1389 insertions(+), 1148 deletions(-) diff --git a/fieldsolver/ldz_electric_field.cpp b/fieldsolver/ldz_electric_field.cpp index d8a30c490..3a7da5977 100644 --- a/fieldsolver/ldz_electric_field.cpp +++ b/fieldsolver/ldz_electric_field.cpp @@ -26,7 +26,7 @@ #include "ldz_electric_field.hpp" #ifndef NDEBUG - #define DEBUG_FSOLVER +#define DEBUG_FSOLVER #endif namespace pc = physicalconstants; @@ -37,10 +37,10 @@ using namespace std; * Computes the correct combination of speeds to determine the CFL limits. * * It should be in-plane, but we use the complete wave speeds from calculateWaveSpeed??(). - * + * * At the moment it computes the geometric mean of both bulk velocity components * and takes the maximum of that plus either the magnetosonic or the whistler speed. - * + * * \sa calculateWaveSpeedYZ calculateWaveSpeedXY calculateWaveSpeedZX * * \param v0 Flow in first direction @@ -49,28 +49,22 @@ using namespace std; * \param vS Sound speed * \param vW Whistler speed */ -Real calculateCflSpeed( - const Real& v0, - const Real& v1, - const Real& vA, - const Real& vS, - const Real& vW -) { - const Real v = sqrt(v0*v0 + v1*v1); - const Real vMS = sqrt(vA*vA + vS*vS); +Real calculateCflSpeed(const Real& v0, const Real& v1, const Real& vA, const Real& vS, const Real& vW) { + const Real v = sqrt(v0 * v0 + v1 * v1); + const Real vMS = sqrt(vA * vA + vS * vS); return max(v + vMS, v + vW); } /*! \brief Low-level helper function. - * + * * Computes the magnetosonic speed in the YZ plane. Used in upwinding the electric field X component, * at the interface between cell (i,j,k) and (nbi,nbj,nbk). - * + * * Expects that the correct RHO and B fields are being passed, depending on the * stage of the Runge-Kutta time stepping method. - * + * * If fields are not propagated, returns 0.0 as there is no information propagating. - * + * * \param perBGrid fsGrid holding the perturbed B quantities * \param momentsGrid fsGrid holding the moment quantities * \param dPerBGrid fsGrid holding the derivatives of perturbed B @@ -92,67 +86,76 @@ Real calculateCflSpeed( * \param ret_vS Sound speed returned * \param ret_vW Whistler speed returned */ -void calculateWaveSpeedYZ( - FsGrid< std::array, FS_STENCIL_WIDTH> & perBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, - cint i, - cint j, - cint k, - cint nbi, - cint nbj, - cint nbk, - const Real& By, - const Real& Bz, - const Real& dBydx, - const Real& dBydz, - const Real& dBzdx, - const Real& dBzdy, - const Real& ydir, - const Real& zdir, - const Real& minRhom, - const Real& maxRhom, - Real& ret_vA, - Real& ret_vS, - Real& ret_vW -) { - std::array * perb = perBGrid.get(i,j,k); - std::array * nbr_perb = perBGrid.get(nbi,nbj,nbk); - std::array * moments = momentsGrid.get(i,j,k); - std::array * dmoments = dMomentsGrid.get(i,j,k); - std::array * dperb = dPerBGrid.get(i,j,k); - std::array * nbr_dperb = dPerBGrid.get(nbi,nbj,nbk); - std::array * bgb = BgBGrid.get(i,j,k); - std::array * nbr_bgb = BgBGrid.get(nbi,nbj,nbk); - +void calculateWaveSpeedYZ(FsGrid, FS_STENCIL_WIDTH>& perBGrid, + FsGrid, FS_STENCIL_WIDTH>& momentsGrid, + FsGrid, FS_STENCIL_WIDTH>& dPerBGrid, + FsGrid, FS_STENCIL_WIDTH>& dMomentsGrid, + FsGrid, FS_STENCIL_WIDTH>& BgBGrid, + cint i, + cint j, + cint k, + cint nbi, + cint nbj, + cint nbk, + const Real& By, + const Real& Bz, + const Real& dBydx, + const Real& dBydz, + const Real& dBzdx, + const Real& dBzdy, + const Real& ydir, + const Real& zdir, + const Real& minRhom, + const Real& maxRhom, + Real& ret_vA, + Real& ret_vS, + Real& ret_vW) { + std::array* perb = perBGrid.get(i, j, k); + std::array* nbr_perb = perBGrid.get(nbi, nbj, nbk); + std::array* moments = momentsGrid.get(i, j, k); + std::array* dmoments = dMomentsGrid.get(i, j, k); + std::array* dperb = dPerBGrid.get(i, j, k); + std::array* nbr_dperb = dPerBGrid.get(nbi, nbj, nbk); + std::array* bgb = BgBGrid.get(i, j, k); + std::array* nbr_bgb = BgBGrid.get(nbi, nbj, nbk); + Real A_0, A_X, rhom, p11, p22, p33; - A_0 = HALF*(nbr_perb->at(fsgrids::bfield::PERBX) + nbr_bgb->at(fsgrids::bgbfield::BGBX) + perb->at(fsgrids::bfield::PERBX) + bgb->at(fsgrids::bgbfield::BGBX)); - A_X = (nbr_perb->at(fsgrids::bfield::PERBX) + nbr_bgb->at(fsgrids::bgbfield::BGBX)) - (perb->at(fsgrids::bfield::PERBX) + bgb->at(fsgrids::bgbfield::BGBX)); - rhom = moments->at(fsgrids::moments::RHOM) + ydir*HALF*dmoments->at(fsgrids::dmoments::drhomdy) + zdir*HALF*dmoments->at(fsgrids::dmoments::drhomdz); - p11 = moments->at(fsgrids::moments::P_11) + ydir*HALF*dmoments->at(fsgrids::dmoments::dp11dy) + zdir*HALF*dmoments->at(fsgrids::dmoments::dp11dz); - p22 = moments->at(fsgrids::moments::P_22) + ydir*HALF*dmoments->at(fsgrids::dmoments::dp22dy) + zdir*HALF*dmoments->at(fsgrids::dmoments::dp22dz); - p33 = moments->at(fsgrids::moments::P_33) + ydir*HALF*dmoments->at(fsgrids::dmoments::dp33dy) + zdir*HALF*dmoments->at(fsgrids::dmoments::dp33dz); - + A_0 = HALF * (nbr_perb->at(fsgrids::bfield::PERBX) + nbr_bgb->at(fsgrids::bgbfield::BGBX) + + perb->at(fsgrids::bfield::PERBX) + bgb->at(fsgrids::bgbfield::BGBX)); + A_X = (nbr_perb->at(fsgrids::bfield::PERBX) + nbr_bgb->at(fsgrids::bgbfield::BGBX)) - + (perb->at(fsgrids::bfield::PERBX) + bgb->at(fsgrids::bgbfield::BGBX)); + rhom = moments->at(fsgrids::moments::RHOM) + ydir * HALF * dmoments->at(fsgrids::dmoments::drhomdy) + + zdir * HALF * dmoments->at(fsgrids::dmoments::drhomdz); + p11 = moments->at(fsgrids::moments::P_11) + ydir * HALF * dmoments->at(fsgrids::dmoments::dp11dy) + + zdir * HALF * dmoments->at(fsgrids::dmoments::dp11dz); + p22 = moments->at(fsgrids::moments::P_22) + ydir * HALF * dmoments->at(fsgrids::dmoments::dp22dy) + + zdir * HALF * dmoments->at(fsgrids::dmoments::dp22dz); + p33 = moments->at(fsgrids::moments::P_33) + ydir * HALF * dmoments->at(fsgrids::dmoments::dp33dy) + + zdir * HALF * dmoments->at(fsgrids::dmoments::dp33dz); + if (rhom < minRhom) { rhom = minRhom; } else if (rhom > maxRhom) { rhom = maxRhom; } - - const Real A_Y = nbr_dperb->at(fsgrids::dperb::dPERBxdy) + nbr_bgb->at(fsgrids::bgbfield::dBGBxdy) + dperb->at(fsgrids::dperb::dPERBxdy) + bgb->at(fsgrids::bgbfield::dBGBxdy); - const Real A_XY = nbr_dperb->at(fsgrids::dperb::dPERBxdy) + nbr_bgb->at(fsgrids::bgbfield::dBGBxdy) - (dperb->at(fsgrids::dperb::dPERBxdy) + bgb->at(fsgrids::bgbfield::dBGBxdy)); - const Real A_Z = nbr_dperb->at(fsgrids::dperb::dPERBxdz) + nbr_bgb->at(fsgrids::bgbfield::dBGBxdz) + dperb->at(fsgrids::dperb::dPERBxdz) + bgb->at(fsgrids::bgbfield::dBGBxdz); - const Real A_XZ = nbr_dperb->at(fsgrids::dperb::dPERBxdz) + nbr_bgb->at(fsgrids::bgbfield::dBGBxdz) - (dperb->at(fsgrids::dperb::dPERBxdz) + bgb->at(fsgrids::bgbfield::dBGBxdz)); - - const Real Bx2 = (A_0 + ydir*HALF*A_Y + zdir*HALF*A_Z)*(A_0 + ydir*HALF*A_Y + zdir*HALF*A_Z) - + TWELWTH*(A_X + ydir*HALF*A_XY + zdir*HALF*A_XZ)*(A_X + ydir*HALF*A_XY + zdir*HALF*A_XZ); // OK - const Real By2 = (By + zdir*HALF*dBydz)*(By + zdir*HALF*dBydz) + TWELWTH*dBydx*dBydx; // OK - const Real Bz2 = (Bz + ydir*HALF*dBzdy)*(Bz + ydir*HALF*dBzdy) + TWELWTH*dBzdx*dBzdx; // OK - + + const Real A_Y = nbr_dperb->at(fsgrids::dperb::dPERBxdy) + nbr_bgb->at(fsgrids::bgbfield::dBGBxdy) + + dperb->at(fsgrids::dperb::dPERBxdy) + bgb->at(fsgrids::bgbfield::dBGBxdy); + const Real A_XY = nbr_dperb->at(fsgrids::dperb::dPERBxdy) + nbr_bgb->at(fsgrids::bgbfield::dBGBxdy) - + (dperb->at(fsgrids::dperb::dPERBxdy) + bgb->at(fsgrids::bgbfield::dBGBxdy)); + const Real A_Z = nbr_dperb->at(fsgrids::dperb::dPERBxdz) + nbr_bgb->at(fsgrids::bgbfield::dBGBxdz) + + dperb->at(fsgrids::dperb::dPERBxdz) + bgb->at(fsgrids::bgbfield::dBGBxdz); + const Real A_XZ = nbr_dperb->at(fsgrids::dperb::dPERBxdz) + nbr_bgb->at(fsgrids::bgbfield::dBGBxdz) - + (dperb->at(fsgrids::dperb::dPERBxdz) + bgb->at(fsgrids::bgbfield::dBGBxdz)); + + const Real Bx2 = (A_0 + ydir * HALF * A_Y + zdir * HALF * A_Z) * (A_0 + ydir * HALF * A_Y + zdir * HALF * A_Z) + + TWELWTH * (A_X + ydir * HALF * A_XY + zdir * HALF * A_XZ) * + (A_X + ydir * HALF * A_XY + zdir * HALF * A_XZ); // OK + const Real By2 = (By + zdir * HALF * dBydz) * (By + zdir * HALF * dBydz) + TWELWTH * dBydx * dBydx; // OK + const Real Bz2 = (Bz + ydir * HALF * dBzdy) * (Bz + ydir * HALF * dBzdy) + TWELWTH * dBzdx * dBzdx; // OK + const Real Bmag2 = Bx2 + By2 + Bz2; - + p11 = p11 < 0.0 ? 0.0 : p11; p22 = p22 < 0.0 ? 0.0 : p22; p33 = p33 < 0.0 ? 0.0 : p33; @@ -172,28 +175,35 @@ void calculateWaveSpeedYZ( // and // http://iopscience.iop.org/article/10.1088/0253-6102/43/2/026/meta (Alfven waves) // for details. - const Real vA2 = divideIfNonZero(Bmag2, pc::MU_0*rhom); // Alfven speed - const Real vS2 = divideIfNonZero(p11+p22+p33, 2.0*rhom); // sound speed, adiabatic coefficient 3/2, P=1/3*trace in sound speed -// const Real vW = Parameters::ohmHallTerm > 0 ? divideIfNonZero(2.0*M_PI*vA2*pc::MASS_PROTON, perBGrid.DX*pc::CHARGE*sqrt(Bmag2)) : 0.0; // whistler speed - const Real vW = Parameters::ohmHallTerm > 0 ? - sqrt(vA2) * (1 + divideIfNonZero(2*M_PI*M_PI*pc::MASS_PROTON*pc::MASS_PROTON, perBGrid.DX*perBGrid.DX*rhom*pc::CHARGE*pc::CHARGE*pc::MU_0) - / sqrt(1 + divideIfNonZero( M_PI*M_PI*pc::MASS_PROTON*pc::MASS_PROTON, perBGrid.DX*perBGrid.DX*rhom*pc::CHARGE*pc::CHARGE*pc::MU_0))) - : 0.0; // whistler speed - + const Real vA2 = divideIfNonZero(Bmag2, pc::MU_0 * rhom); // Alfven speed + const Real vS2 = divideIfNonZero(p11 + p22 + p33, + 2.0 * rhom); // sound speed, adiabatic coefficient 3/2, P=1/3*trace in sound speed + // const Real vW = Parameters::ohmHallTerm > 0 ? divideIfNonZero(2.0*M_PI*vA2*pc::MASS_PROTON, + // perBGrid.DX*pc::CHARGE*sqrt(Bmag2)) : 0.0; // whistler speed + const Real vW = + Parameters::ohmHallTerm > 0 + ? sqrt(vA2) * (1 + divideIfNonZero(2 * M_PI * M_PI * pc::MASS_PROTON * pc::MASS_PROTON, + perBGrid.DX * perBGrid.DX * rhom * pc::CHARGE * pc::CHARGE * pc::MU_0) / + sqrt(1 + divideIfNonZero(M_PI * M_PI * pc::MASS_PROTON * pc::MASS_PROTON, + perBGrid.DX * perBGrid.DX * rhom * pc::CHARGE * pc::CHARGE * + pc::MU_0))) + : 0.0; // whistler speed + ret_vA = sqrt(vA2); ret_vS = sqrt(vS2); ret_vW = vW; } /*! \brief Low-level helper function. - * + * * Computes the magnetosonic speed in the XZ plane. Used in upwinding the electric field Y component, * at the interface between cell (i,j,k) and (nbi,nbj,nbk). - * - * Expects that the correct RHO and B fields are being passed, depending on the stage of the Runge-Kutta time stepping method. - * + * + * Expects that the correct RHO and B fields are being passed, depending on the stage of the Runge-Kutta time stepping + * method. + * * If fields are not propagated, returns 0.0 as there is no information propagating. - * + * * \param perBGrid fsGrid holding the perturbed B quantities * \param momentsGrid fsGrid holding the moment quantities * \param dPerBGrid fsGrid holding the derivatives of perturbed B @@ -215,71 +225,80 @@ void calculateWaveSpeedYZ( * \param ret_vS Sound speed returned * \param ret_vW Whistler speed returned */ -void calculateWaveSpeedXZ( - FsGrid< std::array, FS_STENCIL_WIDTH> & perBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, - cint i, - cint j, - cint k, - cint nbi, - cint nbj, - cint nbk, - const Real& Bx, - const Real& Bz, - const Real& dBxdy, - const Real& dBxdz, - const Real& dBzdx, - const Real& dBzdy, - const Real& xdir, - const Real& zdir, - const Real& minRhom, - const Real& maxRhom, - Real& ret_vA, - Real& ret_vS, - Real& ret_vW -) { - std::array * perb = perBGrid.get(i,j,k); - std::array * nbr_perb = perBGrid.get(nbi,nbj,nbk); - std::array * moments = momentsGrid.get(i,j,k); - std::array * dmoments = dMomentsGrid.get(i,j,k); - std::array * dperb = dPerBGrid.get(i,j,k); - std::array * nbr_dperb = dPerBGrid.get(nbi,nbj,nbk); - std::array * bgb = BgBGrid.get(i,j,k); - std::array * nbr_bgb = BgBGrid.get(nbi,nbj,nbk); - +void calculateWaveSpeedXZ(FsGrid, FS_STENCIL_WIDTH>& perBGrid, + FsGrid, FS_STENCIL_WIDTH>& momentsGrid, + FsGrid, FS_STENCIL_WIDTH>& dPerBGrid, + FsGrid, FS_STENCIL_WIDTH>& dMomentsGrid, + FsGrid, FS_STENCIL_WIDTH>& BgBGrid, + cint i, + cint j, + cint k, + cint nbi, + cint nbj, + cint nbk, + const Real& Bx, + const Real& Bz, + const Real& dBxdy, + const Real& dBxdz, + const Real& dBzdx, + const Real& dBzdy, + const Real& xdir, + const Real& zdir, + const Real& minRhom, + const Real& maxRhom, + Real& ret_vA, + Real& ret_vS, + Real& ret_vW) { + std::array* perb = perBGrid.get(i, j, k); + std::array* nbr_perb = perBGrid.get(nbi, nbj, nbk); + std::array* moments = momentsGrid.get(i, j, k); + std::array* dmoments = dMomentsGrid.get(i, j, k); + std::array* dperb = dPerBGrid.get(i, j, k); + std::array* nbr_dperb = dPerBGrid.get(nbi, nbj, nbk); + std::array* bgb = BgBGrid.get(i, j, k); + std::array* nbr_bgb = BgBGrid.get(nbi, nbj, nbk); + Real B_0, B_Y, rhom, p11, p22, p33; - B_0 = HALF*(nbr_perb->at(fsgrids::bfield::PERBY) + nbr_bgb->at(fsgrids::bgbfield::BGBY) + perb->at(fsgrids::bfield::PERBY) + bgb->at(fsgrids::bgbfield::BGBY)); - B_Y = (nbr_perb->at(fsgrids::bfield::PERBY) + nbr_bgb->at(fsgrids::bgbfield::BGBY)) - (perb->at(fsgrids::bfield::PERBY) + bgb->at(fsgrids::bgbfield::BGBY)); - rhom = moments->at(fsgrids::moments::RHOM) + xdir*HALF*dmoments->at(fsgrids::dmoments::drhomdx) + zdir*HALF*dmoments->at(fsgrids::dmoments::drhomdz); - p11 = moments->at(fsgrids::moments::P_11) + xdir*HALF*dmoments->at(fsgrids::dmoments::dp11dx) + zdir*HALF*dmoments->at(fsgrids::dmoments::dp11dz); - p22 = moments->at(fsgrids::moments::P_22) + xdir*HALF*dmoments->at(fsgrids::dmoments::dp22dx) + zdir*HALF*dmoments->at(fsgrids::dmoments::dp22dz); - p33 = moments->at(fsgrids::moments::P_33) + xdir*HALF*dmoments->at(fsgrids::dmoments::dp33dx) + zdir*HALF*dmoments->at(fsgrids::dmoments::dp33dz); - + B_0 = HALF * (nbr_perb->at(fsgrids::bfield::PERBY) + nbr_bgb->at(fsgrids::bgbfield::BGBY) + + perb->at(fsgrids::bfield::PERBY) + bgb->at(fsgrids::bgbfield::BGBY)); + B_Y = (nbr_perb->at(fsgrids::bfield::PERBY) + nbr_bgb->at(fsgrids::bgbfield::BGBY)) - + (perb->at(fsgrids::bfield::PERBY) + bgb->at(fsgrids::bgbfield::BGBY)); + rhom = moments->at(fsgrids::moments::RHOM) + xdir * HALF * dmoments->at(fsgrids::dmoments::drhomdx) + + zdir * HALF * dmoments->at(fsgrids::dmoments::drhomdz); + p11 = moments->at(fsgrids::moments::P_11) + xdir * HALF * dmoments->at(fsgrids::dmoments::dp11dx) + + zdir * HALF * dmoments->at(fsgrids::dmoments::dp11dz); + p22 = moments->at(fsgrids::moments::P_22) + xdir * HALF * dmoments->at(fsgrids::dmoments::dp22dx) + + zdir * HALF * dmoments->at(fsgrids::dmoments::dp22dz); + p33 = moments->at(fsgrids::moments::P_33) + xdir * HALF * dmoments->at(fsgrids::dmoments::dp33dx) + + zdir * HALF * dmoments->at(fsgrids::dmoments::dp33dz); + if (rhom < minRhom) { rhom = minRhom; } else if (rhom > maxRhom) { rhom = maxRhom; } - - const Real B_X = nbr_dperb->at(fsgrids::dperb::dPERBydx) + nbr_bgb->at(fsgrids::bgbfield::dBGBydx) + dperb->at(fsgrids::dperb::dPERBydx) + bgb->at(fsgrids::bgbfield::dBGBydx); - const Real B_XY = nbr_dperb->at(fsgrids::dperb::dPERBydx) + nbr_bgb->at(fsgrids::bgbfield::dBGBydx) - (dperb->at(fsgrids::dperb::dPERBydx) + bgb->at(fsgrids::bgbfield::dBGBydx)); - const Real B_Z = nbr_dperb->at(fsgrids::dperb::dPERBydz) + nbr_bgb->at(fsgrids::bgbfield::dBGBydz) + dperb->at(fsgrids::dperb::dPERBydz) + bgb->at(fsgrids::bgbfield::dBGBydz); - const Real B_YZ = nbr_dperb->at(fsgrids::dperb::dPERBydz) + nbr_bgb->at(fsgrids::bgbfield::dBGBydz) - (dperb->at(fsgrids::dperb::dPERBydz) + bgb->at(fsgrids::bgbfield::dBGBydz)); - - const Real By2 = (B_0 + xdir*HALF*B_X + zdir*HALF*B_Z)*(B_0 + xdir*HALF*B_X + zdir*HALF*B_Z) - + TWELWTH*(B_Y + xdir*HALF*B_XY + zdir*HALF*B_YZ)*(B_Y + xdir*HALF*B_XY + zdir*HALF*B_YZ); // OK - const Real Bx2 = (Bx + zdir*HALF*dBxdz)*(Bx + zdir*HALF*dBxdz) + TWELWTH*dBxdy*dBxdy; // OK - const Real Bz2 = (Bz + xdir*HALF*dBzdx)*(Bz + xdir*HALF*dBzdx) + TWELWTH*dBzdy*dBzdy; // OK - + + const Real B_X = nbr_dperb->at(fsgrids::dperb::dPERBydx) + nbr_bgb->at(fsgrids::bgbfield::dBGBydx) + + dperb->at(fsgrids::dperb::dPERBydx) + bgb->at(fsgrids::bgbfield::dBGBydx); + const Real B_XY = nbr_dperb->at(fsgrids::dperb::dPERBydx) + nbr_bgb->at(fsgrids::bgbfield::dBGBydx) - + (dperb->at(fsgrids::dperb::dPERBydx) + bgb->at(fsgrids::bgbfield::dBGBydx)); + const Real B_Z = nbr_dperb->at(fsgrids::dperb::dPERBydz) + nbr_bgb->at(fsgrids::bgbfield::dBGBydz) + + dperb->at(fsgrids::dperb::dPERBydz) + bgb->at(fsgrids::bgbfield::dBGBydz); + const Real B_YZ = nbr_dperb->at(fsgrids::dperb::dPERBydz) + nbr_bgb->at(fsgrids::bgbfield::dBGBydz) - + (dperb->at(fsgrids::dperb::dPERBydz) + bgb->at(fsgrids::bgbfield::dBGBydz)); + + const Real By2 = (B_0 + xdir * HALF * B_X + zdir * HALF * B_Z) * (B_0 + xdir * HALF * B_X + zdir * HALF * B_Z) + + TWELWTH * (B_Y + xdir * HALF * B_XY + zdir * HALF * B_YZ) * + (B_Y + xdir * HALF * B_XY + zdir * HALF * B_YZ); // OK + const Real Bx2 = (Bx + zdir * HALF * dBxdz) * (Bx + zdir * HALF * dBxdz) + TWELWTH * dBxdy * dBxdy; // OK + const Real Bz2 = (Bz + xdir * HALF * dBzdx) * (Bz + xdir * HALF * dBzdx) + TWELWTH * dBzdy * dBzdy; // OK + const Real Bmag2 = Bx2 + By2 + Bz2; - + p11 = p11 < 0.0 ? 0.0 : p11; p22 = p22 < 0.0 ? 0.0 : p22; p33 = p33 < 0.0 ? 0.0 : p33; - + // Effective wave speeds for advection and CFL calculation // Note that these are calculated as if the plasma is purely made up of hydrogen, which // is a reasonable approximation if it is proton-dominant. @@ -295,28 +314,35 @@ void calculateWaveSpeedXZ( // and // http://iopscience.iop.org/article/10.1088/0253-6102/43/2/026/meta (Alfven waves) // for details. - const Real vA2 = divideIfNonZero(Bmag2, pc::MU_0*rhom); // Alfven speed - const Real vS2 = divideIfNonZero(p11+p22+p33, 2.0*rhom); // sound speed, adiabatic coefficient 3/2, P=1/3*trace in sound speed -// const Real vW = Parameters::ohmHallTerm > 0 ? divideIfNonZero(2.0*M_PI*vA2*pc::MASS_PROTON, perBGrid.DX*pc::CHARGE*sqrt(Bmag2)) : 0.0; // whistler speed - const Real vW = Parameters::ohmHallTerm > 0 ? - sqrt(vA2) * (1 + divideIfNonZero(2*M_PI*M_PI*pc::MASS_PROTON*pc::MASS_PROTON, perBGrid.DX*perBGrid.DX*rhom*pc::CHARGE*pc::CHARGE*pc::MU_0) - / sqrt(1 + divideIfNonZero( M_PI*M_PI*pc::MASS_PROTON*pc::MASS_PROTON, perBGrid.DX*perBGrid.DX*rhom*pc::CHARGE*pc::CHARGE*pc::MU_0))) - : 0.0; // whistler speed - + const Real vA2 = divideIfNonZero(Bmag2, pc::MU_0 * rhom); // Alfven speed + const Real vS2 = divideIfNonZero(p11 + p22 + p33, + 2.0 * rhom); // sound speed, adiabatic coefficient 3/2, P=1/3*trace in sound speed + // const Real vW = Parameters::ohmHallTerm > 0 ? divideIfNonZero(2.0*M_PI*vA2*pc::MASS_PROTON, + // perBGrid.DX*pc::CHARGE*sqrt(Bmag2)) : 0.0; // whistler speed + const Real vW = + Parameters::ohmHallTerm > 0 + ? sqrt(vA2) * (1 + divideIfNonZero(2 * M_PI * M_PI * pc::MASS_PROTON * pc::MASS_PROTON, + perBGrid.DX * perBGrid.DX * rhom * pc::CHARGE * pc::CHARGE * pc::MU_0) / + sqrt(1 + divideIfNonZero(M_PI * M_PI * pc::MASS_PROTON * pc::MASS_PROTON, + perBGrid.DX * perBGrid.DX * rhom * pc::CHARGE * pc::CHARGE * + pc::MU_0))) + : 0.0; // whistler speed + ret_vA = sqrt(vA2); ret_vS = sqrt(vS2); ret_vW = vW; } /*! \brief Low-level helper function. - * + * * Computes the magnetosonic speed in the XY plane. Used in upwinding the electric field Z component, * at the interface between cell (i,j,k) and (nbi,nbj,nbk). - * - * Expects that the correct RHO and B fields are being passed, depending on the stage of the Runge-Kutta time stepping method. - * + * + * Expects that the correct RHO and B fields are being passed, depending on the stage of the Runge-Kutta time stepping + * method. + * * If fields are not propagated, returns 0.0 as there is no information propagating. - * + * * \param perBGrid fsGrid holding the perturbed B quantities * \param momentsGrid fsGrid holding the moment quantities * \param dPerBGrid fsGrid holding the derivatives of perturbed B @@ -338,71 +364,80 @@ void calculateWaveSpeedXZ( * \param ret_vS Sound speed returned * \param ret_vW Whistler speed returned */ -void calculateWaveSpeedXY( - FsGrid< std::array, FS_STENCIL_WIDTH> & perBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, - cint i, - cint j, - cint k, - cint nbi, - cint nbj, - cint nbk, - const Real& Bx, - const Real& By, - const Real& dBxdy, - const Real& dBxdz, - const Real& dBydx, - const Real& dBydz, - const Real& xdir, - const Real& ydir, - const Real& minRhom, - const Real& maxRhom, - Real& ret_vA, - Real& ret_vS, - Real& ret_vW -) { - std::array * perb = perBGrid.get(i,j,k); - std::array * nbr_perb = perBGrid.get(nbi,nbj,nbk); - std::array * moments = momentsGrid.get(i,j,k); - std::array * dmoments = dMomentsGrid.get(i,j,k); - std::array * dperb = dPerBGrid.get(i,j,k); - std::array * nbr_dperb = dPerBGrid.get(nbi,nbj,nbk); - std::array * bgb = BgBGrid.get(i,j,k); - std::array * nbr_bgb = BgBGrid.get(nbi,nbj,nbk); - +void calculateWaveSpeedXY(FsGrid, FS_STENCIL_WIDTH>& perBGrid, + FsGrid, FS_STENCIL_WIDTH>& momentsGrid, + FsGrid, FS_STENCIL_WIDTH>& dPerBGrid, + FsGrid, FS_STENCIL_WIDTH>& dMomentsGrid, + FsGrid, FS_STENCIL_WIDTH>& BgBGrid, + cint i, + cint j, + cint k, + cint nbi, + cint nbj, + cint nbk, + const Real& Bx, + const Real& By, + const Real& dBxdy, + const Real& dBxdz, + const Real& dBydx, + const Real& dBydz, + const Real& xdir, + const Real& ydir, + const Real& minRhom, + const Real& maxRhom, + Real& ret_vA, + Real& ret_vS, + Real& ret_vW) { + std::array* perb = perBGrid.get(i, j, k); + std::array* nbr_perb = perBGrid.get(nbi, nbj, nbk); + std::array* moments = momentsGrid.get(i, j, k); + std::array* dmoments = dMomentsGrid.get(i, j, k); + std::array* dperb = dPerBGrid.get(i, j, k); + std::array* nbr_dperb = dPerBGrid.get(nbi, nbj, nbk); + std::array* bgb = BgBGrid.get(i, j, k); + std::array* nbr_bgb = BgBGrid.get(nbi, nbj, nbk); + Real C_0, C_Z, rhom, p11, p22, p33; - C_0 = HALF*(nbr_perb->at(fsgrids::bfield::PERBZ) + nbr_bgb->at(fsgrids::bgbfield::BGBZ) + perb->at(fsgrids::bfield::PERBZ) + bgb->at(fsgrids::bgbfield::BGBZ)); - C_Z = (nbr_perb->at(fsgrids::bfield::PERBZ) + nbr_bgb->at(fsgrids::bgbfield::BGBZ)) - (perb->at(fsgrids::bfield::PERBZ) + bgb->at(fsgrids::bgbfield::BGBZ)); - rhom = moments->at(fsgrids::moments::RHOM) + xdir*HALF*dmoments->at(fsgrids::dmoments::drhomdx) + ydir*HALF*dmoments->at(fsgrids::dmoments::drhomdy); - p11 = moments->at(fsgrids::moments::P_11) + xdir*HALF*dmoments->at(fsgrids::dmoments::dp11dx) + ydir*HALF*dmoments->at(fsgrids::dmoments::dp11dy); - p22 = moments->at(fsgrids::moments::P_22) + xdir*HALF*dmoments->at(fsgrids::dmoments::dp22dx) + ydir*HALF*dmoments->at(fsgrids::dmoments::dp22dy); - p33 = moments->at(fsgrids::moments::P_33) + xdir*HALF*dmoments->at(fsgrids::dmoments::dp33dx) + ydir*HALF*dmoments->at(fsgrids::dmoments::dp33dy); - + C_0 = HALF * (nbr_perb->at(fsgrids::bfield::PERBZ) + nbr_bgb->at(fsgrids::bgbfield::BGBZ) + + perb->at(fsgrids::bfield::PERBZ) + bgb->at(fsgrids::bgbfield::BGBZ)); + C_Z = (nbr_perb->at(fsgrids::bfield::PERBZ) + nbr_bgb->at(fsgrids::bgbfield::BGBZ)) - + (perb->at(fsgrids::bfield::PERBZ) + bgb->at(fsgrids::bgbfield::BGBZ)); + rhom = moments->at(fsgrids::moments::RHOM) + xdir * HALF * dmoments->at(fsgrids::dmoments::drhomdx) + + ydir * HALF * dmoments->at(fsgrids::dmoments::drhomdy); + p11 = moments->at(fsgrids::moments::P_11) + xdir * HALF * dmoments->at(fsgrids::dmoments::dp11dx) + + ydir * HALF * dmoments->at(fsgrids::dmoments::dp11dy); + p22 = moments->at(fsgrids::moments::P_22) + xdir * HALF * dmoments->at(fsgrids::dmoments::dp22dx) + + ydir * HALF * dmoments->at(fsgrids::dmoments::dp22dy); + p33 = moments->at(fsgrids::moments::P_33) + xdir * HALF * dmoments->at(fsgrids::dmoments::dp33dx) + + ydir * HALF * dmoments->at(fsgrids::dmoments::dp33dy); + if (rhom < minRhom) { rhom = minRhom; } else if (rhom > maxRhom) { rhom = maxRhom; } - - const Real C_X = nbr_dperb->at(fsgrids::dperb::dPERBzdx) + nbr_bgb->at(fsgrids::bgbfield::dBGBzdx) + dperb->at(fsgrids::dperb::dPERBzdx) + bgb->at(fsgrids::bgbfield::dBGBzdx); - const Real C_XZ = nbr_dperb->at(fsgrids::dperb::dPERBzdx) + nbr_bgb->at(fsgrids::bgbfield::dBGBzdx) - (dperb->at(fsgrids::dperb::dPERBzdx) + bgb->at(fsgrids::bgbfield::dBGBzdx)); - const Real C_Y = nbr_dperb->at(fsgrids::dperb::dPERBzdy) + nbr_bgb->at(fsgrids::bgbfield::dBGBzdy) + dperb->at(fsgrids::dperb::dPERBzdy) + bgb->at(fsgrids::bgbfield::dBGBzdy); - const Real C_YZ = nbr_dperb->at(fsgrids::dperb::dPERBzdy) + nbr_bgb->at(fsgrids::bgbfield::dBGBzdy) - (dperb->at(fsgrids::dperb::dPERBzdy) + bgb->at(fsgrids::bgbfield::dBGBzdy)); - - const Real Bz2 = (C_0 + xdir*HALF*C_X + ydir*HALF*C_Y)*(C_0 + xdir*HALF*C_X + ydir*HALF*C_Y) - + TWELWTH*(C_Z + xdir*HALF*C_XZ + ydir*HALF*C_YZ)*(C_Z + xdir*HALF*C_XZ + ydir*HALF*C_YZ); - const Real Bx2 = (Bx + ydir*HALF*dBxdy)*(Bx + ydir*HALF*dBxdy) + TWELWTH*dBxdz*dBxdz; - const Real By2 = (By + xdir*HALF*dBydx)*(By + xdir*HALF*dBydx) + TWELWTH*dBydz*dBydz; - + + const Real C_X = nbr_dperb->at(fsgrids::dperb::dPERBzdx) + nbr_bgb->at(fsgrids::bgbfield::dBGBzdx) + + dperb->at(fsgrids::dperb::dPERBzdx) + bgb->at(fsgrids::bgbfield::dBGBzdx); + const Real C_XZ = nbr_dperb->at(fsgrids::dperb::dPERBzdx) + nbr_bgb->at(fsgrids::bgbfield::dBGBzdx) - + (dperb->at(fsgrids::dperb::dPERBzdx) + bgb->at(fsgrids::bgbfield::dBGBzdx)); + const Real C_Y = nbr_dperb->at(fsgrids::dperb::dPERBzdy) + nbr_bgb->at(fsgrids::bgbfield::dBGBzdy) + + dperb->at(fsgrids::dperb::dPERBzdy) + bgb->at(fsgrids::bgbfield::dBGBzdy); + const Real C_YZ = nbr_dperb->at(fsgrids::dperb::dPERBzdy) + nbr_bgb->at(fsgrids::bgbfield::dBGBzdy) - + (dperb->at(fsgrids::dperb::dPERBzdy) + bgb->at(fsgrids::bgbfield::dBGBzdy)); + + const Real Bz2 = + (C_0 + xdir * HALF * C_X + ydir * HALF * C_Y) * (C_0 + xdir * HALF * C_X + ydir * HALF * C_Y) + + TWELWTH * (C_Z + xdir * HALF * C_XZ + ydir * HALF * C_YZ) * (C_Z + xdir * HALF * C_XZ + ydir * HALF * C_YZ); + const Real Bx2 = (Bx + ydir * HALF * dBxdy) * (Bx + ydir * HALF * dBxdy) + TWELWTH * dBxdz * dBxdz; + const Real By2 = (By + xdir * HALF * dBydx) * (By + xdir * HALF * dBydx) + TWELWTH * dBydz * dBydz; + const Real Bmag2 = Bx2 + By2 + Bz2; - + p11 = p11 < 0.0 ? 0.0 : p11; p22 = p22 < 0.0 ? 0.0 : p22; p33 = p33 < 0.0 ? 0.0 : p33; - + // Effective wave speeds for advection and CFL calculation // Note that these are calculated as if the plasma is purely made up of hydrogen, which // is a reasonable approximation if it is proton-dominant. @@ -418,27 +453,36 @@ void calculateWaveSpeedXY( // and // http://iopscience.iop.org/article/10.1088/0253-6102/43/2/026/meta (Alfven waves) // for details. - const Real vA2 = divideIfNonZero(Bmag2, pc::MU_0*rhom); // Alfven speed - const Real vS2 = divideIfNonZero(p11+p22+p33, 2.0*rhom); // sound speed, adiabatic coefficient 3/2, P=1/3*trace in sound speed -// const Real vW = Parameters::ohmHallTerm > 0 ? divideIfNonZero(2.0*M_PI*vA2*pc::MASS_PROTON, perBGrid.DX*pc::CHARGE*sqrt(Bmag2)) : 0.0; // whistler speed - const Real vW = Parameters::ohmHallTerm > 0 ? - sqrt(vA2) * (1 + divideIfNonZero(2*M_PI*M_PI*pc::MASS_PROTON*pc::MASS_PROTON, perBGrid.DX*perBGrid.DX*rhom*pc::CHARGE*pc::CHARGE*pc::MU_0) - / sqrt(1 + divideIfNonZero( M_PI*M_PI*pc::MASS_PROTON*pc::MASS_PROTON, perBGrid.DX*perBGrid.DX*rhom*pc::CHARGE*pc::CHARGE*pc::MU_0))) - : 0.0; // whistler speed - + const Real vA2 = divideIfNonZero(Bmag2, pc::MU_0 * rhom); // Alfven speed + const Real vS2 = divideIfNonZero(p11 + p22 + p33, + 2.0 * rhom); // sound speed, adiabatic coefficient 3/2, P=1/3*trace in sound speed + // const Real vW = Parameters::ohmHallTerm > 0 ? divideIfNonZero(2.0*M_PI*vA2*pc::MASS_PROTON, + // perBGrid.DX*pc::CHARGE*sqrt(Bmag2)) : 0.0; // whistler speed + const Real vW = + Parameters::ohmHallTerm > 0 + ? sqrt(vA2) * (1 + divideIfNonZero(2 * M_PI * M_PI * pc::MASS_PROTON * pc::MASS_PROTON, + perBGrid.DX * perBGrid.DX * rhom * pc::CHARGE * pc::CHARGE * pc::MU_0) / + sqrt(1 + divideIfNonZero(M_PI * M_PI * pc::MASS_PROTON * pc::MASS_PROTON, + perBGrid.DX * perBGrid.DX * rhom * pc::CHARGE * pc::CHARGE * + pc::MU_0))) + : 0.0; // whistler speed + ret_vA = sqrt(vA2); ret_vS = sqrt(vS2); ret_vW = vW; } /*! \brief Low-level electric field propagation function. - * - * Computes the upwinded electric field X component along the cell's corresponding edge as the cross product of B and V in the YZ plane. Also includes the calculation of the maximally allowed time step. - * - * Expects that the correct RHO and B fields are being passed, depending on the stage of the Runge-Kutta time stepping method. - * - * Note that the background B field is excluded from the diffusive term calculations because they are equivalent to a current term and the background field is curl-free. - * + * + * Computes the upwinded electric field X component along the cell's corresponding edge as the cross product of B and V + * in the YZ plane. Also includes the calculation of the maximally allowed time step. + * + * Expects that the correct RHO and B fields are being passed, depending on the stage of the Runge-Kutta time stepping + * method. + * + * Note that the background B field is excluded from the diffusive term calculations because they are equivalent to a + * current term and the background field is curl-free. + * * \param perBGrid fsGrid holding the perturbed B quantities * \param EGrid fsGrid holding the electric field * \param EHallGrid fsGrid holding the Hall contributions to the electric field @@ -452,96 +496,91 @@ void calculateWaveSpeedXY( * \param RKCase Element in the enum defining the Runge-Kutta method steps */ void calculateEdgeElectricFieldX( - FsGrid< std::array, FS_STENCIL_WIDTH> & perBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EHallGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, - FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, - cint i, - cint j, - cint k, - cint& RKCase -) { - #ifdef DEBUG_FSOLVER + FsGrid, FS_STENCIL_WIDTH>& perBGrid, + FsGrid, FS_STENCIL_WIDTH>& EGrid, + FsGrid, FS_STENCIL_WIDTH>& EHallGrid, + FsGrid, FS_STENCIL_WIDTH>& EGradPeGrid, + FsGrid, FS_STENCIL_WIDTH>& momentsGrid, + FsGrid, FS_STENCIL_WIDTH>& dPerBGrid, + FsGrid, FS_STENCIL_WIDTH>& dMomentsGrid, + FsGrid, FS_STENCIL_WIDTH>& BgBGrid, + FsGrid& technicalGrid, + cint i, + cint j, + cint k, + cint& RKCase) { +#ifdef DEBUG_FSOLVER bool ok = true; - if (technicalGrid.get(i,j,k) == NULL) ok = false; - if (technicalGrid.get(i,j-1,k) == NULL) ok = false; - if (technicalGrid.get(i,j-1,k-1) == NULL) ok = false; - if (technicalGrid.get(i,j,k-1) == NULL) ok = false; + if (technicalGrid.get(i, j, k) == NULL) + ok = false; + if (technicalGrid.get(i, j - 1, k) == NULL) + ok = false; + if (technicalGrid.get(i, j - 1, k - 1) == NULL) + ok = false; + if (technicalGrid.get(i, j, k - 1) == NULL) + ok = false; if (ok == false) { cerr << "NULL pointer in " << __FILE__ << ":" << __LINE__ << std::endl; exit(1); } - #endif +#endif // An edge has four neighbouring spatial cells. Calculate // electric field in each of the four cells per edge. - Real ay_pos,ay_neg; // Max. characteristic velocities to y-direction - Real az_pos,az_neg; // Max. characteristic velocities to z-direction - Real Vy0,Vz0; // Reconstructed V - Real vA, vS, vW; // Alfven, sound, whistler speed - Real maxV = 0.0; // Max velocity for CFL purposes - Real c_y, c_z; // Wave speeds to yz-directions + Real ay_pos, ay_neg; // Max. characteristic velocities to y-direction + Real az_pos, az_neg; // Max. characteristic velocities to z-direction + Real Vy0, Vz0; // Reconstructed V + Real vA, vS, vW; // Alfven, sound, whistler speed + Real maxV = 0.0; // Max velocity for CFL purposes + Real c_y, c_z; // Wave speeds to yz-directions // Get values at all four neighbours, result is written to SW. - std::array * perb_SW = perBGrid.get(i ,j ,k ); - std::array * perb_SE = perBGrid.get(i ,j-1,k ); - std::array * perb_NE = perBGrid.get(i ,j-1,k-1); - std::array * perb_NW = perBGrid.get(i ,j ,k-1); - std::array * bgb_SW = BgBGrid.get(i,j ,k ); - std::array * bgb_SE = BgBGrid.get(i,j-1,k ); - std::array * bgb_NE = BgBGrid.get(i,j-1,k-1); - std::array * bgb_NW = BgBGrid.get(i,j ,k-1); - std::array * moments_SW = momentsGrid.get(i ,j ,k ); - std::array * moments_SE = momentsGrid.get(i ,j-1,k ); - std::array * moments_NE = momentsGrid.get(i ,j-1,k-1); - std::array * moments_NW = momentsGrid.get(i ,j ,k-1); - std::array * dmoments_SW = dMomentsGrid.get(i ,j ,k ); - std::array * dmoments_SE = dMomentsGrid.get(i ,j-1,k ); - std::array * dmoments_NE = dMomentsGrid.get(i ,j-1,k-1); - std::array * dmoments_NW = dMomentsGrid.get(i ,j ,k-1); - std::array * dperb_SW = dPerBGrid.get(i ,j ,k ); - std::array * dperb_SE = dPerBGrid.get(i ,j-1,k ); - std::array * dperb_NE = dPerBGrid.get(i ,j-1,k-1); - std::array * dperb_NW = dPerBGrid.get(i ,j ,k-1); - - std::array * efield_SW = EGrid.get(i,j,k); - + std::array* perb_SW = perBGrid.get(i, j, k); + std::array* perb_SE = perBGrid.get(i, j - 1, k); + std::array* perb_NE = perBGrid.get(i, j - 1, k - 1); + std::array* perb_NW = perBGrid.get(i, j, k - 1); + std::array* bgb_SW = BgBGrid.get(i, j, k); + std::array* bgb_SE = BgBGrid.get(i, j - 1, k); + std::array* bgb_NE = BgBGrid.get(i, j - 1, k - 1); + std::array* bgb_NW = BgBGrid.get(i, j, k - 1); + std::array* moments_SW = momentsGrid.get(i, j, k); + std::array* moments_SE = momentsGrid.get(i, j - 1, k); + std::array* moments_NE = momentsGrid.get(i, j - 1, k - 1); + std::array* moments_NW = momentsGrid.get(i, j, k - 1); + std::array* dmoments_SW = dMomentsGrid.get(i, j, k); + std::array* dmoments_SE = dMomentsGrid.get(i, j - 1, k); + std::array* dmoments_NE = dMomentsGrid.get(i, j - 1, k - 1); + std::array* dmoments_NW = dMomentsGrid.get(i, j, k - 1); + std::array* dperb_SW = dPerBGrid.get(i, j, k); + std::array* dperb_SE = dPerBGrid.get(i, j - 1, k); + std::array* dperb_NE = dPerBGrid.get(i, j - 1, k - 1); + std::array* dperb_NW = dPerBGrid.get(i, j, k - 1); + + std::array* efield_SW = EGrid.get(i, j, k); + Real By_S, Bz_W, Bz_E, By_N, perBy_S, perBz_W, perBz_E, perBy_N; Real minRhom = std::numeric_limits::max(); Real maxRhom = std::numeric_limits::min(); - By_S = perb_SW->at(fsgrids::bfield::PERBY)+bgb_SW->at(fsgrids::bgbfield::BGBY); - Bz_W = perb_SW->at(fsgrids::bfield::PERBZ)+bgb_SW->at(fsgrids::bgbfield::BGBZ); - Bz_E = perb_SE->at(fsgrids::bfield::PERBZ)+bgb_SE->at(fsgrids::bgbfield::BGBZ); - By_N = perb_NW->at(fsgrids::bfield::PERBY)+bgb_NW->at(fsgrids::bgbfield::BGBY); + By_S = perb_SW->at(fsgrids::bfield::PERBY) + bgb_SW->at(fsgrids::bgbfield::BGBY); + Bz_W = perb_SW->at(fsgrids::bfield::PERBZ) + bgb_SW->at(fsgrids::bgbfield::BGBZ); + Bz_E = perb_SE->at(fsgrids::bfield::PERBZ) + bgb_SE->at(fsgrids::bgbfield::BGBZ); + By_N = perb_NW->at(fsgrids::bfield::PERBY) + bgb_NW->at(fsgrids::bgbfield::BGBY); perBy_S = perb_SW->at(fsgrids::bfield::PERBY); perBz_W = perb_SW->at(fsgrids::bfield::PERBZ); perBz_E = perb_SE->at(fsgrids::bfield::PERBZ); perBy_N = perb_NW->at(fsgrids::bfield::PERBY); - Vy0 = moments_SW->at(fsgrids::moments::VY); - Vz0 = moments_SW->at(fsgrids::moments::VZ); + Vy0 = moments_SW->at(fsgrids::moments::VY); + Vz0 = moments_SW->at(fsgrids::moments::VZ); minRhom = min(minRhom, - min(moments_SW->at(fsgrids::moments::RHOM), - min(moments_SE->at(fsgrids::moments::RHOM), - min(moments_NW->at(fsgrids::moments::RHOM), - moments_NE->at(fsgrids::moments::RHOM)) - ) - ) - ); + min(moments_SW->at(fsgrids::moments::RHOM), + min(moments_SE->at(fsgrids::moments::RHOM), + min(moments_NW->at(fsgrids::moments::RHOM), moments_NE->at(fsgrids::moments::RHOM))))); maxRhom = max(maxRhom, - max(moments_SW->at(fsgrids::moments::RHOM), - max(moments_SE->at(fsgrids::moments::RHOM), - max(moments_NW->at(fsgrids::moments::RHOM), - moments_NE->at(fsgrids::moments::RHOM)) - ) - ) - ); - + max(moments_SW->at(fsgrids::moments::RHOM), + max(moments_SE->at(fsgrids::moments::RHOM), + max(moments_NW->at(fsgrids::moments::RHOM), moments_NE->at(fsgrids::moments::RHOM))))); + creal dBydx_S = dperb_SW->at(fsgrids::dperb::dPERBydx) + bgb_SW->at(fsgrids::bgbfield::dBGBydx); creal dBydz_S = dperb_SW->at(fsgrids::dperb::dPERBydz) + bgb_SW->at(fsgrids::bgbfield::dBGBydz); creal dBzdx_W = dperb_SW->at(fsgrids::dperb::dPERBzdx) + bgb_SW->at(fsgrids::bgbfield::dBGBzdx); @@ -557,348 +596,418 @@ void calculateEdgeElectricFieldX( // Ex and characteristic speeds on this cell: // 1st order terms: - Real Ex_SW = By_S*Vz0 - Bz_W*Vy0; + Real Ex_SW = By_S * Vz0 - Bz_W * Vy0; // Resistive term if (Parameters::resistivity > 0) { - Ex_SW += Parameters::resistivity * - sqrt((bgb_SW->at(fsgrids::bgbfield::BGBX)+perb_SW->at(fsgrids::bfield::PERBX))* - (bgb_SW->at(fsgrids::bgbfield::BGBX)+perb_SW->at(fsgrids::bfield::PERBX)) + - (bgb_SW->at(fsgrids::bgbfield::BGBY)+perb_SW->at(fsgrids::bfield::PERBY))* - (bgb_SW->at(fsgrids::bgbfield::BGBY)+perb_SW->at(fsgrids::bfield::PERBY)) + - (bgb_SW->at(fsgrids::bgbfield::BGBZ)+perb_SW->at(fsgrids::bfield::PERBZ))* - (bgb_SW->at(fsgrids::bgbfield::BGBZ)+perb_SW->at(fsgrids::bfield::PERBZ)) - ) / - moments_SW->at(fsgrids::moments::RHOQ) / - physicalconstants::MU_0 * - (dperb_SW->at(fsgrids::dperb::dPERBzdy)/technicalGrid.DY - dperb_SW->at(fsgrids::dperb::dPERBydz)/technicalGrid.DZ); + Ex_SW += Parameters::resistivity * + sqrt((bgb_SW->at(fsgrids::bgbfield::BGBX) + perb_SW->at(fsgrids::bfield::PERBX)) * + (bgb_SW->at(fsgrids::bgbfield::BGBX) + perb_SW->at(fsgrids::bfield::PERBX)) + + (bgb_SW->at(fsgrids::bgbfield::BGBY) + perb_SW->at(fsgrids::bfield::PERBY)) * + (bgb_SW->at(fsgrids::bgbfield::BGBY) + perb_SW->at(fsgrids::bfield::PERBY)) + + (bgb_SW->at(fsgrids::bgbfield::BGBZ) + perb_SW->at(fsgrids::bfield::PERBZ)) * + (bgb_SW->at(fsgrids::bgbfield::BGBZ) + perb_SW->at(fsgrids::bfield::PERBZ))) / + moments_SW->at(fsgrids::moments::RHOQ) / physicalconstants::MU_0 * + (dperb_SW->at(fsgrids::dperb::dPERBzdy) / technicalGrid.DY - + dperb_SW->at(fsgrids::dperb::dPERBydz) / technicalGrid.DZ); } - + // Hall term - if(Parameters::ohmHallTerm > 0) { - Ex_SW += EHallGrid.get(i,j,k)->at(fsgrids::ehall::EXHALL_000_100); + if (Parameters::ohmHallTerm > 0) { + Ex_SW += EHallGrid.get(i, j, k)->at(fsgrids::ehall::EXHALL_000_100); } - + // Electron pressure gradient term - if(Parameters::ohmGradPeTerm > 0) { - Ex_SW += EGradPeGrid.get(i,j,k)->at(fsgrids::egradpe::EXGRADPE); + if (Parameters::ohmGradPeTerm > 0) { + Ex_SW += EGradPeGrid.get(i, j, k)->at(fsgrids::egradpe::EXGRADPE); } - #ifndef FS_1ST_ORDER_SPACE - // 2nd order terms: - Ex_SW += +HALF*((By_S - HALF*dBydz_S)*(-dmoments_SW->at(fsgrids::dmoments::dVzdy) - dmoments_SW->at(fsgrids::dmoments::dVzdz)) - dBydz_S*Vz0 + SIXTH*dBydx_S*dmoments_SW->at(fsgrids::dmoments::dVzdx)); - Ex_SW += -HALF*((Bz_W - HALF*dBzdy_W)*(-dmoments_SW->at(fsgrids::dmoments::dVydy) - dmoments_SW->at(fsgrids::dmoments::dVydz)) - dBzdy_W*Vy0 + SIXTH*dBzdx_W*dmoments_SW->at(fsgrids::dmoments::dVydx)); - #endif - calculateWaveSpeedYZ( - perBGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - i , j, k, - i+1, j, k, - By_S, Bz_W, dBydx_S, dBydz_S, dBzdx_W, dBzdy_W, MINUS, MINUS, minRhom, maxRhom, vA, vS, vW - ); - c_y = min(Parameters::maxWaveVelocity,sqrt(vA*vA + vS*vS) + vW); +#ifndef FS_1ST_ORDER_SPACE + // 2nd order terms: + Ex_SW += +HALF * ((By_S - HALF * dBydz_S) * + (-dmoments_SW->at(fsgrids::dmoments::dVzdy) - dmoments_SW->at(fsgrids::dmoments::dVzdz)) - + dBydz_S * Vz0 + SIXTH * dBydx_S * dmoments_SW->at(fsgrids::dmoments::dVzdx)); + Ex_SW += -HALF * ((Bz_W - HALF * dBzdy_W) * + (-dmoments_SW->at(fsgrids::dmoments::dVydy) - dmoments_SW->at(fsgrids::dmoments::dVydz)) - + dBzdy_W * Vy0 + SIXTH * dBzdx_W * dmoments_SW->at(fsgrids::dmoments::dVydx)); +#endif + calculateWaveSpeedYZ(perBGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + i, + j, + k, + i + 1, + j, + k, + By_S, + Bz_W, + dBydx_S, + dBydz_S, + dBzdx_W, + dBzdy_W, + MINUS, + MINUS, + minRhom, + maxRhom, + vA, + vS, + vW); + c_y = min(Parameters::maxWaveVelocity, sqrt(vA * vA + vS * vS) + vW); c_z = c_y; - ay_neg = max(ZERO,-Vy0 + c_y); - ay_pos = max(ZERO,+Vy0 + c_y); - az_neg = max(ZERO,-Vz0 + c_z); - az_pos = max(ZERO,+Vz0 + c_z); + ay_neg = max(ZERO, -Vy0 + c_y); + ay_pos = max(ZERO, +Vy0 + c_y); + az_neg = max(ZERO, -Vz0 + c_z); + az_pos = max(ZERO, +Vz0 + c_z); maxV = max(maxV, calculateCflSpeed(Vy0, Vz0, vA, vS, vW)); // Ex and characteristic speeds on j-1 neighbour: - Vy0 = moments_SE->at(fsgrids::moments::VY); - Vz0 = moments_SE->at(fsgrids::moments::VZ); - + Vy0 = moments_SE->at(fsgrids::moments::VY); + Vz0 = moments_SE->at(fsgrids::moments::VZ); + // 1st order terms: - Real Ex_SE = By_S*Vz0 - Bz_E*Vy0; - + Real Ex_SE = By_S * Vz0 - Bz_E * Vy0; + // Resistive term if (Parameters::resistivity > 0) { - Ex_SE += Parameters::resistivity * - sqrt((bgb_SE->at(fsgrids::bgbfield::BGBX)+perb_SE->at(fsgrids::bfield::PERBX))* - (bgb_SE->at(fsgrids::bgbfield::BGBX)+perb_SE->at(fsgrids::bfield::PERBX)) + - (bgb_SE->at(fsgrids::bgbfield::BGBY)+perb_SE->at(fsgrids::bfield::PERBY))* - (bgb_SE->at(fsgrids::bgbfield::BGBY)+perb_SE->at(fsgrids::bfield::PERBY)) + - (bgb_SE->at(fsgrids::bgbfield::BGBZ)+perb_SE->at(fsgrids::bfield::PERBZ))* - (bgb_SE->at(fsgrids::bgbfield::BGBZ)+perb_SE->at(fsgrids::bfield::PERBZ)) - ) / - moments_SE->at(fsgrids::moments::RHOQ) / - physicalconstants::MU_0 * - (dperb_SE->at(fsgrids::dperb::dPERBzdy)/technicalGrid.DY - dperb_SE->at(fsgrids::dperb::dPERBydz)/technicalGrid.DZ); + Ex_SE += Parameters::resistivity * + sqrt((bgb_SE->at(fsgrids::bgbfield::BGBX) + perb_SE->at(fsgrids::bfield::PERBX)) * + (bgb_SE->at(fsgrids::bgbfield::BGBX) + perb_SE->at(fsgrids::bfield::PERBX)) + + (bgb_SE->at(fsgrids::bgbfield::BGBY) + perb_SE->at(fsgrids::bfield::PERBY)) * + (bgb_SE->at(fsgrids::bgbfield::BGBY) + perb_SE->at(fsgrids::bfield::PERBY)) + + (bgb_SE->at(fsgrids::bgbfield::BGBZ) + perb_SE->at(fsgrids::bfield::PERBZ)) * + (bgb_SE->at(fsgrids::bgbfield::BGBZ) + perb_SE->at(fsgrids::bfield::PERBZ))) / + moments_SE->at(fsgrids::moments::RHOQ) / physicalconstants::MU_0 * + (dperb_SE->at(fsgrids::dperb::dPERBzdy) / technicalGrid.DY - + dperb_SE->at(fsgrids::dperb::dPERBydz) / technicalGrid.DZ); } // Hall term - if(Parameters::ohmHallTerm > 0) { - Ex_SE += EHallGrid.get(i,j-1,k)->at(fsgrids::ehall::EXHALL_010_110); + if (Parameters::ohmHallTerm > 0) { + Ex_SE += EHallGrid.get(i, j - 1, k)->at(fsgrids::ehall::EXHALL_010_110); } - + // Electron pressure gradient term - if(Parameters::ohmGradPeTerm > 0) { - Ex_SE += EGradPeGrid.get(i,j-1,k)->at(fsgrids::egradpe::EXGRADPE); + if (Parameters::ohmGradPeTerm > 0) { + Ex_SE += EGradPeGrid.get(i, j - 1, k)->at(fsgrids::egradpe::EXGRADPE); } - - #ifndef FS_1ST_ORDER_SPACE - // 2nd order terms: - Ex_SE += +HALF*((By_S - HALF*dBydz_S)*(+dmoments_SE->at(fsgrids::dmoments::dVzdy) - dmoments_SE->at(fsgrids::dmoments::dVzdz)) - dBydz_S*Vz0 + SIXTH*dBydx_S*dmoments_SE->at(fsgrids::dmoments::dVzdx)); - Ex_SE += -HALF*((Bz_E + HALF*dBzdy_E)*(+dmoments_SE->at(fsgrids::dmoments::dVydy) - dmoments_SE->at(fsgrids::dmoments::dVydz)) + dBzdy_E*Vy0 + SIXTH*dBzdx_E*dmoments_SE->at(fsgrids::dmoments::dVydx)); - #endif - - calculateWaveSpeedYZ( - perBGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - i , j-1, k, - i+1, j-1, k, - By_S, Bz_E, dBydx_S, dBydz_S, dBzdx_E, dBzdy_E, PLUS, MINUS, minRhom, maxRhom, vA, vS, vW - ); - c_y = min(Parameters::maxWaveVelocity,sqrt(vA*vA + vS*vS) + vW); + +#ifndef FS_1ST_ORDER_SPACE + // 2nd order terms: + Ex_SE += +HALF * ((By_S - HALF * dBydz_S) * + (+dmoments_SE->at(fsgrids::dmoments::dVzdy) - dmoments_SE->at(fsgrids::dmoments::dVzdz)) - + dBydz_S * Vz0 + SIXTH * dBydx_S * dmoments_SE->at(fsgrids::dmoments::dVzdx)); + Ex_SE += -HALF * ((Bz_E + HALF * dBzdy_E) * + (+dmoments_SE->at(fsgrids::dmoments::dVydy) - dmoments_SE->at(fsgrids::dmoments::dVydz)) + + dBzdy_E * Vy0 + SIXTH * dBzdx_E * dmoments_SE->at(fsgrids::dmoments::dVydx)); +#endif + + calculateWaveSpeedYZ(perBGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + i, + j - 1, + k, + i + 1, + j - 1, + k, + By_S, + Bz_E, + dBydx_S, + dBydz_S, + dBzdx_E, + dBzdy_E, + PLUS, + MINUS, + minRhom, + maxRhom, + vA, + vS, + vW); + c_y = min(Parameters::maxWaveVelocity, sqrt(vA * vA + vS * vS) + vW); c_z = c_y; - ay_neg = max(ay_neg,-Vy0 + c_y); - ay_pos = max(ay_pos,+Vy0 + c_y); - az_neg = max(az_neg,-Vz0 + c_z); - az_pos = max(az_pos,+Vz0 + c_z); + ay_neg = max(ay_neg, -Vy0 + c_y); + ay_pos = max(ay_pos, +Vy0 + c_y); + az_neg = max(az_neg, -Vz0 + c_z); + az_pos = max(az_pos, +Vz0 + c_z); maxV = max(maxV, calculateCflSpeed(Vy0, Vz0, vA, vS, vW)); // Ex and characteristic speeds on k-1 neighbour: - Vy0 = moments_NW->at(fsgrids::moments::VY); - Vz0 = moments_NW->at(fsgrids::moments::VZ); - + Vy0 = moments_NW->at(fsgrids::moments::VY); + Vz0 = moments_NW->at(fsgrids::moments::VZ); + // 1st order terms: - Real Ex_NW = By_N*Vz0 - Bz_W*Vy0; - + Real Ex_NW = By_N * Vz0 - Bz_W * Vy0; + // Resistive term if (Parameters::resistivity > 0) { - Ex_NW += Parameters::resistivity * - sqrt((bgb_NW->at(fsgrids::bgbfield::BGBX)+perb_NW->at(fsgrids::bfield::PERBX))* - (bgb_NW->at(fsgrids::bgbfield::BGBX)+perb_NW->at(fsgrids::bfield::PERBX)) + - (bgb_NW->at(fsgrids::bgbfield::BGBY)+perb_NW->at(fsgrids::bfield::PERBY))* - (bgb_NW->at(fsgrids::bgbfield::BGBY)+perb_NW->at(fsgrids::bfield::PERBY)) + - (bgb_NW->at(fsgrids::bgbfield::BGBZ)+perb_NW->at(fsgrids::bfield::PERBZ))* - (bgb_NW->at(fsgrids::bgbfield::BGBZ)+perb_NW->at(fsgrids::bfield::PERBZ)) - ) / - moments_NW->at(fsgrids::moments::RHOQ) / - physicalconstants::MU_0 * - (dperb_NW->at(fsgrids::dperb::dPERBzdy)/technicalGrid.DY - dperb_NW->at(fsgrids::dperb::dPERBydz)/technicalGrid.DZ); + Ex_NW += Parameters::resistivity * + sqrt((bgb_NW->at(fsgrids::bgbfield::BGBX) + perb_NW->at(fsgrids::bfield::PERBX)) * + (bgb_NW->at(fsgrids::bgbfield::BGBX) + perb_NW->at(fsgrids::bfield::PERBX)) + + (bgb_NW->at(fsgrids::bgbfield::BGBY) + perb_NW->at(fsgrids::bfield::PERBY)) * + (bgb_NW->at(fsgrids::bgbfield::BGBY) + perb_NW->at(fsgrids::bfield::PERBY)) + + (bgb_NW->at(fsgrids::bgbfield::BGBZ) + perb_NW->at(fsgrids::bfield::PERBZ)) * + (bgb_NW->at(fsgrids::bgbfield::BGBZ) + perb_NW->at(fsgrids::bfield::PERBZ))) / + moments_NW->at(fsgrids::moments::RHOQ) / physicalconstants::MU_0 * + (dperb_NW->at(fsgrids::dperb::dPERBzdy) / technicalGrid.DY - + dperb_NW->at(fsgrids::dperb::dPERBydz) / technicalGrid.DZ); } - + // Hall term - if(Parameters::ohmHallTerm > 0) { - Ex_NW += EHallGrid.get(i,j,k-1)->at(fsgrids::ehall::EXHALL_001_101); + if (Parameters::ohmHallTerm > 0) { + Ex_NW += EHallGrid.get(i, j, k - 1)->at(fsgrids::ehall::EXHALL_001_101); } - + // Electron pressure gradient term - if(Parameters::ohmGradPeTerm > 0) { - Ex_NW += EGradPeGrid.get(i,j,k-1)->at(fsgrids::egradpe::EXGRADPE); + if (Parameters::ohmGradPeTerm > 0) { + Ex_NW += EGradPeGrid.get(i, j, k - 1)->at(fsgrids::egradpe::EXGRADPE); } - - #ifndef FS_1ST_ORDER_SPACE - // 2nd order terms: - Ex_NW += +HALF*((By_N + HALF*dBydz_N)*(-dmoments_NW->at(fsgrids::dmoments::dVzdy) + dmoments_NW->at(fsgrids::dmoments::dVzdz)) + dBydz_N*Vz0 + SIXTH*dBydx_N*dmoments_NW->at(fsgrids::dmoments::dVzdx)); - Ex_NW += -HALF*((Bz_W - HALF*dBzdy_W)*(-dmoments_NW->at(fsgrids::dmoments::dVydy) + dmoments_NW->at(fsgrids::dmoments::dVydz)) - dBzdy_W*Vy0 + SIXTH*dBzdx_W*dmoments_NW->at(fsgrids::dmoments::dVydx)); - #endif - - calculateWaveSpeedYZ( - perBGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - i , j, k-1, - i+1, j, k-1, - By_N, Bz_W, dBydx_N, dBydz_N, dBzdx_W, dBzdy_W, MINUS, PLUS, minRhom, maxRhom, vA, vS, vW - ); - c_y = min(Parameters::maxWaveVelocity,sqrt(vA*vA + vS*vS) + vW); + +#ifndef FS_1ST_ORDER_SPACE + // 2nd order terms: + Ex_NW += +HALF * ((By_N + HALF * dBydz_N) * + (-dmoments_NW->at(fsgrids::dmoments::dVzdy) + dmoments_NW->at(fsgrids::dmoments::dVzdz)) + + dBydz_N * Vz0 + SIXTH * dBydx_N * dmoments_NW->at(fsgrids::dmoments::dVzdx)); + Ex_NW += -HALF * ((Bz_W - HALF * dBzdy_W) * + (-dmoments_NW->at(fsgrids::dmoments::dVydy) + dmoments_NW->at(fsgrids::dmoments::dVydz)) - + dBzdy_W * Vy0 + SIXTH * dBzdx_W * dmoments_NW->at(fsgrids::dmoments::dVydx)); +#endif + + calculateWaveSpeedYZ(perBGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + i, + j, + k - 1, + i + 1, + j, + k - 1, + By_N, + Bz_W, + dBydx_N, + dBydz_N, + dBzdx_W, + dBzdy_W, + MINUS, + PLUS, + minRhom, + maxRhom, + vA, + vS, + vW); + c_y = min(Parameters::maxWaveVelocity, sqrt(vA * vA + vS * vS) + vW); c_z = c_y; - ay_neg = max(ay_neg,-Vy0 + c_y); - ay_pos = max(ay_pos,+Vy0 + c_y); - az_neg = max(az_neg,-Vz0 + c_z); - az_pos = max(az_pos,+Vz0 + c_z); + ay_neg = max(ay_neg, -Vy0 + c_y); + ay_pos = max(ay_pos, +Vy0 + c_y); + az_neg = max(az_neg, -Vz0 + c_z); + az_pos = max(az_pos, +Vz0 + c_z); maxV = max(maxV, calculateCflSpeed(Vy0, Vz0, vA, vS, vW)); // Ex and characteristic speeds on j-1,k-1 neighbour: Vy0 = moments_NE->at(fsgrids::moments::VY); Vz0 = moments_NE->at(fsgrids::moments::VZ); - + // 1st order terms: - Real Ex_NE = By_N*Vz0 - Bz_E*Vy0; + Real Ex_NE = By_N * Vz0 - Bz_E * Vy0; // Resistive term if (Parameters::resistivity > 0) { Ex_NE += Parameters::resistivity * - sqrt((bgb_NE->at(fsgrids::bgbfield::BGBX)+perb_NE->at(fsgrids::bfield::PERBX))* - (bgb_NE->at(fsgrids::bgbfield::BGBX)+perb_NE->at(fsgrids::bfield::PERBX)) + - (bgb_NE->at(fsgrids::bgbfield::BGBY)+perb_NE->at(fsgrids::bfield::PERBY))* - (bgb_NE->at(fsgrids::bgbfield::BGBY)+perb_NE->at(fsgrids::bfield::PERBY)) + - (bgb_NE->at(fsgrids::bgbfield::BGBZ)+perb_NE->at(fsgrids::bfield::PERBZ))* - (bgb_NE->at(fsgrids::bgbfield::BGBZ)+perb_NE->at(fsgrids::bfield::PERBZ)) - ) / - moments_NE->at(fsgrids::moments::RHOQ) / - physicalconstants::MU_0 * - (dperb_NE->at(fsgrids::dperb::dPERBzdy)/technicalGrid.DY - dperb_NE->at(fsgrids::dperb::dPERBydz)/technicalGrid.DZ); + sqrt((bgb_NE->at(fsgrids::bgbfield::BGBX) + perb_NE->at(fsgrids::bfield::PERBX)) * + (bgb_NE->at(fsgrids::bgbfield::BGBX) + perb_NE->at(fsgrids::bfield::PERBX)) + + (bgb_NE->at(fsgrids::bgbfield::BGBY) + perb_NE->at(fsgrids::bfield::PERBY)) * + (bgb_NE->at(fsgrids::bgbfield::BGBY) + perb_NE->at(fsgrids::bfield::PERBY)) + + (bgb_NE->at(fsgrids::bgbfield::BGBZ) + perb_NE->at(fsgrids::bfield::PERBZ)) * + (bgb_NE->at(fsgrids::bgbfield::BGBZ) + perb_NE->at(fsgrids::bfield::PERBZ))) / + moments_NE->at(fsgrids::moments::RHOQ) / physicalconstants::MU_0 * + (dperb_NE->at(fsgrids::dperb::dPERBzdy) / technicalGrid.DY - + dperb_NE->at(fsgrids::dperb::dPERBydz) / technicalGrid.DZ); } // Hall term - if(Parameters::ohmHallTerm > 0) { - Ex_NE += EHallGrid.get(i,j-1,k-1)->at(fsgrids::ehall::EXHALL_011_111); + if (Parameters::ohmHallTerm > 0) { + Ex_NE += EHallGrid.get(i, j - 1, k - 1)->at(fsgrids::ehall::EXHALL_011_111); } - + // Electron pressure gradient term - if(Parameters::ohmGradPeTerm > 0) { - Ex_NE += EGradPeGrid.get(i,j-1,k-1)->at(fsgrids::egradpe::EXGRADPE); + if (Parameters::ohmGradPeTerm > 0) { + Ex_NE += EGradPeGrid.get(i, j - 1, k - 1)->at(fsgrids::egradpe::EXGRADPE); } - - #ifndef FS_1ST_ORDER_SPACE - // 2nd order terms: - Ex_NE += +HALF*((By_N + HALF*dBydz_N)*(+dmoments_NE->at(fsgrids::dmoments::dVzdy) + dmoments_NE->at(fsgrids::dmoments::dVzdz)) + dBydz_N*Vz0 + SIXTH*dBydx_N*dmoments_NE->at(fsgrids::dmoments::dVzdx)); - Ex_NE += -HALF*((Bz_E + HALF*dBzdy_E)*(+dmoments_NE->at(fsgrids::dmoments::dVydy) + dmoments_NE->at(fsgrids::dmoments::dVydz)) + dBzdy_E*Vy0 + SIXTH*dBzdx_E*dmoments_NE->at(fsgrids::dmoments::dVydx)); - #endif - - calculateWaveSpeedYZ( - perBGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - i ,j-1,k-1, - i+1,j-1,k-1, - By_N, Bz_E, dBydx_N, dBydz_N, dBzdx_E, dBzdy_E, PLUS, PLUS, minRhom, maxRhom, vA, vS, vW - ); - c_y = min(Parameters::maxWaveVelocity,sqrt(vA*vA + vS*vS) + vW); + +#ifndef FS_1ST_ORDER_SPACE + // 2nd order terms: + Ex_NE += +HALF * ((By_N + HALF * dBydz_N) * + (+dmoments_NE->at(fsgrids::dmoments::dVzdy) + dmoments_NE->at(fsgrids::dmoments::dVzdz)) + + dBydz_N * Vz0 + SIXTH * dBydx_N * dmoments_NE->at(fsgrids::dmoments::dVzdx)); + Ex_NE += -HALF * ((Bz_E + HALF * dBzdy_E) * + (+dmoments_NE->at(fsgrids::dmoments::dVydy) + dmoments_NE->at(fsgrids::dmoments::dVydz)) + + dBzdy_E * Vy0 + SIXTH * dBzdx_E * dmoments_NE->at(fsgrids::dmoments::dVydx)); +#endif + + calculateWaveSpeedYZ(perBGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + i, + j - 1, + k - 1, + i + 1, + j - 1, + k - 1, + By_N, + Bz_E, + dBydx_N, + dBydz_N, + dBzdx_E, + dBzdy_E, + PLUS, + PLUS, + minRhom, + maxRhom, + vA, + vS, + vW); + c_y = min(Parameters::maxWaveVelocity, sqrt(vA * vA + vS * vS) + vW); c_z = c_y; - ay_neg = max(ay_neg,-Vy0 + c_y); - ay_pos = max(ay_pos,+Vy0 + c_y); - az_neg = max(az_neg,-Vz0 + c_z); - az_pos = max(az_pos,+Vz0 + c_z); + ay_neg = max(ay_neg, -Vy0 + c_y); + ay_pos = max(ay_pos, +Vy0 + c_y); + az_neg = max(az_neg, -Vz0 + c_z); + az_pos = max(az_pos, +Vz0 + c_z); maxV = max(maxV, calculateCflSpeed(Vy0, Vz0, vA, vS, vW)); // Calculate properly upwinded edge-averaged Ex: - efield_SW->at(fsgrids::efield::EX) = ay_pos*az_pos*Ex_NE + ay_pos*az_neg*Ex_SE + ay_neg*az_pos*Ex_NW + ay_neg*az_neg*Ex_SW; - efield_SW->at(fsgrids::efield::EX) /= ((ay_pos+ay_neg)*(az_pos+az_neg)+EPS); + efield_SW->at(fsgrids::efield::EX) = + ay_pos * az_pos * Ex_NE + ay_pos * az_neg * Ex_SE + ay_neg * az_pos * Ex_NW + ay_neg * az_neg * Ex_SW; + efield_SW->at(fsgrids::efield::EX) /= ((ay_pos + ay_neg) * (az_pos + az_neg) + EPS); if (Parameters::fieldSolverDiffusiveEterms) { #ifdef FS_1ST_ORDER_SPACE // 1st order diffusive terms: - efield_SW->at(fsgrids::efield::EX) -= az_pos*az_neg/(az_pos+az_neg+EPS)*(perBy_S-perBy_N); - efield_SW->at(fsgrids::efield::EX) += ay_pos*ay_neg/(ay_pos+ay_neg+EPS)*(perBz_W-perBz_E); + efield_SW->at(fsgrids::efield::EX) -= az_pos * az_neg / (az_pos + az_neg + EPS) * (perBy_S - perBy_N); + efield_SW->at(fsgrids::efield::EX) += ay_pos * ay_neg / (ay_pos + ay_neg + EPS) * (perBz_W - perBz_E); #else // 2nd order diffusive terms - efield_SW->at(fsgrids::efield::EX) -= az_pos*az_neg/(az_pos+az_neg+EPS)*((perBy_S-HALF*dperBydz_S) - (perBy_N+HALF*dperBydz_N)); - efield_SW->at(fsgrids::efield::EX) += ay_pos*ay_neg/(ay_pos+ay_neg+EPS)*((perBz_W-HALF*dperBzdy_W) - (perBz_E+HALF*dperBzdy_E)); + efield_SW->at(fsgrids::efield::EX) -= + az_pos * az_neg / (az_pos + az_neg + EPS) * ((perBy_S - HALF * dperBydz_S) - (perBy_N + HALF * dperBydz_N)); + efield_SW->at(fsgrids::efield::EX) += + ay_pos * ay_neg / (ay_pos + ay_neg + EPS) * ((perBz_W - HALF * dperBzdy_W) - (perBz_E + HALF * dperBzdy_E)); #endif } - + if ((RKCase == RK_ORDER1) || (RKCase == RK_ORDER2_STEP2)) { - //compute maximum timestep for fieldsolver in this cell (CFL=1) - Real min_dx=std::numeric_limits::max(); - min_dx=min(min_dx,technicalGrid.DY); - min_dx=min(min_dx,technicalGrid.DZ); - //update max allowed timestep for field propagation in this cell, which is the minimum of CFL=1 timesteps - if (maxV != ZERO) technicalGrid.get(i,j,k)->maxFsDt = min(technicalGrid.get(i,j,k)->maxFsDt,min_dx/maxV); + // compute maximum timestep for fieldsolver in this cell (CFL=1) + Real min_dx = std::numeric_limits::max(); + min_dx = min(min_dx, technicalGrid.DY); + min_dx = min(min_dx, technicalGrid.DZ); + // update max allowed timestep for field propagation in this cell, which is the minimum of CFL=1 timesteps + if (maxV != ZERO) + technicalGrid.get(i, j, k)->maxFsDt = min(technicalGrid.get(i, j, k)->maxFsDt, min_dx / maxV); } } /*! \brief Low-level electric field propagation function. - * - * Computes the upwinded electric field Y component along the cell's corresponding edge as the cross product of B and V in the XZ plane. Also includes the calculation of the maximally allowed time step. - * - * Expects that the correct RHO and B fields are being passed, depending on the stage of the Runge-Kutta time stepping method. - * - * Note that the background B field is excluded from the diffusive term calculations because they are equivalent to a current term and the background field is curl-free. - * + * + * Computes the upwinded electric field Y component along the cell's corresponding edge as the cross product of B and V + * in the XZ plane. Also includes the calculation of the maximally allowed time step. + * + * Expects that the correct RHO and B fields are being passed, depending on the stage of the Runge-Kutta time stepping + * method. + * + * Note that the background B field is excluded from the diffusive term calculations because they are equivalent to a + * current term and the background field is curl-free. + * * \param RKCase Element in the enum defining the Runge-Kutta method steps */ void calculateEdgeElectricFieldY( - FsGrid< std::array, FS_STENCIL_WIDTH> & perBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EHallGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, - FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, - cint i, - cint j, - cint k, - cint& RKCase -) { - #ifdef DEBUG_FSOLVER + FsGrid, FS_STENCIL_WIDTH>& perBGrid, + FsGrid, FS_STENCIL_WIDTH>& EGrid, + FsGrid, FS_STENCIL_WIDTH>& EHallGrid, + FsGrid, FS_STENCIL_WIDTH>& EGradPeGrid, + FsGrid, FS_STENCIL_WIDTH>& momentsGrid, + FsGrid, FS_STENCIL_WIDTH>& dPerBGrid, + FsGrid, FS_STENCIL_WIDTH>& dMomentsGrid, + FsGrid, FS_STENCIL_WIDTH>& BgBGrid, + FsGrid& technicalGrid, + cint i, + cint j, + cint k, + cint& RKCase) { +#ifdef DEBUG_FSOLVER bool ok = true; - if (technicalGrid.get(i,j,k) == NULL) ok = false; - if (technicalGrid.get(i,j,k-1) == NULL) ok = false; - if (technicalGrid.get(i-1,j,k) == NULL) ok = false; - if (technicalGrid.get(i-1,j,k-1) == NULL) ok = false; + if (technicalGrid.get(i, j, k) == NULL) + ok = false; + if (technicalGrid.get(i, j, k - 1) == NULL) + ok = false; + if (technicalGrid.get(i - 1, j, k) == NULL) + ok = false; + if (technicalGrid.get(i - 1, j, k - 1) == NULL) + ok = false; if (ok == false) { cerr << "NULL pointer in " << __FILE__ << ":" << __LINE__ << std::endl; exit(1); } - #endif - +#endif + // An edge has four neighbouring spatial cells. Calculate // electric field in each of the four cells per edge. - Real ax_pos,ax_neg; // Max. characteristic velocities to x-direction - Real az_pos,az_neg; // Max. characteristic velocities to z-direction - Real Vx0,Vz0; // Reconstructed V - Real vA, vS, vW; // Alfven, sound, whistler speed - Real maxV = 0.0; // Max velocity for CFL purposes - Real c_x,c_z; // Wave speeds to xz-directions - - std::array * perb_SW = perBGrid.get(i ,j ,k ); - std::array * perb_SE = perBGrid.get(i ,j ,k-1); - std::array * perb_NW = perBGrid.get(i-1,j ,k ); - std::array * perb_NE = perBGrid.get(i-1,j ,k-1); - std::array * bgb_SW = BgBGrid.get(i ,j ,k ); - std::array * bgb_SE = BgBGrid.get(i ,j ,k-1); - std::array * bgb_NW = BgBGrid.get(i-1,j ,k ); - std::array * bgb_NE = BgBGrid.get(i-1,j ,k-1); - std::array * moments_SW = momentsGrid.get(i ,j ,k ); - std::array * moments_SE = momentsGrid.get(i ,j ,k-1); - std::array * moments_NW = momentsGrid.get(i-1,j ,k ); - std::array * moments_NE = momentsGrid.get(i-1,j ,k-1); - std::array * dmoments_SW = dMomentsGrid.get(i ,j ,k ); - std::array * dmoments_SE = dMomentsGrid.get(i ,j ,k-1); - std::array * dmoments_NW = dMomentsGrid.get(i-1,j ,k ); - std::array * dmoments_NE = dMomentsGrid.get(i-1,j ,k-1); - std::array * dperb_SW = dPerBGrid.get(i ,j ,k ); - std::array * dperb_SE = dPerBGrid.get(i ,j ,k-1); - std::array * dperb_NW = dPerBGrid.get(i-1,j ,k ); - std::array * dperb_NE = dPerBGrid.get(i-1,j ,k-1); - - std::array * efield_SW = EGrid.get(i,j,k); - + Real ax_pos, ax_neg; // Max. characteristic velocities to x-direction + Real az_pos, az_neg; // Max. characteristic velocities to z-direction + Real Vx0, Vz0; // Reconstructed V + Real vA, vS, vW; // Alfven, sound, whistler speed + Real maxV = 0.0; // Max velocity for CFL purposes + Real c_x, c_z; // Wave speeds to xz-directions + + std::array* perb_SW = perBGrid.get(i, j, k); + std::array* perb_SE = perBGrid.get(i, j, k - 1); + std::array* perb_NW = perBGrid.get(i - 1, j, k); + std::array* perb_NE = perBGrid.get(i - 1, j, k - 1); + std::array* bgb_SW = BgBGrid.get(i, j, k); + std::array* bgb_SE = BgBGrid.get(i, j, k - 1); + std::array* bgb_NW = BgBGrid.get(i - 1, j, k); + std::array* bgb_NE = BgBGrid.get(i - 1, j, k - 1); + std::array* moments_SW = momentsGrid.get(i, j, k); + std::array* moments_SE = momentsGrid.get(i, j, k - 1); + std::array* moments_NW = momentsGrid.get(i - 1, j, k); + std::array* moments_NE = momentsGrid.get(i - 1, j, k - 1); + std::array* dmoments_SW = dMomentsGrid.get(i, j, k); + std::array* dmoments_SE = dMomentsGrid.get(i, j, k - 1); + std::array* dmoments_NW = dMomentsGrid.get(i - 1, j, k); + std::array* dmoments_NE = dMomentsGrid.get(i - 1, j, k - 1); + std::array* dperb_SW = dPerBGrid.get(i, j, k); + std::array* dperb_SE = dPerBGrid.get(i, j, k - 1); + std::array* dperb_NW = dPerBGrid.get(i - 1, j, k); + std::array* dperb_NE = dPerBGrid.get(i - 1, j, k - 1); + + std::array* efield_SW = EGrid.get(i, j, k); + // Fetch required plasma parameters: Real Bz_S, Bx_W, Bx_E, Bz_N, perBz_S, perBx_W, perBx_E, perBz_N; Real minRhom = std::numeric_limits::max(); Real maxRhom = std::numeric_limits::min(); - Bz_S = perb_SW->at(fsgrids::bfield::PERBZ)+bgb_SW->at(fsgrids::bgbfield::BGBZ); - Bx_W = perb_SW->at(fsgrids::bfield::PERBX)+bgb_SW->at(fsgrids::bgbfield::BGBX); - Bx_E = perb_SE->at(fsgrids::bfield::PERBX)+bgb_SE->at(fsgrids::bgbfield::BGBX); - Bz_N = perb_NW->at(fsgrids::bfield::PERBZ)+bgb_NW->at(fsgrids::bgbfield::BGBZ); + Bz_S = perb_SW->at(fsgrids::bfield::PERBZ) + bgb_SW->at(fsgrids::bgbfield::BGBZ); + Bx_W = perb_SW->at(fsgrids::bfield::PERBX) + bgb_SW->at(fsgrids::bgbfield::BGBX); + Bx_E = perb_SE->at(fsgrids::bfield::PERBX) + bgb_SE->at(fsgrids::bgbfield::BGBX); + Bz_N = perb_NW->at(fsgrids::bfield::PERBZ) + bgb_NW->at(fsgrids::bgbfield::BGBZ); perBz_S = perb_SW->at(fsgrids::bfield::PERBZ); perBx_W = perb_SW->at(fsgrids::bfield::PERBX); perBx_E = perb_SE->at(fsgrids::bfield::PERBX); perBz_N = perb_NW->at(fsgrids::bfield::PERBZ); - Vx0 = moments_SW->at(fsgrids::moments::VX); - Vz0 = moments_SW->at(fsgrids::moments::VZ); + Vx0 = moments_SW->at(fsgrids::moments::VX); + Vz0 = moments_SW->at(fsgrids::moments::VZ); minRhom = min(minRhom, - min(moments_SW->at(fsgrids::moments::RHOM), - min(moments_SE->at(fsgrids::moments::RHOM), - min(moments_NW->at(fsgrids::moments::RHOM), - moments_NE->at(fsgrids::moments::RHOM)) - ) - ) - ); + min(moments_SW->at(fsgrids::moments::RHOM), + min(moments_SE->at(fsgrids::moments::RHOM), + min(moments_NW->at(fsgrids::moments::RHOM), moments_NE->at(fsgrids::moments::RHOM))))); maxRhom = max(maxRhom, - max(moments_SW->at(fsgrids::moments::RHOM), - max(moments_SE->at(fsgrids::moments::RHOM), - max(moments_NW->at(fsgrids::moments::RHOM), - moments_NE->at(fsgrids::moments::RHOM)) - ) - ) - ); - + max(moments_SW->at(fsgrids::moments::RHOM), + max(moments_SE->at(fsgrids::moments::RHOM), + max(moments_NW->at(fsgrids::moments::RHOM), moments_NE->at(fsgrids::moments::RHOM))))); + creal dBxdy_W = dperb_SW->at(fsgrids::dperb::dPERBxdy) + bgb_SW->at(fsgrids::bgbfield::dBGBxdy); creal dBxdz_W = dperb_SW->at(fsgrids::dperb::dPERBxdz) + bgb_SW->at(fsgrids::bgbfield::dBGBxdz); creal dBzdx_S = dperb_SW->at(fsgrids::dperb::dPERBzdx) + bgb_SW->at(fsgrids::bgbfield::dBGBzdx); @@ -911,353 +1020,424 @@ void calculateEdgeElectricFieldY( creal dperBzdx_N = dperb_NW->at(fsgrids::dperb::dPERBzdx); creal dperBxdz_W = dperb_SW->at(fsgrids::dperb::dPERBxdz); creal dperBxdz_E = dperb_SE->at(fsgrids::dperb::dPERBxdz); - + // Ey and characteristic speeds on this cell: // 1st order terms: - Real Ey_SW = Bz_S*Vx0 - Bx_W*Vz0; - + Real Ey_SW = Bz_S * Vx0 - Bx_W * Vz0; + // Resistive term if (Parameters::resistivity > 0) { Ey_SW += Parameters::resistivity * - sqrt((bgb_SW->at(fsgrids::bgbfield::BGBX)+perb_SW->at(fsgrids::bfield::PERBX))* - (bgb_SW->at(fsgrids::bgbfield::BGBX)+perb_SW->at(fsgrids::bfield::PERBX)) + - (bgb_SW->at(fsgrids::bgbfield::BGBY)+perb_SW->at(fsgrids::bfield::PERBY))* - (bgb_SW->at(fsgrids::bgbfield::BGBY)+perb_SW->at(fsgrids::bfield::PERBY)) + - (bgb_SW->at(fsgrids::bgbfield::BGBZ)+perb_SW->at(fsgrids::bfield::PERBZ))* - (bgb_SW->at(fsgrids::bgbfield::BGBZ)+perb_SW->at(fsgrids::bfield::PERBZ)) - ) / - moments_SW->at(fsgrids::moments::RHOQ) / - physicalconstants::MU_0 * - (dperb_SW->at(fsgrids::dperb::dPERBxdz)/technicalGrid.DZ - dperb_SW->at(fsgrids::dperb::dPERBzdx)/technicalGrid.DX); + sqrt((bgb_SW->at(fsgrids::bgbfield::BGBX) + perb_SW->at(fsgrids::bfield::PERBX)) * + (bgb_SW->at(fsgrids::bgbfield::BGBX) + perb_SW->at(fsgrids::bfield::PERBX)) + + (bgb_SW->at(fsgrids::bgbfield::BGBY) + perb_SW->at(fsgrids::bfield::PERBY)) * + (bgb_SW->at(fsgrids::bgbfield::BGBY) + perb_SW->at(fsgrids::bfield::PERBY)) + + (bgb_SW->at(fsgrids::bgbfield::BGBZ) + perb_SW->at(fsgrids::bfield::PERBZ)) * + (bgb_SW->at(fsgrids::bgbfield::BGBZ) + perb_SW->at(fsgrids::bfield::PERBZ))) / + moments_SW->at(fsgrids::moments::RHOQ) / physicalconstants::MU_0 * + (dperb_SW->at(fsgrids::dperb::dPERBxdz) / technicalGrid.DZ - + dperb_SW->at(fsgrids::dperb::dPERBzdx) / technicalGrid.DX); } // Hall term if (Parameters::ohmHallTerm > 0) { - Ey_SW += EHallGrid.get(i,j,k)->at(fsgrids::ehall::EYHALL_000_010); + Ey_SW += EHallGrid.get(i, j, k)->at(fsgrids::ehall::EYHALL_000_010); } - + // Electron pressure gradient term - if(Parameters::ohmGradPeTerm > 0) { - Ey_SW += EGradPeGrid.get(i,j,k)->at(fsgrids::egradpe::EYGRADPE); + if (Parameters::ohmGradPeTerm > 0) { + Ey_SW += EGradPeGrid.get(i, j, k)->at(fsgrids::egradpe::EYGRADPE); } - - #ifndef FS_1ST_ORDER_SPACE - // 2nd order terms - Ey_SW += +HALF*((Bz_S - HALF*dBzdx_S)*(-dmoments_SW->at(fsgrids::dmoments::dVxdx) - dmoments_SW->at(fsgrids::dmoments::dVxdz)) - dBzdx_S*Vx0 + SIXTH*dBzdy_S*dmoments_SW->at(fsgrids::dmoments::dVxdy)); - Ey_SW += -HALF*((Bx_W - HALF*dBxdz_W)*(-dmoments_SW->at(fsgrids::dmoments::dVzdx) - dmoments_SW->at(fsgrids::dmoments::dVzdz)) - dBxdz_W*Vz0 + SIXTH*dBxdy_W*dmoments_SW->at(fsgrids::dmoments::dVzdy)); - #endif - - calculateWaveSpeedXZ( - perBGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - i, j, k, - i, j+1, k, - Bx_W, Bz_S, dBxdy_W, dBxdz_W, dBzdx_S, dBzdy_S, MINUS, MINUS, minRhom, maxRhom, vA, vS, vW - ); - c_z = min(Parameters::maxWaveVelocity,sqrt(vA*vA + vS*vS) + vW); + +#ifndef FS_1ST_ORDER_SPACE + // 2nd order terms + Ey_SW += +HALF * ((Bz_S - HALF * dBzdx_S) * + (-dmoments_SW->at(fsgrids::dmoments::dVxdx) - dmoments_SW->at(fsgrids::dmoments::dVxdz)) - + dBzdx_S * Vx0 + SIXTH * dBzdy_S * dmoments_SW->at(fsgrids::dmoments::dVxdy)); + Ey_SW += -HALF * ((Bx_W - HALF * dBxdz_W) * + (-dmoments_SW->at(fsgrids::dmoments::dVzdx) - dmoments_SW->at(fsgrids::dmoments::dVzdz)) - + dBxdz_W * Vz0 + SIXTH * dBxdy_W * dmoments_SW->at(fsgrids::dmoments::dVzdy)); +#endif + + calculateWaveSpeedXZ(perBGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + i, + j, + k, + i, + j + 1, + k, + Bx_W, + Bz_S, + dBxdy_W, + dBxdz_W, + dBzdx_S, + dBzdy_S, + MINUS, + MINUS, + minRhom, + maxRhom, + vA, + vS, + vW); + c_z = min(Parameters::maxWaveVelocity, sqrt(vA * vA + vS * vS) + vW); c_x = c_z; - az_neg = max(ZERO,-Vz0 + c_z); - az_pos = max(ZERO,+Vz0 + c_z); - ax_neg = max(ZERO,-Vx0 + c_x); - ax_pos = max(ZERO,+Vx0 + c_x); + az_neg = max(ZERO, -Vz0 + c_z); + az_pos = max(ZERO, +Vz0 + c_z); + ax_neg = max(ZERO, -Vx0 + c_x); + ax_pos = max(ZERO, +Vx0 + c_x); maxV = max(maxV, calculateCflSpeed(Vz0, Vx0, vA, vS, vW)); // Ey and characteristic speeds on k-1 neighbour: - Vx0 = moments_SE->at(fsgrids::moments::VX); - Vz0 = moments_SE->at(fsgrids::moments::VZ); + Vx0 = moments_SE->at(fsgrids::moments::VX); + Vz0 = moments_SE->at(fsgrids::moments::VZ); // 1st order terms: - Real Ey_SE = Bz_S*Vx0 - Bx_E*Vz0; + Real Ey_SE = Bz_S * Vx0 - Bx_E * Vz0; // Resistive term if (Parameters::resistivity > 0) { Ey_SE += Parameters::resistivity * - sqrt((bgb_SE->at(fsgrids::bgbfield::BGBX)+perb_SE->at(fsgrids::bfield::PERBX))* - (bgb_SE->at(fsgrids::bgbfield::BGBX)+perb_SE->at(fsgrids::bfield::PERBX)) + - (bgb_SE->at(fsgrids::bgbfield::BGBY)+perb_SE->at(fsgrids::bfield::PERBY))* - (bgb_SE->at(fsgrids::bgbfield::BGBY)+perb_SE->at(fsgrids::bfield::PERBY)) + - (bgb_SE->at(fsgrids::bgbfield::BGBZ)+perb_SE->at(fsgrids::bfield::PERBZ))* - (bgb_SE->at(fsgrids::bgbfield::BGBZ)+perb_SE->at(fsgrids::bfield::PERBZ)) - ) / - moments_SE->at(fsgrids::moments::RHOQ) / - physicalconstants::MU_0 * - (dperb_SE->at(fsgrids::dperb::dPERBxdz)/technicalGrid.DZ - dperb_SE->at(fsgrids::dperb::dPERBzdx)/technicalGrid.DX); + sqrt((bgb_SE->at(fsgrids::bgbfield::BGBX) + perb_SE->at(fsgrids::bfield::PERBX)) * + (bgb_SE->at(fsgrids::bgbfield::BGBX) + perb_SE->at(fsgrids::bfield::PERBX)) + + (bgb_SE->at(fsgrids::bgbfield::BGBY) + perb_SE->at(fsgrids::bfield::PERBY)) * + (bgb_SE->at(fsgrids::bgbfield::BGBY) + perb_SE->at(fsgrids::bfield::PERBY)) + + (bgb_SE->at(fsgrids::bgbfield::BGBZ) + perb_SE->at(fsgrids::bfield::PERBZ)) * + (bgb_SE->at(fsgrids::bgbfield::BGBZ) + perb_SE->at(fsgrids::bfield::PERBZ))) / + moments_SE->at(fsgrids::moments::RHOQ) / physicalconstants::MU_0 * + (dperb_SE->at(fsgrids::dperb::dPERBxdz) / technicalGrid.DZ - + dperb_SE->at(fsgrids::dperb::dPERBzdx) / technicalGrid.DX); } // Hall term if (Parameters::ohmHallTerm > 0) { - Ey_SE += EHallGrid.get(i,j,k-1)->at(fsgrids::ehall::EYHALL_001_011); + Ey_SE += EHallGrid.get(i, j, k - 1)->at(fsgrids::ehall::EYHALL_001_011); } - + // Electron pressure gradient term - if(Parameters::ohmGradPeTerm > 0) { - Ey_SE += EGradPeGrid.get(i,j,k-1)->at(fsgrids::egradpe::EYGRADPE); + if (Parameters::ohmGradPeTerm > 0) { + Ey_SE += EGradPeGrid.get(i, j, k - 1)->at(fsgrids::egradpe::EYGRADPE); } - - #ifndef FS_1ST_ORDER_SPACE - // 2nd order terms: - Ey_SE += +HALF*((Bz_S - HALF*dBzdx_S)*(-dmoments_SE->at(fsgrids::dmoments::dVxdx) + dmoments_SE->at(fsgrids::dmoments::dVxdz)) - dBzdx_S*Vx0 + SIXTH*dBzdy_S*dmoments_SE->at(fsgrids::dmoments::dVxdy)); - Ey_SE += -HALF*((Bx_E + HALF*dBxdz_E)*(-dmoments_SE->at(fsgrids::dmoments::dVzdx) + dmoments_SE->at(fsgrids::dmoments::dVzdz)) + dBxdz_E*Vz0 + SIXTH*dBxdy_E*dmoments_SE->at(fsgrids::dmoments::dVzdy)); - #endif - - calculateWaveSpeedXZ( - perBGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - i, j , k-1, - i, j+1, k-1, - Bx_E, Bz_S, dBxdy_E, dBxdz_E, dBzdx_S, dBzdy_S, MINUS, PLUS, minRhom, maxRhom, vA, vS, vW - ); - c_z = min(Parameters::maxWaveVelocity,sqrt(vA*vA + vS*vS) + vW); + +#ifndef FS_1ST_ORDER_SPACE + // 2nd order terms: + Ey_SE += +HALF * ((Bz_S - HALF * dBzdx_S) * + (-dmoments_SE->at(fsgrids::dmoments::dVxdx) + dmoments_SE->at(fsgrids::dmoments::dVxdz)) - + dBzdx_S * Vx0 + SIXTH * dBzdy_S * dmoments_SE->at(fsgrids::dmoments::dVxdy)); + Ey_SE += -HALF * ((Bx_E + HALF * dBxdz_E) * + (-dmoments_SE->at(fsgrids::dmoments::dVzdx) + dmoments_SE->at(fsgrids::dmoments::dVzdz)) + + dBxdz_E * Vz0 + SIXTH * dBxdy_E * dmoments_SE->at(fsgrids::dmoments::dVzdy)); +#endif + + calculateWaveSpeedXZ(perBGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + i, + j, + k - 1, + i, + j + 1, + k - 1, + Bx_E, + Bz_S, + dBxdy_E, + dBxdz_E, + dBzdx_S, + dBzdy_S, + MINUS, + PLUS, + minRhom, + maxRhom, + vA, + vS, + vW); + c_z = min(Parameters::maxWaveVelocity, sqrt(vA * vA + vS * vS) + vW); c_x = c_z; - az_neg = max(az_neg,-Vz0 + c_z); - az_pos = max(az_pos,+Vz0 + c_z); - ax_neg = max(ax_neg,-Vx0 + c_x); - ax_pos = max(ax_pos,+Vx0 + c_x); + az_neg = max(az_neg, -Vz0 + c_z); + az_pos = max(az_pos, +Vz0 + c_z); + ax_neg = max(ax_neg, -Vx0 + c_x); + ax_pos = max(ax_pos, +Vx0 + c_x); maxV = max(maxV, calculateCflSpeed(Vz0, Vx0, vA, vS, vW)); - + // Ey and characteristic speeds on i-1 neighbour: - Vz0 = moments_NW->at(fsgrids::moments::VZ); - Vx0 = moments_NW->at(fsgrids::moments::VX); - + Vz0 = moments_NW->at(fsgrids::moments::VZ); + Vx0 = moments_NW->at(fsgrids::moments::VX); + // 1st order terms: - Real Ey_NW = Bz_N*Vx0 - Bx_W*Vz0; + Real Ey_NW = Bz_N * Vx0 - Bx_W * Vz0; // Resistive term if (Parameters::resistivity > 0) { Ey_NW += Parameters::resistivity * - sqrt((bgb_NW->at(fsgrids::bgbfield::BGBX)+perb_NW->at(fsgrids::bfield::PERBX))* - (bgb_NW->at(fsgrids::bgbfield::BGBX)+perb_NW->at(fsgrids::bfield::PERBX)) + - (bgb_NW->at(fsgrids::bgbfield::BGBY)+perb_NW->at(fsgrids::bfield::PERBY))* - (bgb_NW->at(fsgrids::bgbfield::BGBY)+perb_NW->at(fsgrids::bfield::PERBY)) + - (bgb_NW->at(fsgrids::bgbfield::BGBZ)+perb_NW->at(fsgrids::bfield::PERBZ))* - (bgb_NW->at(fsgrids::bgbfield::BGBZ)+perb_NW->at(fsgrids::bfield::PERBZ)) - ) / - moments_NW->at(fsgrids::moments::RHOQ) / - physicalconstants::MU_0 * - (dperb_NW->at(fsgrids::dperb::dPERBxdz)/technicalGrid.DZ - dperb_NW->at(fsgrids::dperb::dPERBzdx)/technicalGrid.DX); + sqrt((bgb_NW->at(fsgrids::bgbfield::BGBX) + perb_NW->at(fsgrids::bfield::PERBX)) * + (bgb_NW->at(fsgrids::bgbfield::BGBX) + perb_NW->at(fsgrids::bfield::PERBX)) + + (bgb_NW->at(fsgrids::bgbfield::BGBY) + perb_NW->at(fsgrids::bfield::PERBY)) * + (bgb_NW->at(fsgrids::bgbfield::BGBY) + perb_NW->at(fsgrids::bfield::PERBY)) + + (bgb_NW->at(fsgrids::bgbfield::BGBZ) + perb_NW->at(fsgrids::bfield::PERBZ)) * + (bgb_NW->at(fsgrids::bgbfield::BGBZ) + perb_NW->at(fsgrids::bfield::PERBZ))) / + moments_NW->at(fsgrids::moments::RHOQ) / physicalconstants::MU_0 * + (dperb_NW->at(fsgrids::dperb::dPERBxdz) / technicalGrid.DZ - + dperb_NW->at(fsgrids::dperb::dPERBzdx) / technicalGrid.DX); } // Hall term - if(Parameters::ohmHallTerm > 0) { - Ey_NW += EHallGrid.get(i-1,j,k)->at(fsgrids::ehall::EYHALL_100_110); + if (Parameters::ohmHallTerm > 0) { + Ey_NW += EHallGrid.get(i - 1, j, k)->at(fsgrids::ehall::EYHALL_100_110); } - + // Electron pressure gradient term - if(Parameters::ohmGradPeTerm > 0) { - Ey_NW += EGradPeGrid.get(i-1,j,k)->at(fsgrids::egradpe::EYGRADPE); + if (Parameters::ohmGradPeTerm > 0) { + Ey_NW += EGradPeGrid.get(i - 1, j, k)->at(fsgrids::egradpe::EYGRADPE); } - - #ifndef FS_1ST_ORDER_SPACE - // 2nd order terms: - Ey_NW += +HALF*((Bz_N + HALF*dBzdx_N)*(+dmoments_NW->at(fsgrids::dmoments::dVxdx) - dmoments_NW->at(fsgrids::dmoments::dVxdz)) + dBzdx_N*Vx0 + SIXTH*dBzdy_N*dmoments_NW->at(fsgrids::dmoments::dVxdy)); - Ey_NW += -HALF*((Bx_W - HALF*dBxdz_W)*(+dmoments_NW->at(fsgrids::dmoments::dVzdx) - dmoments_NW->at(fsgrids::dmoments::dVzdz)) - dBxdz_W*Vz0 + SIXTH*dBxdy_W*dmoments_NW->at(fsgrids::dmoments::dVzdy)); - #endif - - calculateWaveSpeedXZ( - perBGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - i-1,j ,k, - i-1,j+1,k, - Bx_W, Bz_N, dBxdy_W, dBxdz_W, dBzdx_N, dBzdy_N, PLUS, MINUS, minRhom, maxRhom, vA, vS, vW - ); - c_z = min(Parameters::maxWaveVelocity,sqrt(vA*vA + vS*vS) + vW); + +#ifndef FS_1ST_ORDER_SPACE + // 2nd order terms: + Ey_NW += +HALF * ((Bz_N + HALF * dBzdx_N) * + (+dmoments_NW->at(fsgrids::dmoments::dVxdx) - dmoments_NW->at(fsgrids::dmoments::dVxdz)) + + dBzdx_N * Vx0 + SIXTH * dBzdy_N * dmoments_NW->at(fsgrids::dmoments::dVxdy)); + Ey_NW += -HALF * ((Bx_W - HALF * dBxdz_W) * + (+dmoments_NW->at(fsgrids::dmoments::dVzdx) - dmoments_NW->at(fsgrids::dmoments::dVzdz)) - + dBxdz_W * Vz0 + SIXTH * dBxdy_W * dmoments_NW->at(fsgrids::dmoments::dVzdy)); +#endif + + calculateWaveSpeedXZ(perBGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + i - 1, + j, + k, + i - 1, + j + 1, + k, + Bx_W, + Bz_N, + dBxdy_W, + dBxdz_W, + dBzdx_N, + dBzdy_N, + PLUS, + MINUS, + minRhom, + maxRhom, + vA, + vS, + vW); + c_z = min(Parameters::maxWaveVelocity, sqrt(vA * vA + vS * vS) + vW); c_x = c_z; - az_neg = max(az_neg,-Vz0 + c_z); - az_pos = max(az_pos,+Vz0 + c_z); - ax_neg = max(ax_neg,-Vx0 + c_x); - ax_pos = max(ax_pos,+Vx0 + c_x); + az_neg = max(az_neg, -Vz0 + c_z); + az_pos = max(az_pos, +Vz0 + c_z); + ax_neg = max(ax_neg, -Vx0 + c_x); + ax_pos = max(ax_pos, +Vx0 + c_x); maxV = max(maxV, calculateCflSpeed(Vz0, Vx0, vA, vS, vW)); // Ey and characteristic speeds on i-1,k-1 neighbour: Vz0 = moments_NE->at(fsgrids::moments::VZ); Vx0 = moments_NE->at(fsgrids::moments::VX); - + // 1st order terms: - Real Ey_NE = Bz_N*Vx0 - Bx_E*Vz0; - + Real Ey_NE = Bz_N * Vx0 - Bx_E * Vz0; + // Resistive term if (Parameters::resistivity > 0) { Ey_NE += Parameters::resistivity * - sqrt((bgb_NE->at(fsgrids::bgbfield::BGBX)+perb_NE->at(fsgrids::bfield::PERBX))* - (bgb_NE->at(fsgrids::bgbfield::BGBX)+perb_NE->at(fsgrids::bfield::PERBX)) + - (bgb_NE->at(fsgrids::bgbfield::BGBY)+perb_NE->at(fsgrids::bfield::PERBY))* - (bgb_NE->at(fsgrids::bgbfield::BGBY)+perb_NE->at(fsgrids::bfield::PERBY)) + - (bgb_NE->at(fsgrids::bgbfield::BGBZ)+perb_NE->at(fsgrids::bfield::PERBZ))* - (bgb_NE->at(fsgrids::bgbfield::BGBZ)+perb_NE->at(fsgrids::bfield::PERBZ)) - ) / - moments_NE->at(fsgrids::moments::RHOQ) / - physicalconstants::MU_0 * - (dperb_NE->at(fsgrids::dperb::dPERBxdz)/technicalGrid.DZ - dperb_NE->at(fsgrids::dperb::dPERBzdx)/technicalGrid.DX); + sqrt((bgb_NE->at(fsgrids::bgbfield::BGBX) + perb_NE->at(fsgrids::bfield::PERBX)) * + (bgb_NE->at(fsgrids::bgbfield::BGBX) + perb_NE->at(fsgrids::bfield::PERBX)) + + (bgb_NE->at(fsgrids::bgbfield::BGBY) + perb_NE->at(fsgrids::bfield::PERBY)) * + (bgb_NE->at(fsgrids::bgbfield::BGBY) + perb_NE->at(fsgrids::bfield::PERBY)) + + (bgb_NE->at(fsgrids::bgbfield::BGBZ) + perb_NE->at(fsgrids::bfield::PERBZ)) * + (bgb_NE->at(fsgrids::bgbfield::BGBZ) + perb_NE->at(fsgrids::bfield::PERBZ))) / + moments_NE->at(fsgrids::moments::RHOQ) / physicalconstants::MU_0 * + (dperb_NE->at(fsgrids::dperb::dPERBxdz) / technicalGrid.DZ - + dperb_NE->at(fsgrids::dperb::dPERBzdx) / technicalGrid.DX); } // Hall term - if(Parameters::ohmHallTerm > 0) { - Ey_NE += EHallGrid.get(i-1,j,k-1)->at(fsgrids::ehall::EYHALL_101_111); + if (Parameters::ohmHallTerm > 0) { + Ey_NE += EHallGrid.get(i - 1, j, k - 1)->at(fsgrids::ehall::EYHALL_101_111); } - + // Electron pressure gradient term - if(Parameters::ohmGradPeTerm > 0) { - Ey_NE += EGradPeGrid.get(i-1,j,k-1)->at(fsgrids::egradpe::EYGRADPE); + if (Parameters::ohmGradPeTerm > 0) { + Ey_NE += EGradPeGrid.get(i - 1, j, k - 1)->at(fsgrids::egradpe::EYGRADPE); } - - #ifndef FS_1ST_ORDER_SPACE - // 2nd order terms: - Ey_NE += +HALF*((Bz_N + HALF*dBzdx_N)*(+dmoments_NE->at(fsgrids::dmoments::dVxdx) + dmoments_NE->at(fsgrids::dmoments::dVxdz)) + dBzdx_N*Vx0 + SIXTH*dBzdy_N*dmoments_NE->at(fsgrids::dmoments::dVxdy)); - Ey_NE += -HALF*((Bx_E + HALF*dBxdz_E)*(+dmoments_NE->at(fsgrids::dmoments::dVzdx) + dmoments_NE->at(fsgrids::dmoments::dVzdz)) + dBxdz_E*Vz0 + SIXTH*dBxdy_E*dmoments_NE->at(fsgrids::dmoments::dVzdy)); - #endif - - calculateWaveSpeedXZ( - perBGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - i-1,j ,k-1, - i-1,j+1,k-1, - Bx_E, Bz_N, dBxdy_E, dBxdz_E, dBzdx_N, dBzdy_N, PLUS, PLUS, minRhom, maxRhom, vA, vS, vW - ); - c_z = min(Parameters::maxWaveVelocity,sqrt(vA*vA + vS*vS) + vW); + +#ifndef FS_1ST_ORDER_SPACE + // 2nd order terms: + Ey_NE += +HALF * ((Bz_N + HALF * dBzdx_N) * + (+dmoments_NE->at(fsgrids::dmoments::dVxdx) + dmoments_NE->at(fsgrids::dmoments::dVxdz)) + + dBzdx_N * Vx0 + SIXTH * dBzdy_N * dmoments_NE->at(fsgrids::dmoments::dVxdy)); + Ey_NE += -HALF * ((Bx_E + HALF * dBxdz_E) * + (+dmoments_NE->at(fsgrids::dmoments::dVzdx) + dmoments_NE->at(fsgrids::dmoments::dVzdz)) + + dBxdz_E * Vz0 + SIXTH * dBxdy_E * dmoments_NE->at(fsgrids::dmoments::dVzdy)); +#endif + + calculateWaveSpeedXZ(perBGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + i - 1, + j, + k - 1, + i - 1, + j + 1, + k - 1, + Bx_E, + Bz_N, + dBxdy_E, + dBxdz_E, + dBzdx_N, + dBzdy_N, + PLUS, + PLUS, + minRhom, + maxRhom, + vA, + vS, + vW); + c_z = min(Parameters::maxWaveVelocity, sqrt(vA * vA + vS * vS) + vW); c_x = c_z; - az_neg = max(az_neg,-Vz0 + c_z); - az_pos = max(az_pos,+Vz0 + c_z); - ax_neg = max(ax_neg,-Vx0 + c_x); - ax_pos = max(ax_pos,+Vx0 + c_x); + az_neg = max(az_neg, -Vz0 + c_z); + az_pos = max(az_pos, +Vz0 + c_z); + ax_neg = max(ax_neg, -Vx0 + c_x); + ax_pos = max(ax_pos, +Vx0 + c_x); maxV = max(maxV, calculateCflSpeed(Vz0, Vx0, vA, vS, vW)); // Calculate properly upwinded edge-averaged Ey: - efield_SW->at(fsgrids::efield::EY) = az_pos*ax_pos*Ey_NE + az_pos*ax_neg*Ey_SE + az_neg*ax_pos*Ey_NW + az_neg*ax_neg*Ey_SW; - efield_SW->at(fsgrids::efield::EY) /= ((az_pos+az_neg)*(ax_pos+ax_neg)+EPS); + efield_SW->at(fsgrids::efield::EY) = + az_pos * ax_pos * Ey_NE + az_pos * ax_neg * Ey_SE + az_neg * ax_pos * Ey_NW + az_neg * ax_neg * Ey_SW; + efield_SW->at(fsgrids::efield::EY) /= ((az_pos + az_neg) * (ax_pos + ax_neg) + EPS); if (Parameters::fieldSolverDiffusiveEterms) { #ifdef FS_1ST_ORDER_SPACE - efield_SW->at(fsgrids::efield::EY) -= ax_pos*ax_neg/(ax_pos+ax_neg+EPS)*(perBz_S-perBz_N); - efield_SW->at(fsgrids::efield::EY) += az_pos*az_neg/(az_pos+az_neg+EPS)*(perBx_W-perBx_E); + efield_SW->at(fsgrids::efield::EY) -= ax_pos * ax_neg / (ax_pos + ax_neg + EPS) * (perBz_S - perBz_N); + efield_SW->at(fsgrids::efield::EY) += az_pos * az_neg / (az_pos + az_neg + EPS) * (perBx_W - perBx_E); #else - efield_SW->at(fsgrids::efield::EY) -= ax_pos*ax_neg/(ax_pos+ax_neg+EPS)*((perBz_S-HALF*dperBzdx_S) - (perBz_N+HALF*dperBzdx_N)); - efield_SW->at(fsgrids::efield::EY) += az_pos*az_neg/(az_pos+az_neg+EPS)*((perBx_W-HALF*dperBxdz_W) - (perBx_E+HALF*dperBxdz_E)); + efield_SW->at(fsgrids::efield::EY) -= + ax_pos * ax_neg / (ax_pos + ax_neg + EPS) * ((perBz_S - HALF * dperBzdx_S) - (perBz_N + HALF * dperBzdx_N)); + efield_SW->at(fsgrids::efield::EY) += + az_pos * az_neg / (az_pos + az_neg + EPS) * ((perBx_W - HALF * dperBxdz_W) - (perBx_E + HALF * dperBxdz_E)); #endif } - + if ((RKCase == RK_ORDER1) || (RKCase == RK_ORDER2_STEP2)) { - //compute maximum timestep for fieldsolver in this cell (CFL=1) - Real min_dx=std::numeric_limits::max();; - min_dx=min(min_dx,technicalGrid.DX); - min_dx=min(min_dx,technicalGrid.DZ); - //update max allowed timestep for field propagation in this cell, which is the minimum of CFL=1 timesteps - if (maxV!=ZERO) technicalGrid.get(i,j,k)->maxFsDt=min(technicalGrid.get(i,j,k)->maxFsDt,min_dx/maxV); + // compute maximum timestep for fieldsolver in this cell (CFL=1) + Real min_dx = std::numeric_limits::max(); + ; + min_dx = min(min_dx, technicalGrid.DX); + min_dx = min(min_dx, technicalGrid.DZ); + // update max allowed timestep for field propagation in this cell, which is the minimum of CFL=1 timesteps + if (maxV != ZERO) + technicalGrid.get(i, j, k)->maxFsDt = min(technicalGrid.get(i, j, k)->maxFsDt, min_dx / maxV); } } /*! \brief Low-level electric field propagation function. * - * Computes the upwinded electric field Z component along the cell's corresponding edge as the cross product of B and V in the XY plane. Also includes the calculation of the maximally allowed time step. - * - * Expects that the correct RHO and B fields are being passed, depending on the stage of the Runge-Kutta time stepping method. - * - * Note that the background B field is excluded from the diffusive term calculations because they are equivalent to a current term and the background field is curl-free. - * + * Computes the upwinded electric field Z component along the cell's corresponding edge as the cross product of B and V + * in the XY plane. Also includes the calculation of the maximally allowed time step. + * + * Expects that the correct RHO and B fields are being passed, depending on the stage of the Runge-Kutta time stepping + * method. + * + * Note that the background B field is excluded from the diffusive term calculations because they are equivalent to a + * current term and the background field is curl-free. + * * \param RKCase Element in the enum defining the Runge-Kutta method steps */ void calculateEdgeElectricFieldZ( - FsGrid< std::array, FS_STENCIL_WIDTH> & perBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EHallGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, - FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, - cint i, - cint j, - cint k, - cint& RKCase -) { - #ifdef DEBUG_FSOLVER + FsGrid, FS_STENCIL_WIDTH>& perBGrid, + FsGrid, FS_STENCIL_WIDTH>& EGrid, + FsGrid, FS_STENCIL_WIDTH>& EHallGrid, + FsGrid, FS_STENCIL_WIDTH>& EGradPeGrid, + FsGrid, FS_STENCIL_WIDTH>& momentsGrid, + FsGrid, FS_STENCIL_WIDTH>& dPerBGrid, + FsGrid, FS_STENCIL_WIDTH>& dMomentsGrid, + FsGrid, FS_STENCIL_WIDTH>& BgBGrid, + FsGrid& technicalGrid, + cint i, + cint j, + cint k, + cint& RKCase) { +#ifdef DEBUG_FSOLVER bool ok = true; - if (technicalGrid.get(i,j,k) == NULL) ok = false; - if (technicalGrid.get(i-1,j,k) == NULL) ok = false; - if (technicalGrid.get(i-1,j-1,k) == NULL) ok = false; - if (technicalGrid.get(i,j-1,k) == NULL) ok = false; + if (technicalGrid.get(i, j, k) == NULL) + ok = false; + if (technicalGrid.get(i - 1, j, k) == NULL) + ok = false; + if (technicalGrid.get(i - 1, j - 1, k) == NULL) + ok = false; + if (technicalGrid.get(i, j - 1, k) == NULL) + ok = false; if (ok == false) { cerr << "NULL pointer in " << __FILE__ << ":" << __LINE__ << std::endl; exit(1); } - #endif +#endif - // An edge has four neighbouring spatial cells. Calculate + // An edge has four neighbouring spatial cells. Calculate // electric field in each of the four cells per edge. - Real ax_pos,ax_neg; // Max. characteristic velocities to x-direction - Real ay_pos,ay_neg; // Max. characteristic velocities to y-direction - Real Vx0,Vy0; // Reconstructed V - Real vA, vS, vW; // Alfven, sound, whistler speed - Real maxV = 0.0; // Max velocity for CFL purposes - Real c_x,c_y; // Characteristic speeds to xy-directions - + Real ax_pos, ax_neg; // Max. characteristic velocities to x-direction + Real ay_pos, ay_neg; // Max. characteristic velocities to y-direction + Real Vx0, Vy0; // Reconstructed V + Real vA, vS, vW; // Alfven, sound, whistler speed + Real maxV = 0.0; // Max velocity for CFL purposes + Real c_x, c_y; // Characteristic speeds to xy-directions + // Get read-only pointers to NE,NW,SE,SW states (SW is rw, result is written there): - std::array * perb_SW = perBGrid.get(i ,j ,k ); - std::array * perb_SE = perBGrid.get(i-1,j ,k ); - std::array * perb_NE = perBGrid.get(i-1,j-1,k ); - std::array * perb_NW = perBGrid.get(i ,j-1,k ); - std::array * bgb_SW = BgBGrid.get(i ,j ,k ); - std::array * bgb_SE = BgBGrid.get(i-1,j ,k ); - std::array * bgb_NE = BgBGrid.get(i-1,j-1,k ); - std::array * bgb_NW = BgBGrid.get(i ,j-1,k ); - std::array * moments_SW = momentsGrid.get(i ,j ,k ); - std::array * moments_SE = momentsGrid.get(i-1,j ,k ); - std::array * moments_NE = momentsGrid.get(i-1,j-1,k ); - std::array * moments_NW = momentsGrid.get(i ,j-1,k ); - std::array * dmoments_SW = dMomentsGrid.get(i ,j ,k ); - std::array * dmoments_SE = dMomentsGrid.get(i-1,j ,k ); - std::array * dmoments_NE = dMomentsGrid.get(i-1,j-1,k ); - std::array * dmoments_NW = dMomentsGrid.get(i ,j-1,k ); - std::array * dperb_SW = dPerBGrid.get(i ,j ,k ); - std::array * dperb_SE = dPerBGrid.get(i-1,j ,k ); - std::array * dperb_NE = dPerBGrid.get(i-1,j-1,k ); - std::array * dperb_NW = dPerBGrid.get(i ,j-1,k ); - - std::array * efield_SW = EGrid.get(i,j,k); - + std::array* perb_SW = perBGrid.get(i, j, k); + std::array* perb_SE = perBGrid.get(i - 1, j, k); + std::array* perb_NE = perBGrid.get(i - 1, j - 1, k); + std::array* perb_NW = perBGrid.get(i, j - 1, k); + std::array* bgb_SW = BgBGrid.get(i, j, k); + std::array* bgb_SE = BgBGrid.get(i - 1, j, k); + std::array* bgb_NE = BgBGrid.get(i - 1, j - 1, k); + std::array* bgb_NW = BgBGrid.get(i, j - 1, k); + std::array* moments_SW = momentsGrid.get(i, j, k); + std::array* moments_SE = momentsGrid.get(i - 1, j, k); + std::array* moments_NE = momentsGrid.get(i - 1, j - 1, k); + std::array* moments_NW = momentsGrid.get(i, j - 1, k); + std::array* dmoments_SW = dMomentsGrid.get(i, j, k); + std::array* dmoments_SE = dMomentsGrid.get(i - 1, j, k); + std::array* dmoments_NE = dMomentsGrid.get(i - 1, j - 1, k); + std::array* dmoments_NW = dMomentsGrid.get(i, j - 1, k); + std::array* dperb_SW = dPerBGrid.get(i, j, k); + std::array* dperb_SE = dPerBGrid.get(i - 1, j, k); + std::array* dperb_NE = dPerBGrid.get(i - 1, j - 1, k); + std::array* dperb_NW = dPerBGrid.get(i, j - 1, k); + + std::array* efield_SW = EGrid.get(i, j, k); + // Fetch needed plasma parameters/derivatives from the four cells: Real Bx_S, By_W, By_E, Bx_N, perBx_S, perBy_W, perBy_E, perBx_N; Real minRhom = std::numeric_limits::max(); Real maxRhom = std::numeric_limits::min(); - Bx_S = perb_SW->at(fsgrids::bfield::PERBX) + bgb_SW->at(fsgrids::bgbfield::BGBX); - By_W = perb_SW->at(fsgrids::bfield::PERBY) + bgb_SW->at(fsgrids::bgbfield::BGBY); - By_E = perb_SE->at(fsgrids::bfield::PERBY) + bgb_SE->at(fsgrids::bgbfield::BGBY); - Bx_N = perb_NW->at(fsgrids::bfield::PERBX) + bgb_NW->at(fsgrids::bgbfield::BGBX); - perBx_S = perb_SW->at(fsgrids::bfield::PERBX); - perBy_W = perb_SW->at(fsgrids::bfield::PERBY); - perBy_E = perb_SE->at(fsgrids::bfield::PERBY); - perBx_N = perb_NW->at(fsgrids::bfield::PERBX); - Vx0 = moments_SW->at(fsgrids::moments::VX); - Vy0 = moments_SW->at(fsgrids::moments::VY); + Bx_S = perb_SW->at(fsgrids::bfield::PERBX) + bgb_SW->at(fsgrids::bgbfield::BGBX); + By_W = perb_SW->at(fsgrids::bfield::PERBY) + bgb_SW->at(fsgrids::bgbfield::BGBY); + By_E = perb_SE->at(fsgrids::bfield::PERBY) + bgb_SE->at(fsgrids::bgbfield::BGBY); + Bx_N = perb_NW->at(fsgrids::bfield::PERBX) + bgb_NW->at(fsgrids::bgbfield::BGBX); + perBx_S = perb_SW->at(fsgrids::bfield::PERBX); + perBy_W = perb_SW->at(fsgrids::bfield::PERBY); + perBy_E = perb_SE->at(fsgrids::bfield::PERBY); + perBx_N = perb_NW->at(fsgrids::bfield::PERBX); + Vx0 = moments_SW->at(fsgrids::moments::VX); + Vy0 = moments_SW->at(fsgrids::moments::VY); minRhom = min(minRhom, - min(moments_SW->at(fsgrids::moments::RHOM), - min(moments_SE->at(fsgrids::moments::RHOM), - min(moments_NW->at(fsgrids::moments::RHOM), - moments_NE->at(fsgrids::moments::RHOM)) - ) - ) - ); + min(moments_SW->at(fsgrids::moments::RHOM), + min(moments_SE->at(fsgrids::moments::RHOM), + min(moments_NW->at(fsgrids::moments::RHOM), moments_NE->at(fsgrids::moments::RHOM))))); maxRhom = max(maxRhom, - max(moments_SW->at(fsgrids::moments::RHOM), - max(moments_SE->at(fsgrids::moments::RHOM), - max(moments_NW->at(fsgrids::moments::RHOM), - moments_NE->at(fsgrids::moments::RHOM)) - ) - ) - ); - + max(moments_SW->at(fsgrids::moments::RHOM), + max(moments_SE->at(fsgrids::moments::RHOM), + max(moments_NW->at(fsgrids::moments::RHOM), moments_NE->at(fsgrids::moments::RHOM))))); + creal dBxdy_S = dperb_SW->at(fsgrids::dperb::dPERBxdy) + bgb_SW->at(fsgrids::bgbfield::dBGBxdy); creal dBxdz_S = dperb_SW->at(fsgrids::dperb::dPERBxdz) + bgb_SW->at(fsgrids::bgbfield::dBGBxdz); creal dBydx_W = dperb_SW->at(fsgrids::dperb::dPERBydx) + bgb_SW->at(fsgrids::bgbfield::dBGBydx); @@ -1270,258 +1450,331 @@ void calculateEdgeElectricFieldZ( creal dperBxdy_N = dperb_NW->at(fsgrids::dperb::dPERBxdy); creal dperBydx_W = dperb_SW->at(fsgrids::dperb::dPERBydx); creal dperBydx_E = dperb_SE->at(fsgrids::dperb::dPERBydx); - + // Ez and characteristic speeds on SW cell: // 1st order terms: - Real Ez_SW = Bx_S*Vy0 - By_W*Vx0; - + Real Ez_SW = Bx_S * Vy0 - By_W * Vx0; + // Resistive term if (Parameters::resistivity > 0) { - Ez_SW += Parameters::resistivity * - sqrt((bgb_SW->at(fsgrids::bgbfield::BGBX)+perb_SW->at(fsgrids::bfield::PERBX))* - (bgb_SW->at(fsgrids::bgbfield::BGBX)+perb_SW->at(fsgrids::bfield::PERBX)) + - (bgb_SW->at(fsgrids::bgbfield::BGBY)+perb_SW->at(fsgrids::bfield::PERBY))* - (bgb_SW->at(fsgrids::bgbfield::BGBY)+perb_SW->at(fsgrids::bfield::PERBY)) + - (bgb_SW->at(fsgrids::bgbfield::BGBZ)+perb_SW->at(fsgrids::bfield::PERBZ))* - (bgb_SW->at(fsgrids::bgbfield::BGBZ)+perb_SW->at(fsgrids::bfield::PERBZ)) - ) / - moments_SW->at(fsgrids::moments::RHOQ) / - physicalconstants::MU_0 * - (dperb_SW->at(fsgrids::dperb::dPERBydx)/technicalGrid.DX - dperb_SW->at(fsgrids::dperb::dPERBxdy)/technicalGrid.DY); + Ez_SW += Parameters::resistivity * + sqrt((bgb_SW->at(fsgrids::bgbfield::BGBX) + perb_SW->at(fsgrids::bfield::PERBX)) * + (bgb_SW->at(fsgrids::bgbfield::BGBX) + perb_SW->at(fsgrids::bfield::PERBX)) + + (bgb_SW->at(fsgrids::bgbfield::BGBY) + perb_SW->at(fsgrids::bfield::PERBY)) * + (bgb_SW->at(fsgrids::bgbfield::BGBY) + perb_SW->at(fsgrids::bfield::PERBY)) + + (bgb_SW->at(fsgrids::bgbfield::BGBZ) + perb_SW->at(fsgrids::bfield::PERBZ)) * + (bgb_SW->at(fsgrids::bgbfield::BGBZ) + perb_SW->at(fsgrids::bfield::PERBZ))) / + moments_SW->at(fsgrids::moments::RHOQ) / physicalconstants::MU_0 * + (dperb_SW->at(fsgrids::dperb::dPERBydx) / technicalGrid.DX - + dperb_SW->at(fsgrids::dperb::dPERBxdy) / technicalGrid.DY); } - + // Hall term if (Parameters::ohmHallTerm > 0) { - Ez_SW += EHallGrid.get(i,j,k)->at(fsgrids::ehall::EZHALL_000_001); + Ez_SW += EHallGrid.get(i, j, k)->at(fsgrids::ehall::EZHALL_000_001); } - + // Electron pressure gradient term - if(Parameters::ohmGradPeTerm > 0) { - Ez_SW += EGradPeGrid.get(i,j,k)->at(fsgrids::egradpe::EZGRADPE); + if (Parameters::ohmGradPeTerm > 0) { + Ez_SW += EGradPeGrid.get(i, j, k)->at(fsgrids::egradpe::EZGRADPE); } - - #ifndef FS_1ST_ORDER_SPACE - // 2nd order terms: - Ez_SW += +HALF*((Bx_S - HALF*dBxdy_S)*(-dmoments_SW->at(fsgrids::dmoments::dVydx) - dmoments_SW->at(fsgrids::dmoments::dVydy)) - dBxdy_S*Vy0 + SIXTH*dBxdz_S*dmoments_SW->at(fsgrids::dmoments::dVydz)); - Ez_SW += -HALF*((By_W - HALF*dBydx_W)*(-dmoments_SW->at(fsgrids::dmoments::dVxdx) - dmoments_SW->at(fsgrids::dmoments::dVxdy)) - dBydx_W*Vx0 + SIXTH*dBydz_W*dmoments_SW->at(fsgrids::dmoments::dVxdz)); - #endif - - // Calculate maximum wave speed (fast magnetosonic speed) on SW cell. In order + +#ifndef FS_1ST_ORDER_SPACE + // 2nd order terms: + Ez_SW += +HALF * ((Bx_S - HALF * dBxdy_S) * + (-dmoments_SW->at(fsgrids::dmoments::dVydx) - dmoments_SW->at(fsgrids::dmoments::dVydy)) - + dBxdy_S * Vy0 + SIXTH * dBxdz_S * dmoments_SW->at(fsgrids::dmoments::dVydz)); + Ez_SW += -HALF * ((By_W - HALF * dBydx_W) * + (-dmoments_SW->at(fsgrids::dmoments::dVxdx) - dmoments_SW->at(fsgrids::dmoments::dVxdy)) - + dBydx_W * Vx0 + SIXTH * dBydz_W * dmoments_SW->at(fsgrids::dmoments::dVxdz)); +#endif + + // Calculate maximum wave speed (fast magnetosonic speed) on SW cell. In order // to get Alfven speed we need to calculate some reconstruction coeff. for Bz: - calculateWaveSpeedXY( - perBGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - i, j, k, - i, j, k+1, - Bx_S, By_W, dBxdy_S, dBxdz_S, dBydx_W, dBydz_W, MINUS, MINUS, minRhom, maxRhom, vA, vS, vW - ); - c_x = min(Parameters::maxWaveVelocity,sqrt(vA*vA + vS*vS) + vW); + calculateWaveSpeedXY(perBGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + i, + j, + k, + i, + j, + k + 1, + Bx_S, + By_W, + dBxdy_S, + dBxdz_S, + dBydx_W, + dBydz_W, + MINUS, + MINUS, + minRhom, + maxRhom, + vA, + vS, + vW); + c_x = min(Parameters::maxWaveVelocity, sqrt(vA * vA + vS * vS) + vW); c_y = c_x; - ax_neg = max(ZERO,-Vx0 + c_x); - ax_pos = max(ZERO,+Vx0 + c_x); - ay_neg = max(ZERO,-Vy0 + c_y); - ay_pos = max(ZERO,+Vy0 + c_y); + ax_neg = max(ZERO, -Vx0 + c_x); + ax_pos = max(ZERO, +Vx0 + c_x); + ay_neg = max(ZERO, -Vy0 + c_y); + ay_pos = max(ZERO, +Vy0 + c_y); maxV = max(maxV, calculateCflSpeed(Vx0, Vy0, vA, vS, vW)); // Ez and characteristic speeds on SE (i-1) cell: - Vx0 = moments_SE->at(fsgrids::moments::VX); - Vy0 = moments_SE->at(fsgrids::moments::VY); - + Vx0 = moments_SE->at(fsgrids::moments::VX); + Vy0 = moments_SE->at(fsgrids::moments::VY); + // 1st order terms: - Real Ez_SE = Bx_S*Vy0 - By_E*Vx0; + Real Ez_SE = Bx_S * Vy0 - By_E * Vx0; // Resistive term if (Parameters::resistivity > 0) { Ez_SE += Parameters::resistivity * - sqrt((bgb_SE->at(fsgrids::bgbfield::BGBX)+perb_SE->at(fsgrids::bfield::PERBX))* - (bgb_SE->at(fsgrids::bgbfield::BGBX)+perb_SE->at(fsgrids::bfield::PERBX)) + - (bgb_SE->at(fsgrids::bgbfield::BGBY)+perb_SE->at(fsgrids::bfield::PERBY))* - (bgb_SE->at(fsgrids::bgbfield::BGBY)+perb_SE->at(fsgrids::bfield::PERBY)) + - (bgb_SE->at(fsgrids::bgbfield::BGBZ)+perb_SE->at(fsgrids::bfield::PERBZ))* - (bgb_SE->at(fsgrids::bgbfield::BGBZ)+perb_SE->at(fsgrids::bfield::PERBZ)) - ) / - moments_SE->at(fsgrids::moments::RHOQ) / - physicalconstants::MU_0 * - (dperb_SE->at(fsgrids::dperb::dPERBydx)/technicalGrid.DX - dperb_SE->at(fsgrids::dperb::dPERBxdy)/technicalGrid.DY); + sqrt((bgb_SE->at(fsgrids::bgbfield::BGBX) + perb_SE->at(fsgrids::bfield::PERBX)) * + (bgb_SE->at(fsgrids::bgbfield::BGBX) + perb_SE->at(fsgrids::bfield::PERBX)) + + (bgb_SE->at(fsgrids::bgbfield::BGBY) + perb_SE->at(fsgrids::bfield::PERBY)) * + (bgb_SE->at(fsgrids::bgbfield::BGBY) + perb_SE->at(fsgrids::bfield::PERBY)) + + (bgb_SE->at(fsgrids::bgbfield::BGBZ) + perb_SE->at(fsgrids::bfield::PERBZ)) * + (bgb_SE->at(fsgrids::bgbfield::BGBZ) + perb_SE->at(fsgrids::bfield::PERBZ))) / + moments_SE->at(fsgrids::moments::RHOQ) / physicalconstants::MU_0 * + (dperb_SE->at(fsgrids::dperb::dPERBydx) / technicalGrid.DX - + dperb_SE->at(fsgrids::dperb::dPERBxdy) / technicalGrid.DY); } - + // Hall term if (Parameters::ohmHallTerm > 0) { - Ez_SE += EHallGrid.get(i-1,j,k)->at(fsgrids::ehall::EZHALL_100_101); + Ez_SE += EHallGrid.get(i - 1, j, k)->at(fsgrids::ehall::EZHALL_100_101); } - + // Electron pressure gradient term - if(Parameters::ohmGradPeTerm > 0) { - Ez_SE += EGradPeGrid.get(i-1,j,k)->at(fsgrids::egradpe::EZGRADPE); + if (Parameters::ohmGradPeTerm > 0) { + Ez_SE += EGradPeGrid.get(i - 1, j, k)->at(fsgrids::egradpe::EZGRADPE); } - - #ifndef FS_1ST_ORDER_SPACE - // 2nd order terms: - Ez_SE += +HALF*((Bx_S - HALF*dBxdy_S)*(+dmoments_SE->at(fsgrids::dmoments::dVydx) - dmoments_SE->at(fsgrids::dmoments::dVydy)) - dBxdy_S*Vy0 + SIXTH*dBxdz_S*dmoments_SE->at(fsgrids::dmoments::dVydz)); - Ez_SE += -HALF*((By_E + HALF*dBydx_E)*(+dmoments_SE->at(fsgrids::dmoments::dVxdx) - dmoments_SE->at(fsgrids::dmoments::dVxdy)) + dBydx_E*Vx0 + SIXTH*dBydz_E*dmoments_SE->at(fsgrids::dmoments::dVxdz)); - #endif - - calculateWaveSpeedXY( - perBGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - i-1,j ,k, - i-1,j ,k+1, - Bx_S, By_E, dBxdy_S, dBxdz_S, dBydx_E, dBydz_E, PLUS, MINUS, minRhom, maxRhom, vA, vS, vW - ); - c_x = min(Parameters::maxWaveVelocity,sqrt(vA*vA + vS*vS) + vW); + +#ifndef FS_1ST_ORDER_SPACE + // 2nd order terms: + Ez_SE += +HALF * ((Bx_S - HALF * dBxdy_S) * + (+dmoments_SE->at(fsgrids::dmoments::dVydx) - dmoments_SE->at(fsgrids::dmoments::dVydy)) - + dBxdy_S * Vy0 + SIXTH * dBxdz_S * dmoments_SE->at(fsgrids::dmoments::dVydz)); + Ez_SE += -HALF * ((By_E + HALF * dBydx_E) * + (+dmoments_SE->at(fsgrids::dmoments::dVxdx) - dmoments_SE->at(fsgrids::dmoments::dVxdy)) + + dBydx_E * Vx0 + SIXTH * dBydz_E * dmoments_SE->at(fsgrids::dmoments::dVxdz)); +#endif + + calculateWaveSpeedXY(perBGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + i - 1, + j, + k, + i - 1, + j, + k + 1, + Bx_S, + By_E, + dBxdy_S, + dBxdz_S, + dBydx_E, + dBydz_E, + PLUS, + MINUS, + minRhom, + maxRhom, + vA, + vS, + vW); + c_x = min(Parameters::maxWaveVelocity, sqrt(vA * vA + vS * vS) + vW); c_y = c_x; - ax_neg = max(ax_neg,-Vx0 + c_x); - ax_pos = max(ax_pos,+Vx0 + c_x); - ay_neg = max(ay_neg,-Vy0 + c_y); - ay_pos = max(ay_pos,+Vy0 + c_y); + ax_neg = max(ax_neg, -Vx0 + c_x); + ax_pos = max(ax_pos, +Vx0 + c_x); + ay_neg = max(ay_neg, -Vy0 + c_y); + ay_pos = max(ay_pos, +Vy0 + c_y); maxV = max(maxV, calculateCflSpeed(Vx0, Vy0, vA, vS, vW)); // Ez and characteristic speeds on NW (j-1) cell: - Vx0 = moments_NW->at(fsgrids::moments::VX); - Vy0 = moments_NW->at(fsgrids::moments::VY); - + Vx0 = moments_NW->at(fsgrids::moments::VX); + Vy0 = moments_NW->at(fsgrids::moments::VY); + // 1st order terms: - Real Ez_NW = Bx_N*Vy0 - By_W*Vx0; - + Real Ez_NW = Bx_N * Vy0 - By_W * Vx0; + // Resistive term if (Parameters::resistivity > 0) { Ez_NW += Parameters::resistivity * - sqrt((bgb_NW->at(fsgrids::bgbfield::BGBX)+perb_NW->at(fsgrids::bfield::PERBX))* - (bgb_NW->at(fsgrids::bgbfield::BGBX)+perb_NW->at(fsgrids::bfield::PERBX)) + - (bgb_NW->at(fsgrids::bgbfield::BGBY)+perb_NW->at(fsgrids::bfield::PERBY))* - (bgb_NW->at(fsgrids::bgbfield::BGBY)+perb_NW->at(fsgrids::bfield::PERBY)) + - (bgb_NW->at(fsgrids::bgbfield::BGBZ)+perb_NW->at(fsgrids::bfield::PERBZ))* - (bgb_NW->at(fsgrids::bgbfield::BGBZ)+perb_NW->at(fsgrids::bfield::PERBZ)) - ) / - moments_NW->at(fsgrids::moments::RHOQ) / - physicalconstants::MU_0 * - (dperb_NW->at(fsgrids::dperb::dPERBydx)/technicalGrid.DX - dperb_NW->at(fsgrids::dperb::dPERBxdy)/technicalGrid.DY); + sqrt((bgb_NW->at(fsgrids::bgbfield::BGBX) + perb_NW->at(fsgrids::bfield::PERBX)) * + (bgb_NW->at(fsgrids::bgbfield::BGBX) + perb_NW->at(fsgrids::bfield::PERBX)) + + (bgb_NW->at(fsgrids::bgbfield::BGBY) + perb_NW->at(fsgrids::bfield::PERBY)) * + (bgb_NW->at(fsgrids::bgbfield::BGBY) + perb_NW->at(fsgrids::bfield::PERBY)) + + (bgb_NW->at(fsgrids::bgbfield::BGBZ) + perb_NW->at(fsgrids::bfield::PERBZ)) * + (bgb_NW->at(fsgrids::bgbfield::BGBZ) + perb_NW->at(fsgrids::bfield::PERBZ))) / + moments_NW->at(fsgrids::moments::RHOQ) / physicalconstants::MU_0 * + (dperb_NW->at(fsgrids::dperb::dPERBydx) / technicalGrid.DX - + dperb_NW->at(fsgrids::dperb::dPERBxdy) / technicalGrid.DY); } - + // Hall term - if(Parameters::ohmHallTerm > 0) { - Ez_NW += EHallGrid.get(i,j-1,k)->at(fsgrids::ehall::EZHALL_010_011); + if (Parameters::ohmHallTerm > 0) { + Ez_NW += EHallGrid.get(i, j - 1, k)->at(fsgrids::ehall::EZHALL_010_011); } - + // Electron pressure gradient term - if(Parameters::ohmGradPeTerm > 0) { - Ez_NW += EGradPeGrid.get(i,j-1,k)->at(fsgrids::egradpe::EZGRADPE); + if (Parameters::ohmGradPeTerm > 0) { + Ez_NW += EGradPeGrid.get(i, j - 1, k)->at(fsgrids::egradpe::EZGRADPE); } - - #ifndef FS_1ST_ORDER_SPACE - // 2nd order terms: - Ez_NW += +HALF*((Bx_N + HALF*dBxdy_N)*(-dmoments_NW->at(fsgrids::dmoments::dVydx) + dmoments_NW->at(fsgrids::dmoments::dVydy)) + dBxdy_N*Vy0 + SIXTH*dBxdz_N*dmoments_NW->at(fsgrids::dmoments::dVydz)); - Ez_NW += -HALF*((By_W - HALF*dBydx_W)*(-dmoments_NW->at(fsgrids::dmoments::dVxdx) + dmoments_NW->at(fsgrids::dmoments::dVxdy)) - dBydx_W*Vx0 + SIXTH*dBydz_W*dmoments_NW->at(fsgrids::dmoments::dVxdz)); - #endif - - calculateWaveSpeedXY( - perBGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - i, j-1, k, - i, j-1, k+1, - Bx_N, By_W, dBxdy_N, dBxdz_N, dBydx_W, dBydz_W, MINUS, PLUS, minRhom, maxRhom, vA, vS, vW - ); - c_x = min(Parameters::maxWaveVelocity,sqrt(vA*vA + vS*vS) + vW); + +#ifndef FS_1ST_ORDER_SPACE + // 2nd order terms: + Ez_NW += +HALF * ((Bx_N + HALF * dBxdy_N) * + (-dmoments_NW->at(fsgrids::dmoments::dVydx) + dmoments_NW->at(fsgrids::dmoments::dVydy)) + + dBxdy_N * Vy0 + SIXTH * dBxdz_N * dmoments_NW->at(fsgrids::dmoments::dVydz)); + Ez_NW += -HALF * ((By_W - HALF * dBydx_W) * + (-dmoments_NW->at(fsgrids::dmoments::dVxdx) + dmoments_NW->at(fsgrids::dmoments::dVxdy)) - + dBydx_W * Vx0 + SIXTH * dBydz_W * dmoments_NW->at(fsgrids::dmoments::dVxdz)); +#endif + + calculateWaveSpeedXY(perBGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + i, + j - 1, + k, + i, + j - 1, + k + 1, + Bx_N, + By_W, + dBxdy_N, + dBxdz_N, + dBydx_W, + dBydz_W, + MINUS, + PLUS, + minRhom, + maxRhom, + vA, + vS, + vW); + c_x = min(Parameters::maxWaveVelocity, sqrt(vA * vA + vS * vS) + vW); c_y = c_x; - ax_neg = max(ax_neg,-Vx0 + c_x); - ax_pos = max(ax_pos,+Vx0 + c_x); - ay_neg = max(ay_neg,-Vy0 + c_y); - ay_pos = max(ay_pos,+Vy0 + c_y); + ax_neg = max(ax_neg, -Vx0 + c_x); + ax_pos = max(ax_pos, +Vx0 + c_x); + ay_neg = max(ay_neg, -Vy0 + c_y); + ay_pos = max(ay_pos, +Vy0 + c_y); maxV = max(maxV, calculateCflSpeed(Vx0, Vy0, vA, vS, vW)); - + // Ez and characteristic speeds on NE (i-1,j-1) cell: - Vx0 = moments_NE->at(fsgrids::moments::VX); - Vy0 = moments_NE->at(fsgrids::moments::VY); - + Vx0 = moments_NE->at(fsgrids::moments::VX); + Vy0 = moments_NE->at(fsgrids::moments::VY); + // 1st order terms: - Real Ez_NE = Bx_N*Vy0 - By_E*Vx0; - + Real Ez_NE = Bx_N * Vy0 - By_E * Vx0; + // Resistive term if (Parameters::resistivity > 0) { Ez_NE += Parameters::resistivity * - sqrt((bgb_NE->at(fsgrids::bgbfield::BGBX)+perb_NE->at(fsgrids::bfield::PERBX))* - (bgb_NE->at(fsgrids::bgbfield::BGBX)+perb_NE->at(fsgrids::bfield::PERBX)) + - (bgb_NE->at(fsgrids::bgbfield::BGBY)+perb_NE->at(fsgrids::bfield::PERBY))* - (bgb_NE->at(fsgrids::bgbfield::BGBY)+perb_NE->at(fsgrids::bfield::PERBY)) + - (bgb_NE->at(fsgrids::bgbfield::BGBZ)+perb_NE->at(fsgrids::bfield::PERBZ))* - (bgb_NE->at(fsgrids::bgbfield::BGBZ)+perb_NE->at(fsgrids::bfield::PERBZ)) - ) / - moments_NE->at(fsgrids::moments::RHOQ) / - physicalconstants::MU_0 * - (dperb_NE->at(fsgrids::dperb::dPERBydx)/technicalGrid.DX - dperb_NE->at(fsgrids::dperb::dPERBxdy)/technicalGrid.DY); + sqrt((bgb_NE->at(fsgrids::bgbfield::BGBX) + perb_NE->at(fsgrids::bfield::PERBX)) * + (bgb_NE->at(fsgrids::bgbfield::BGBX) + perb_NE->at(fsgrids::bfield::PERBX)) + + (bgb_NE->at(fsgrids::bgbfield::BGBY) + perb_NE->at(fsgrids::bfield::PERBY)) * + (bgb_NE->at(fsgrids::bgbfield::BGBY) + perb_NE->at(fsgrids::bfield::PERBY)) + + (bgb_NE->at(fsgrids::bgbfield::BGBZ) + perb_NE->at(fsgrids::bfield::PERBZ)) * + (bgb_NE->at(fsgrids::bgbfield::BGBZ) + perb_NE->at(fsgrids::bfield::PERBZ))) / + moments_NE->at(fsgrids::moments::RHOQ) / physicalconstants::MU_0 * + (dperb_NE->at(fsgrids::dperb::dPERBydx) / technicalGrid.DX - + dperb_NE->at(fsgrids::dperb::dPERBxdy) / technicalGrid.DY); } - + // Hall term - if(Parameters::ohmHallTerm > 0) { - Ez_NE += EHallGrid.get(i-1,j-1,k)->at(fsgrids::ehall::EZHALL_110_111); + if (Parameters::ohmHallTerm > 0) { + Ez_NE += EHallGrid.get(i - 1, j - 1, k)->at(fsgrids::ehall::EZHALL_110_111); } - + // Electron pressure gradient term - if(Parameters::ohmGradPeTerm > 0) { - Ez_NE += EGradPeGrid.get(i-1,j-1,k)->at(fsgrids::egradpe::EZGRADPE); + if (Parameters::ohmGradPeTerm > 0) { + Ez_NE += EGradPeGrid.get(i - 1, j - 1, k)->at(fsgrids::egradpe::EZGRADPE); } - - #ifndef FS_1ST_ORDER_SPACE - // 2nd order terms: - Ez_NE += +HALF*((Bx_N + HALF*dBxdy_N)*(+dmoments_NE->at(fsgrids::dmoments::dVydx) + dmoments_NE->at(fsgrids::dmoments::dVydy)) + dBxdy_N*Vy0 + SIXTH*dBxdz_N*dmoments_NE->at(fsgrids::dmoments::dVydz)); - Ez_NE += -HALF*((By_E + HALF*dBydx_E)*(+dmoments_NE->at(fsgrids::dmoments::dVxdx) + dmoments_NE->at(fsgrids::dmoments::dVxdy)) + dBydx_E*Vx0 + SIXTH*dBydz_E*dmoments_NE->at(fsgrids::dmoments::dVxdz)); - #endif - - calculateWaveSpeedXY( - perBGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - i-1,j-1,k, - i-1,j-1,k+1, - Bx_N, By_E, dBxdy_N, dBxdz_N, dBydx_E, dBydz_E, PLUS, PLUS, minRhom, maxRhom, vA, vS, vW - ); - c_x = min(Parameters::maxWaveVelocity,sqrt(vA*vA + vS*vS) + vW); + +#ifndef FS_1ST_ORDER_SPACE + // 2nd order terms: + Ez_NE += +HALF * ((Bx_N + HALF * dBxdy_N) * + (+dmoments_NE->at(fsgrids::dmoments::dVydx) + dmoments_NE->at(fsgrids::dmoments::dVydy)) + + dBxdy_N * Vy0 + SIXTH * dBxdz_N * dmoments_NE->at(fsgrids::dmoments::dVydz)); + Ez_NE += -HALF * ((By_E + HALF * dBydx_E) * + (+dmoments_NE->at(fsgrids::dmoments::dVxdx) + dmoments_NE->at(fsgrids::dmoments::dVxdy)) + + dBydx_E * Vx0 + SIXTH * dBydz_E * dmoments_NE->at(fsgrids::dmoments::dVxdz)); +#endif + + calculateWaveSpeedXY(perBGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + i - 1, + j - 1, + k, + i - 1, + j - 1, + k + 1, + Bx_N, + By_E, + dBxdy_N, + dBxdz_N, + dBydx_E, + dBydz_E, + PLUS, + PLUS, + minRhom, + maxRhom, + vA, + vS, + vW); + c_x = min(Parameters::maxWaveVelocity, sqrt(vA * vA + vS * vS) + vW); c_y = c_x; - ax_neg = max(ax_neg,-Vx0 + c_x); - ax_pos = max(ax_pos,+Vx0 + c_x); - ay_neg = max(ay_neg,-Vy0 + c_y); - ay_pos = max(ay_pos,+Vy0 + c_y); + ax_neg = max(ax_neg, -Vx0 + c_x); + ax_pos = max(ax_pos, +Vx0 + c_x); + ay_neg = max(ay_neg, -Vy0 + c_y); + ay_pos = max(ay_pos, +Vy0 + c_y); maxV = max(maxV, calculateCflSpeed(Vx0, Vy0, vA, vS, vW)); // Calculate properly upwinded edge-averaged Ez: - efield_SW->at(fsgrids::efield::EZ) = ax_pos*ay_pos*Ez_NE + ax_pos*ay_neg*Ez_SE + ax_neg*ay_pos*Ez_NW + ax_neg*ay_neg*Ez_SW; - efield_SW->at(fsgrids::efield::EZ) /= ((ax_pos+ax_neg)*(ay_pos+ay_neg)+EPS); + efield_SW->at(fsgrids::efield::EZ) = + ax_pos * ay_pos * Ez_NE + ax_pos * ay_neg * Ez_SE + ax_neg * ay_pos * Ez_NW + ax_neg * ay_neg * Ez_SW; + efield_SW->at(fsgrids::efield::EZ) /= ((ax_pos + ax_neg) * (ay_pos + ay_neg) + EPS); if (Parameters::fieldSolverDiffusiveEterms) { #ifdef FS_1ST_ORDER_SPACE - efield_SW->at(fsgrids::efield::EZ) -= ay_pos*ay_neg/(ay_pos+ay_neg+EPS)*(perBx_S-perBx_N); - efield_SW->at(fsgrids::efield::EZ) += ax_pos*ax_neg/(ax_pos+ax_neg+EPS)*(perBy_W-perBy_E); + efield_SW->at(fsgrids::efield::EZ) -= ay_pos * ay_neg / (ay_pos + ay_neg + EPS) * (perBx_S - perBx_N); + efield_SW->at(fsgrids::efield::EZ) += ax_pos * ax_neg / (ax_pos + ax_neg + EPS) * (perBy_W - perBy_E); #else - efield_SW->at(fsgrids::efield::EZ) -= ay_pos*ay_neg/(ay_pos+ay_neg+EPS)*((perBx_S-HALF*dperBxdy_S) - (perBx_N+HALF*dperBxdy_N)); - efield_SW->at(fsgrids::efield::EZ) += ax_pos*ax_neg/(ax_pos+ax_neg+EPS)*((perBy_W-HALF*dperBydx_W) - (perBy_E+HALF*dperBydx_E)); + efield_SW->at(fsgrids::efield::EZ) -= + ay_pos * ay_neg / (ay_pos + ay_neg + EPS) * ((perBx_S - HALF * dperBxdy_S) - (perBx_N + HALF * dperBxdy_N)); + efield_SW->at(fsgrids::efield::EZ) += + ax_pos * ax_neg / (ax_pos + ax_neg + EPS) * ((perBy_W - HALF * dperBydx_W) - (perBy_E + HALF * dperBydx_E)); #endif } - + if ((RKCase == RK_ORDER1) || (RKCase == RK_ORDER2_STEP2)) { - //compute maximum timestep for fieldsolver in this cell (CFL=1) - Real min_dx=std::numeric_limits::max();; - min_dx=min(min_dx,technicalGrid.DX); - min_dx=min(min_dx,technicalGrid.DY); - //update max allowed timestep for field propagation in this cell, which is the minimum of CFL=1 timesteps - if(maxV!=ZERO) technicalGrid.get(i,j,k)->maxFsDt=min(technicalGrid.get(i,j,k)->maxFsDt,min_dx/maxV); + // compute maximum timestep for fieldsolver in this cell (CFL=1) + Real min_dx = std::numeric_limits::max(); + ; + min_dx = min(min_dx, technicalGrid.DX); + min_dx = min(min_dx, technicalGrid.DY); + // update max allowed timestep for field propagation in this cell, which is the minimum of CFL=1 timesteps + if (maxV != ZERO) + technicalGrid.get(i, j, k)->maxFsDt = min(technicalGrid.get(i, j, k)->maxFsDt, min_dx / maxV); } } /*! \brief Electric field propagation function. - * + * * Calls the general or the system boundary electric field propagation functions. - * + * * \param perBGrid fsGrid holding the perturbed B quantities * \param EGrid fsGrid holding the electric field * \param EHallGrid fsGrid holding the Hall contributions to the electric field @@ -1534,97 +1787,91 @@ void calculateEdgeElectricFieldZ( * \param i,j,k fsGrid cell coordinates for the current cell * \param sysBoundaries System boundary conditions existing * \param RKCase Element in the enum defining the Runge-Kutta method steps - * - * \sa calculateUpwindedElectricFieldSimple calculateEdgeElectricFieldX calculateEdgeElectricFieldY calculateEdgeElectricFieldZ - * + * + * \sa calculateUpwindedElectricFieldSimple calculateEdgeElectricFieldX calculateEdgeElectricFieldY + * calculateEdgeElectricFieldZ + * */ -void calculateElectricField( - FsGrid< std::array, FS_STENCIL_WIDTH> & perBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EHallGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, - FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, - cint i, - cint j, - cint k, - SysBoundary& sysBoundaries, - cint& RKCase -) { - cuint cellSysBoundaryFlag = technicalGrid.get(i,j,k)->sysBoundaryFlag; - - if (cellSysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) return; - - cuint bitfield = technicalGrid.get(i,j,k)->SOLVE; - +void calculateElectricField(FsGrid, FS_STENCIL_WIDTH>& perBGrid, + FsGrid, FS_STENCIL_WIDTH>& EGrid, + FsGrid, FS_STENCIL_WIDTH>& EHallGrid, + FsGrid, FS_STENCIL_WIDTH>& EGradPeGrid, + FsGrid, FS_STENCIL_WIDTH>& momentsGrid, + FsGrid, FS_STENCIL_WIDTH>& dPerBGrid, + FsGrid, FS_STENCIL_WIDTH>& dMomentsGrid, + FsGrid, FS_STENCIL_WIDTH>& BgBGrid, + FsGrid& technicalGrid, + cint i, + cint j, + cint k, + SysBoundary& sysBoundaries, + cint& RKCase) { + cuint cellSysBoundaryFlag = technicalGrid.get(i, j, k)->sysBoundaryFlag; + + if (cellSysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE) + return; + + cuint bitfield = technicalGrid.get(i, j, k)->SOLVE; + if ((bitfield & compute::EX) == compute::EX) { - calculateEdgeElectricFieldX( - perBGrid, - EGrid, - EHallGrid, - EGradPeGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - technicalGrid, - i, - j, - k, - RKCase - ); + calculateEdgeElectricFieldX(perBGrid, + EGrid, + EHallGrid, + EGradPeGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + technicalGrid, + i, + j, + k, + RKCase); } else { sysBoundaries.getSysBoundary(cellSysBoundaryFlag)->fieldSolverBoundaryCondElectricField(EGrid, i, j, k, 0); } - + if ((bitfield & compute::EY) == compute::EY) { - calculateEdgeElectricFieldY( - perBGrid, - EGrid, - EHallGrid, - EGradPeGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - technicalGrid, - i, - j, - k, - RKCase - ); + calculateEdgeElectricFieldY(perBGrid, + EGrid, + EHallGrid, + EGradPeGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + technicalGrid, + i, + j, + k, + RKCase); } else { sysBoundaries.getSysBoundary(cellSysBoundaryFlag)->fieldSolverBoundaryCondElectricField(EGrid, i, j, k, 1); } - + if ((bitfield & compute::EZ) == compute::EZ) { - calculateEdgeElectricFieldZ( - perBGrid, - EGrid, - EHallGrid, - EGradPeGrid, - momentsGrid, - dPerBGrid, - dMomentsGrid, - BgBGrid, - technicalGrid, - i, - j, - k, - RKCase - ); + calculateEdgeElectricFieldZ(perBGrid, + EGrid, + EHallGrid, + EGradPeGrid, + momentsGrid, + dPerBGrid, + dMomentsGrid, + BgBGrid, + technicalGrid, + i, + j, + k, + RKCase); } else { sysBoundaries.getSysBoundary(cellSysBoundaryFlag)->fieldSolverBoundaryCondElectricField(EGrid, i, j, k, 2); } } /*! \brief High-level electric field computation function. - * + * * Transfers the derivatives, calculates the edge electric fields and transfers the new electric fields. - * + * * \param perBGrid fsGrid holding the perturbed B quantities at runge-kutta t=0 * \param perBDt2Grid fsGrid holding the perturbed B quantities at runge-kutta t=0.5 * \param EGrid fsGrid holding the Electric field quantities at runge-kutta t=0 @@ -1639,105 +1886,100 @@ void calculateElectricField( * \param technicalGrid fsGrid holding technical information (such as boundary types) * \param sysBoundaries System boundary conditions existing * \param RKCase Element in the enum defining the Runge-Kutta method steps - * + * * \sa calculateElectricField calculateEdgeElectricFieldX calculateEdgeElectricFieldY calculateEdgeElectricFieldZ */ void calculateUpwindedElectricFieldSimple( - FsGrid< std::array, FS_STENCIL_WIDTH> & perBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & perBDt2Grid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EDt2Grid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EHallGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, - FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, - SysBoundary& sysBoundaries, - cint& RKCase -) { - //const std::array gridDims = technicalGrid.getLocalSize(); + FsGrid, FS_STENCIL_WIDTH>& perBGrid, + FsGrid, FS_STENCIL_WIDTH>& perBDt2Grid, + FsGrid, FS_STENCIL_WIDTH>& EGrid, + FsGrid, FS_STENCIL_WIDTH>& EDt2Grid, + FsGrid, FS_STENCIL_WIDTH>& EHallGrid, + FsGrid, FS_STENCIL_WIDTH>& EGradPeGrid, + FsGrid, FS_STENCIL_WIDTH>& momentsGrid, + FsGrid, FS_STENCIL_WIDTH>& momentsDt2Grid, + FsGrid, FS_STENCIL_WIDTH>& dPerBGrid, + FsGrid, FS_STENCIL_WIDTH>& dMomentsGrid, + FsGrid, FS_STENCIL_WIDTH>& BgBGrid, + FsGrid& technicalGrid, + SysBoundary& sysBoundaries, + cint& RKCase) { + // const std::array gridDims = technicalGrid.getLocalSize(); const int* gridDims = &technicalGrid.getLocalSize()[0]; - const size_t N_cells = gridDims[0]*gridDims[1]*gridDims[2]; - phiprof::Timer upwindedETimer {"Calculate upwinded electric field"}; - int computeTimerID {phiprof::initializeTimer("Electric field compute cells")}; - - phiprof::Timer mpiTimer {"Electric field ghost updates MPI", {"MPI"}}; + const size_t N_cells = gridDims[0] * gridDims[1] * gridDims[2]; + phiprof::Timer upwindedETimer{"Calculate upwinded electric field"}; + int computeTimerID{phiprof::initializeTimer("Electric field compute cells")}; + + phiprof::Timer mpiTimer{"Electric field ghost updates MPI", {"MPI"}}; // Update ghosts if necessary, unless previous terms have already updated them - if(P::ohmHallTerm > 0) { + if (P::ohmHallTerm > 0) { EHallGrid.updateGhostCells(); } - if(P::ohmGradPeTerm > 0) { + if (P::ohmGradPeTerm > 0) { EGradPeGrid.updateGhostCells(); } - if(P::ohmHallTerm == 0) { + if (P::ohmHallTerm == 0) { dPerBGrid.updateGhostCells(); } - if(P::ohmHallTerm == 0 && P::ohmGradPeTerm == 0) { + if (P::ohmHallTerm == 0 && P::ohmGradPeTerm == 0) { dMomentsGrid.updateGhostCells(); } - + mpiTimer.stop(); - - // Calculate upwinded electric field on inner cells - #pragma omp parallel + +// Calculate upwinded electric field on inner cells +#pragma omp parallel { - phiprof::Timer computeTimer {computeTimerID}; - #pragma omp for collapse(2) - for (int k=0; k, FS_STENCIL_WIDTH> & perBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & perBDt2Grid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EDt2Grid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EHallGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & EGradPeGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & momentsDt2Grid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dPerBGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & dMomentsGrid, - FsGrid< std::array, FS_STENCIL_WIDTH> & BgBGrid, - FsGrid< fsgrids::technical, FS_STENCIL_WIDTH> & technicalGrid, - SysBoundary& sysBoundaries, - cint& RKCase -); + FsGrid, FS_STENCIL_WIDTH>& perBGrid, + FsGrid, FS_STENCIL_WIDTH>& perBDt2Grid, + FsGrid, FS_STENCIL_WIDTH>& EGrid, + FsGrid, FS_STENCIL_WIDTH>& EDt2Grid, + FsGrid, FS_STENCIL_WIDTH>& EHallGrid, + FsGrid, FS_STENCIL_WIDTH>& EGradPeGrid, + FsGrid, FS_STENCIL_WIDTH>& momentsGrid, + FsGrid, FS_STENCIL_WIDTH>& momentsDt2Grid, + FsGrid, FS_STENCIL_WIDTH>& dPerBGrid, + FsGrid, FS_STENCIL_WIDTH>& dMomentsGrid, + FsGrid, FS_STENCIL_WIDTH>& BgBGrid, + FsGrid& technicalGrid, + SysBoundary& sysBoundaries, + cint& RKCase); #endif