Add hyperresistivity term to electric field in fieldsolver

Makefile

@@ -225,7 +225,7 @@ else
 endif
 
 # Add field solver objects
-OBJS_FSOLVER = 	ldz_magnetic_field.o ldz_volume.o derivatives.o ldz_electric_field.o ldz_hall.o ldz_gradpe.o
+OBJS_FSOLVER = 	ldz_magnetic_field.o ldz_volume.o derivatives.o ldz_electric_field.o ldz_hall.o ldz_gradpe.o ldz_hyper.o
 
 # Include autogenerated dependency files, if they exist
 -include $(OBJS_FSOLVER:%.o=%.d)

common.h

@@ -321,6 +321,13 @@ namespace fsgrids {
       N_EGRADPE
    };
 
+   enum ehyper {
+      EXHYPER,         /*!< Hyperresistivity term x.*/
+      EYHYPER,         /*!< Hyperresistivity term y.*/
+      EZHYPER,         /*!< Hyperresistivity term z.*/
+      N_EHYPER
+   };
+
    enum moments {
       RHOM, /*!< Overall mass density. Calculated by Vlasov propagator, used to propagate fields.*/
       RHOQ, /*!< Overall charge density. Calculated by Vlasov propagator, used to propagate fields.*/

datareduction/datareducer.h

@@ -60,6 +60,7 @@ class DataReducer {
                       FsGrid< std::array<Real, fsgrids::efield::N_EFIELD>, FS_STENCIL_WIDTH> & EGrid,
                       FsGrid< std::array<Real, fsgrids::ehall::N_EHALL>, FS_STENCIL_WIDTH> & EHallGrid,
                       FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeGrid,
+                      FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperGrid,
                       FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsGrid,
                       FsGrid< std::array<Real, fsgrids::dperb::N_DPERB>, FS_STENCIL_WIDTH> & dPerBGrid,
                       FsGrid< std::array<Real, fsgrids::dmoments::N_DMOMENTS>, FS_STENCIL_WIDTH> & dMomentsGrid,

datareduction/datareductionoperator.cpp

@@ -132,6 +132,7 @@ namespace DRO {
                       FsGrid< std::array<Real, fsgrids::efield::N_EFIELD>, FS_STENCIL_WIDTH> & EGrid,
                       FsGrid< std::array<Real, fsgrids::ehall::N_EHALL>, FS_STENCIL_WIDTH> & EHallGrid,
                       FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeGrid,
+                      FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperGrid,
                       FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsGrid,
                       FsGrid< std::array<Real, fsgrids::dperb::N_DPERB>, FS_STENCIL_WIDTH> & dPerBGrid,
                       FsGrid< std::array<Real, fsgrids::dmoments::N_DMOMENTS>, FS_STENCIL_WIDTH> & dMomentsGrid,
@@ -150,7 +151,7 @@ namespace DRO {
       attribs["variableLaTeX"]=variableLaTeX;
 
       std::vector<double> varBuffer =
-         lambda(perBGrid,EGrid,EHallGrid,EGradPeGrid,momentsGrid,dPerBGrid,dMomentsGrid,BgBGrid,volGrid,technicalGrid);
+         lambda(perBGrid,EGrid,EHallGrid,EGradPeGrid,EHyperGrid,momentsGrid,dPerBGrid,dMomentsGrid,BgBGrid,volGrid,technicalGrid);
 
       std::array<FsGridTools::FsIndex_t,3>& gridSize = technicalGrid.getLocalSize();
       int vectorSize;

datareduction/datareductionoperator.h

@@ -98,6 +98,7 @@ namespace DRO {
                       FsGrid< std::array<Real, fsgrids::efield::N_EFIELD>, FS_STENCIL_WIDTH> & EGrid,
                       FsGrid< std::array<Real, fsgrids::ehall::N_EHALL>, FS_STENCIL_WIDTH> & EHallGrid,
                       FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeGrid,
+                      FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperGrid,
                       FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsGrid,
                       FsGrid< std::array<Real, fsgrids::dperb::N_DPERB>, FS_STENCIL_WIDTH> & dPerBGrid,
                       FsGrid< std::array<Real, fsgrids::dmoments::N_DMOMENTS>, FS_STENCIL_WIDTH> & dMomentsGrid,
@@ -120,6 +121,7 @@ namespace DRO {
                       FsGrid< std::array<Real, fsgrids::efield::N_EFIELD>, FS_STENCIL_WIDTH> & EGrid,
                       FsGrid< std::array<Real, fsgrids::ehall::N_EHALL>, FS_STENCIL_WIDTH> & EHallGrid,
                       FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeGrid,
+                      FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperGrid,
                       FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsGrid,
                       FsGrid< std::array<Real, fsgrids::dperb::N_DPERB>, FS_STENCIL_WIDTH> & dPerBGrid,
                       FsGrid< std::array<Real, fsgrids::dmoments::N_DMOMENTS>, FS_STENCIL_WIDTH> & dMomentsGrid,

fieldsolver/fs_common.h

@@ -70,6 +70,8 @@ bool propagateFields(
    FsGrid< std::array<Real, fsgrids::ehall::N_EHALL>, FS_STENCIL_WIDTH> & EHallGrid,
    FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeGrid,
    FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeDt2Grid,
+   FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperGrid,
+   FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperDt2Grid,
    FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsGrid,
    FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsDt2Grid,
    FsGrid< std::array<Real, fsgrids::dperb::N_DPERB>, FS_STENCIL_WIDTH> & dPerBGrid,

fieldsolver/ldz_electric_field.cpp

@@ -456,6 +456,7 @@ void calculateEdgeElectricFieldX(
    FsGrid< std::array<Real, fsgrids::efield::N_EFIELD>, FS_STENCIL_WIDTH> & EGrid,
    FsGrid< std::array<Real, fsgrids::ehall::N_EHALL>, FS_STENCIL_WIDTH> & EHallGrid,
    FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeGrid,
+   FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperGrid,
    FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsGrid,
    FsGrid< std::array<Real, fsgrids::dperb::N_DPERB>, FS_STENCIL_WIDTH> & dPerBGrid,
    FsGrid< std::array<Real, fsgrids::dmoments::N_DMOMENTS>, FS_STENCIL_WIDTH> & dMomentsGrid,
@@ -584,6 +585,11 @@ void calculateEdgeElectricFieldX(
       Ex_SW += EGradPeGrid.get(i,j,k)->at(fsgrids::egradpe::EXGRADPE);
    }
 
+   // Hyperresistivity term
+   if(Parameters::ohmHyperTerm > 0) {
+      Ex_SW += EHyperGrid.get(i,j,k)->at(fsgrids::ehyper::EXHYPER);
+   }
+
    #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));
@@ -638,6 +644,11 @@ void calculateEdgeElectricFieldX(
    if(Parameters::ohmGradPeTerm > 0) {
       Ex_SE += EGradPeGrid.get(i,j-1,k)->at(fsgrids::egradpe::EXGRADPE);
    }
+
+   // Hyperresistivity term
+   if(Parameters::ohmHyperTerm > 0) {
+      Ex_SE += EHyperGrid.get(i,j-1,k)->at(fsgrids::ehyper::EXHYPER);
+   }
 
    #ifndef FS_1ST_ORDER_SPACE
       // 2nd order terms:
@@ -694,6 +705,11 @@ void calculateEdgeElectricFieldX(
    if(Parameters::ohmGradPeTerm > 0) {
       Ex_NW += EGradPeGrid.get(i,j,k-1)->at(fsgrids::egradpe::EXGRADPE);
    }
+
+   // Hyperresistivity term
+   if(Parameters::ohmHyperTerm > 0) {
+      Ex_NW += EHyperGrid.get(i,j,k-1)->at(fsgrids::ehyper::EXHYPER);
+   }
 
    #ifndef FS_1ST_ORDER_SPACE
       // 2nd order terms:
@@ -750,6 +766,11 @@ void calculateEdgeElectricFieldX(
    if(Parameters::ohmGradPeTerm > 0) {
       Ex_NE += EGradPeGrid.get(i,j-1,k-1)->at(fsgrids::egradpe::EXGRADPE);
    }
+
+   // Hyperresistivity term
+   if(Parameters::ohmHyperTerm > 0) {
+      Ex_NE += EHyperGrid.get(i,j-1,k-1)->at(fsgrids::ehyper::EXHYPER);
+   }
 
    #ifndef FS_1ST_ORDER_SPACE
       // 2nd order terms:
@@ -814,6 +835,7 @@ void calculateEdgeElectricFieldY(
    FsGrid< std::array<Real, fsgrids::efield::N_EFIELD>, FS_STENCIL_WIDTH> & EGrid,
    FsGrid< std::array<Real, fsgrids::ehall::N_EHALL>, FS_STENCIL_WIDTH> & EHallGrid,
    FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeGrid,
+   FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperGrid,
    FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsGrid,
    FsGrid< std::array<Real, fsgrids::dperb::N_DPERB>, FS_STENCIL_WIDTH> & dPerBGrid,
    FsGrid< std::array<Real, fsgrids::dmoments::N_DMOMENTS>, FS_STENCIL_WIDTH> & dMomentsGrid,
@@ -940,6 +962,11 @@ void calculateEdgeElectricFieldY(
    if(Parameters::ohmGradPeTerm > 0) {
       Ey_SW += EGradPeGrid.get(i,j,k)->at(fsgrids::egradpe::EYGRADPE);
    }
+
+   // Hyperresistivity term
+   if(Parameters::ohmHyperTerm > 0) {
+      Ey_SW += EHyperGrid.get(i,j,k)->at(fsgrids::ehyper::EYHYPER);
+   }
 
    #ifndef FS_1ST_ORDER_SPACE
       // 2nd order terms
@@ -996,6 +1023,11 @@ void calculateEdgeElectricFieldY(
    if(Parameters::ohmGradPeTerm > 0) {
       Ey_SE += EGradPeGrid.get(i,j,k-1)->at(fsgrids::egradpe::EYGRADPE);
    }
+
+   // Hyperresistivity term
+   if(Parameters::ohmHyperTerm > 0) {
+      Ey_SE += EHyperGrid.get(i,j,k-1)->at(fsgrids::ehyper::EYHYPER);
+   }
 
    #ifndef FS_1ST_ORDER_SPACE
       // 2nd order terms:
@@ -1052,6 +1084,11 @@ void calculateEdgeElectricFieldY(
    if(Parameters::ohmGradPeTerm > 0) {
       Ey_NW += EGradPeGrid.get(i-1,j,k)->at(fsgrids::egradpe::EYGRADPE);
    }
+
+   // Hyperresistivity term
+   if(Parameters::ohmHyperTerm > 0) {
+      Ey_NW += EHyperGrid.get(i-1,j,k)->at(fsgrids::ehyper::EYHYPER);
+   }
 
    #ifndef FS_1ST_ORDER_SPACE
       // 2nd order terms:
@@ -1108,6 +1145,11 @@ void calculateEdgeElectricFieldY(
    if(Parameters::ohmGradPeTerm > 0) {
       Ey_NE += EGradPeGrid.get(i-1,j,k-1)->at(fsgrids::egradpe::EYGRADPE);
    }
+
+   // Hyperresistivity term
+   if(Parameters::ohmHyperTerm > 0) {
+      Ey_NE += EHyperGrid.get(i-1,j,k-1)->at(fsgrids::ehyper::EYHYPER);
+   }
 
    #ifndef FS_1ST_ORDER_SPACE
       // 2nd order terms:
@@ -1171,6 +1213,7 @@ void calculateEdgeElectricFieldZ(
    FsGrid< std::array<Real, fsgrids::efield::N_EFIELD>, FS_STENCIL_WIDTH> & EGrid,
    FsGrid< std::array<Real, fsgrids::ehall::N_EHALL>, FS_STENCIL_WIDTH> & EHallGrid,
    FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeGrid,
+   FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperGrid,
    FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsGrid,
    FsGrid< std::array<Real, fsgrids::dperb::N_DPERB>, FS_STENCIL_WIDTH> & dPerBGrid,
    FsGrid< std::array<Real, fsgrids::dmoments::N_DMOMENTS>, FS_STENCIL_WIDTH> & dMomentsGrid,
@@ -1299,6 +1342,11 @@ void calculateEdgeElectricFieldZ(
    if(Parameters::ohmGradPeTerm > 0) {
       Ez_SW += EGradPeGrid.get(i,j,k)->at(fsgrids::egradpe::EZGRADPE);
    }
+
+   // Hyperresistivity term
+   if(Parameters::ohmHyperTerm > 0) {
+      Ez_SW += EHyperGrid.get(i,j,k)->at(fsgrids::ehyper::EZHYPER);
+   }
 
    #ifndef FS_1ST_ORDER_SPACE
       // 2nd order terms:
@@ -1357,6 +1405,11 @@ void calculateEdgeElectricFieldZ(
    if(Parameters::ohmGradPeTerm > 0) {
       Ez_SE += EGradPeGrid.get(i-1,j,k)->at(fsgrids::egradpe::EZGRADPE);
    }
+
+   // Hyperresistivity term
+   if(Parameters::ohmHyperTerm > 0) {
+      Ez_SW += EHyperGrid.get(i-1,j,k)->at(fsgrids::ehyper::EZHYPER);
+   }
 
    #ifndef FS_1ST_ORDER_SPACE
       // 2nd order terms:
@@ -1413,6 +1466,11 @@ void calculateEdgeElectricFieldZ(
    if(Parameters::ohmGradPeTerm > 0) {
       Ez_NW += EGradPeGrid.get(i,j-1,k)->at(fsgrids::egradpe::EZGRADPE);
    }
+
+   // Hyperresistivity term
+   if(Parameters::ohmHyperTerm > 0) {
+      Ez_NW += EHyperGrid.get(i,j-1,k)->at(fsgrids::ehyper::EZHYPER);
+   }
 
    #ifndef FS_1ST_ORDER_SPACE
       // 2nd order terms:
@@ -1469,6 +1527,11 @@ void calculateEdgeElectricFieldZ(
    if(Parameters::ohmGradPeTerm > 0) {
       Ez_NE += EGradPeGrid.get(i-1,j-1,k)->at(fsgrids::egradpe::EZGRADPE);
    }
+
+   // Hyperresistivity term
+   if(Parameters::ohmHyperTerm > 0) {
+      Ez_NE += EHyperGrid.get(i-1,j-1,k)->at(fsgrids::ehyper::EZHYPER);
+   }
 
    #ifndef FS_1ST_ORDER_SPACE
       // 2nd order terms:
@@ -1543,6 +1606,7 @@ void calculateElectricField(
    FsGrid< std::array<Real, fsgrids::efield::N_EFIELD>, FS_STENCIL_WIDTH> & EGrid,
    FsGrid< std::array<Real, fsgrids::ehall::N_EHALL>, FS_STENCIL_WIDTH> & EHallGrid,
    FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeGrid,
+   FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperGrid,
    FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsGrid,
    FsGrid< std::array<Real, fsgrids::dperb::N_DPERB>, FS_STENCIL_WIDTH> & dPerBGrid,
    FsGrid< std::array<Real, fsgrids::dmoments::N_DMOMENTS>, FS_STENCIL_WIDTH> & dMomentsGrid,
@@ -1566,6 +1630,7 @@ void calculateElectricField(
          EGrid,
          EHallGrid,
          EGradPeGrid,
+         EHyperGrid,
          momentsGrid,
          dPerBGrid,
          dMomentsGrid,
@@ -1586,6 +1651,7 @@ void calculateElectricField(
          EGrid,
          EHallGrid,
          EGradPeGrid,
+         EHyperGrid,
          momentsGrid,
          dPerBGrid,
          dMomentsGrid,
@@ -1606,6 +1672,7 @@ void calculateElectricField(
          EGrid,
          EHallGrid,
          EGradPeGrid,
+         EHyperGrid,
          momentsGrid,
          dPerBGrid,
          dMomentsGrid,
@@ -1651,6 +1718,8 @@ void calculateUpwindedElectricFieldSimple(
    FsGrid< std::array<Real, fsgrids::ehall::N_EHALL>, FS_STENCIL_WIDTH> & EHallGrid,
    FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeGrid,
    FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeDt2Grid,
+   FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperGrid,
+   FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperDt2Grid,
    FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsGrid,
    FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsDt2Grid,
    FsGrid< std::array<Real, fsgrids::dperb::N_DPERB>, FS_STENCIL_WIDTH> & dPerBGrid,
@@ -1680,7 +1749,14 @@ void calculateUpwindedElectricFieldSimple(
          EGradPeDt2Grid.updateGhostCells();
       }
    }
-   if(P::ohmHallTerm == 0) {
+   if(P::ohmHyperTerm > 0) {
+      if (RKCase == RK_ORDER1 || RKCase == RK_ORDER2_STEP2) {
+         EHyperGrid.updateGhostCells();
+      } else {
+         EHyperDt2Grid.updateGhostCells();
+      }
+   }
+   if(P::ohmHallTerm == 0 && P::ohmHyperTerm == 0) {
       dPerBGrid.updateGhostCells();
    }
    if(P::ohmHallTerm == 0 && P::ohmGradPeTerm == 0 && communicateEGradPeOrMomentsDerivatives) {
@@ -1707,6 +1783,7 @@ void calculateUpwindedElectricFieldSimple(
                      EGrid,
                      EHallGrid,
                      EGradPeGrid,
+                     EHyperGrid,
                      momentsGrid,
                      dPerBGrid,
                      dMomentsGrid,
@@ -1724,6 +1801,7 @@ void calculateUpwindedElectricFieldSimple(
                      EDt2Grid,
                      EHallGrid,
                      EGradPeDt2Grid,
+                     EHyperDt2Grid,
                      momentsDt2Grid,
                      dPerBGrid,
                      dMomentsDt2Grid,

fieldsolver/ldz_electric_field.hpp

@@ -30,6 +30,8 @@ void calculateUpwindedElectricFieldSimple(
    FsGrid< std::array<Real, fsgrids::ehall::N_EHALL>, FS_STENCIL_WIDTH> & EHallGrid,
    FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeGrid,
    FsGrid< std::array<Real, fsgrids::egradpe::N_EGRADPE>, FS_STENCIL_WIDTH> & EGradPeDt2Grid,
+   FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperGrid,
+   FsGrid< std::array<Real, fsgrids::ehyper::N_EHYPER>, FS_STENCIL_WIDTH> & EHyperDt2Grid,
    FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsGrid,
    FsGrid< std::array<Real, fsgrids::moments::N_MOMENTS>, FS_STENCIL_WIDTH> & momentsDt2Grid,
    FsGrid< std::array<Real, fsgrids::dperb::N_DPERB>, FS_STENCIL_WIDTH> & dPerBGrid,