Update the magnetic energy when enforcing consistency in the interface magnetic field

doc/wiki/Runtime-Parameters-related/[Runtime-Parameters]-All.md

@@ -327,7 +327,7 @@ For variables with `Default/Min/Max` labeled as `Depend`, click the parameter na
 | OPT__RES_PHASE                                                                                       |               0 |            None |            None | restriction on phase [0] ##ELBDM ONLY## |
 | [[ OPT__REUSE_MEMORY \| [Runtime-Parameters]-Refinement#OPT__REUSE_MEMORY ]]                         |               2 |               0 |               2 | reuse patch memory to reduce memory fragmentation: (0=off, 1=on, 2=aggressive) [2] |
 | [[ OPT__RHO_INT_SCHEME \| [Runtime-Parameters]-Interpolation#OPT__RHO_INT_SCHEME ]]                  |       INT_CQUAD |               1 |               7 | ghost-zone mass density for the Poisson solver [4] |
-| [[ OPT__SAME_INTERFACE_B \| [Runtime-Parameters]-Hydro#OPT__SAME_INTERFACE_B ]]                      |               0 |            None |            None | ensure B field consistency on the shared interfaces between sibling patches (for debugging) [0] ##MHD ONLY## |
+| [[ OPT__SAME_INTERFACE_B \| [Runtime-Parameters]-Hydro#OPT__SAME_INTERFACE_B ]]                      |              -1 |              -1 |               1 | ensure B field consistency on the shared interfaces between sibling patches (mainly for debugging) (-1=auto, 0=off, 1=on) [-1] ##MHD ONLY## |
 | [[ OPT__SELF_GRAVITY \| [Runtime-Parameters]-Gravity#OPT__SELF_GRAVITY ]]                            |               1 |            None |            None | add self-gravity [1] |
 | [[ OPT__SORT_PATCH_BY_LBIDX \| [Runtime-Parameters]-Miscellaneous#OPT__SORT_PATCH_BY_LBIDX ]]        |          Depend |          Depend |          Depend | sort patches to improve bitwise reproducibility [SERIAL:0, LOAD_BALACNE:1] |
 | [[ OPT__TIMING_BALANCE \| [Runtime-Parameters]-Miscellaneous#OPT__TIMING_BALANCE ]]                  |               0 |            None |            None | record the max/min elapsed time in various code sections for checking load balance [0] |

doc/wiki/Runtime-Parameters-related/[Runtime-Parameters]-Hydro.md

@@ -155,7 +155,8 @@ solutions.
     * **Description:**
 Maximum number of iterations to reduce [MINMOD_COEFF](#MINMOD_COEFF)
 when data reconstruction fails. It improves code stability but
-may break strict conservation.
+may break strict conservation and magnetic field consistency at
+patch interfaces.
     * **Restriction:**
 
 <a name="OPT__1ST_FLUX_CORR"></a>
@@ -196,11 +197,12 @@ Only applicable when enabling the compilation option
 [[--dual | [Installation]-Option-List#--dual]].
 
 <a name="OPT__SAME_INTERFACE_B"></a>
-* #### `OPT__SAME_INTERFACE_B` &ensp; (0=off, 1=on) &ensp; [0]
+* #### `OPT__SAME_INTERFACE_B` &ensp; (-1=set to default, 0=off, 1=on) &ensp; [depend]
     * **Description:**
 Ensure adjacent patches on the same level have exactly the same magnetic field
 on their shared interfaces (including round-off errors).
 This option is mainly for debugging purposes.
+This option is enabled by default when [MINMOD_MAX_ITER](#MINMOD_MAX_ITER) is enabled.
     * **Restriction:**
 
 <a name="OPT__FIXUP_FLUX"></a>

example/test_problem/Template/Input__Parameter

@@ -257,7 +257,7 @@ OPT__1ST_FLUX_CORR           -1           # correct unphysical results (defined
                                           # (<0=auto, 0=off, 1=3D, 2=3D+1D) [-1] ##MHM/MHM_RP/CTU ONLY##
 OPT__1ST_FLUX_CORR_SCHEME    -1           # Riemann solver for OPT__1ST_FLUX_CORR (-1=auto, 0=none, 1=Roe, 2=HLLC, 3=HLLE, 4=HLLD) [-1]
 DUAL_ENERGY_SWITCH            2.0e-2      # apply dual-energy if E_int/E_kin < DUAL_ENERGY_SWITCH [2.0e-2] ##DUAL_ENERGY ONLY##
-OPT__SAME_INTERFACE_B         0           # ensure B field consistency on the shared interfaces between sibling patches (for debugging) [0] ##MHD ONLY##
+OPT__SAME_INTERFACE_B        -1           # ensure B field consistency on the shared interfaces between sibling patches (mainly for debugging) (-1=auto, 0=off, 1=on) [-1] ##MHD ONLY##
 
 
 # fluid solver in ELBDM (MODEL==ELBDM only)

include/Global.h

@@ -133,7 +133,7 @@ extern bool             OPT__FIXUP_ELECTRIC, OPT__CK_INTERFACE_B, OPT__OUTPUT_CC
 extern bool             OPT__OUTPUT_DIVMAG;
 extern int              OPT__CK_DIVERGENCE_B;
 extern double           UNIT_B;
-extern bool             OPT__SAME_INTERFACE_B;
+extern SameInterfaceB_t OPT__SAME_INTERFACE_B;
 
 extern OptInitMagByVecPot_t OPT__INIT_BFIELD_BYVECPOT;
 #endif

include/Prototype.h

@@ -574,7 +574,7 @@ void MHD_AllocateElectricArray( const int lv );
 void MHD_Aux_Check_InterfaceB( const char *comment );
 void MHD_Aux_Check_DivergenceB( const bool Verbose, const char *comment );
 void MHD_FixUp_Electric( const int lv );
-void MHD_SameInterfaceB( const int lv );
+void MHD_SameInterfaceB( const int lv, const int FluSg, const int MagSg );
 void MHD_CopyPatchInterfaceBField( const int lv, const int PID, const int SibID, const int MagSg );
 void MHD_BoundaryCondition_Outflow( real **Array, const int BC_Face, const int NVar, const int GhostSize,
                                     const int ArraySizeX, const int ArraySizeY, const int ArraySizeZ,

include/Typedef.h

@@ -567,6 +567,14 @@ const LoadParaMode_t
    LOAD_HDF5_OUTPUT = 2;
 
 
+// whether to enforce consistency of magnetic field at the patch interface
+typedef int SameInterfaceB_t;
+const SameInterfaceB_t
+   SAME_INTERFACE_B_DEFAULT = -1,
+   SAME_INTERFACE_B_NO      = 0,
+   SAME_INTERFACE_B_YES     = 1;
+
+
 // function pointers
 typedef real (*EoS_GUESS_t)    ( const real Con[], real* const Constant, const double AuxArray_Flt[],
                                  const int AuxArray_Int[], const real *const Table[EOS_NTABLE_MAX] );

src/Init/Init_GAMER.cpp

@@ -246,8 +246,8 @@ void Init_GAMER( int *argc, char ***argv )
 
 // ensure B field consistency on the shared interfaces between sibling patches
 #  if ( MODEL == HYDRO  &&  defined MHD )
-   if ( OPT__SAME_INTERFACE_B )
-   for (int lv=0; lv<NLEVEL; lv++)  MHD_SameInterfaceB( lv );
+   if ( OPT__SAME_INTERFACE_B == SAME_INTERFACE_B_YES )
+   for (int lv=0; lv<NLEVEL; lv++)  MHD_SameInterfaceB( lv, amr->FluSg[lv], amr->MagSg[lv] );
 #  endif
 
 

src/Init/Init_Load_Parameter.cpp

@@ -320,7 +320,7 @@ void Init_Load_Parameter()
    ReadPara->Add( "DUAL_ENERGY_SWITCH",         &DUAL_ENERGY_SWITCH,              2.0e-2,          0.0,           NoMax_double   );
 #  endif
 #  ifdef MHD
-   ReadPara->Add( "OPT__SAME_INTERFACE_B",      &OPT__SAME_INTERFACE_B,           false,           Useless_bool,  Useless_bool   );
+   ReadPara->Add( "OPT__SAME_INTERFACE_B",      &OPT__SAME_INTERFACE_B,          -1,                      -1,     1              );
 #  endif
 #  endif // #if ( MODEL == HYDRO )
 

src/Init/Init_ResetParameter.cpp

@@ -324,6 +324,15 @@ void Init_ResetParameter()
 #  endif
 
 
+#  ifdef MHD
+   if ( OPT__SAME_INTERFACE_B == SAME_INTERFACE_B_DEFAULT )
+   {
+      OPT__SAME_INTERFACE_B = ( MINMOD_MAX_ITER > 0 ) ? SAME_INTERFACE_B_YES : SAME_INTERFACE_B_NO;
+
+      PRINT_RESET_PARA( OPT__SAME_INTERFACE_B, FORMAT_INT, "" );
+   }
+#  endif
+
 // text format parameters
 // --> The current strategy is to read the integer in between % and . to determine the string length.
 //     For example, the format %20.16e will give a length of 20. However, if only checking the string after %,

src/Main/EvolveLevel.cpp

@@ -303,6 +303,18 @@ void EvolveLevel( const int lv, const double dTime_FaLv )
 
       } // if ( OPT__OVERLAP_MPI ) ... else ...
 
+#     if ( MODEL == HYDRO  &&  defined MHD )
+      if ( OPT__SAME_INTERFACE_B == SAME_INTERFACE_B_YES )
+      {
+         TIMING_FUNC(   Buf_GetBufferData( lv, SaveSg_Flu, SaveSg_Mag, NULL_INT, DATA_GENERAL,
+                                           _ENGY, _MAG, Flu_ParaBuf, USELB_YES  ),
+                        Timer_GetBuf[lv][0],   TIMER_ON   );
+
+         TIMING_FUNC(   MHD_SameInterfaceB( lv, SaveSg_Flu, SaveSg_Mag ),
+                        Timer_Flu_Advance[lv],   TIMER_ON   );
+      } // if ( OPT__SAME_INTERFACE_B == SAME_INTERFACE_B_YES )
+#     endif // #if ( MODEL == HYDRO  &&  defined MHD )
+
       amr->FluSg    [lv]             = SaveSg_Flu;
       amr->FluSgTime[lv][SaveSg_Flu] = TimeNew;
 #     ifdef MHD

src/Main/Main.cpp

@@ -123,7 +123,7 @@ bool                 OPT__FIXUP_ELECTRIC, OPT__CK_INTERFACE_B, OPT__OUTPUT_CC_MA
 bool                 OPT__OUTPUT_DIVMAG;
 int                  OPT__CK_DIVERGENCE_B;
 double               UNIT_B;
-bool                 OPT__SAME_INTERFACE_B;
+SameInterfaceB_t     OPT__SAME_INTERFACE_B;
 
 OptInitMagByVecPot_t OPT__INIT_BFIELD_BYVECPOT;
 #endif
@@ -697,14 +697,18 @@ int main( int argc, char *argv[] )
       if ( OPT__CORR_AFTER_ALL_SYNC == CORR_AFTER_SYNC_EVERY_STEP )
       TIMING_FUNC(   Flu_CorrAfterAllSync(),          Timer_Main[6],   TIMER_ON   );
 
+//    it may be unnecessary to call `MHD_SameInterfaceB()` here, since we already call it every sub-step
+//    in EvolveLevel() after the fluid solver (unless something other than the fluid solver can also introduce
+//    inconsistencies in the B field, such as fix-up operations?).
 #     if ( MODEL == HYDRO  &&  defined MHD )
-      if ( OPT__SAME_INTERFACE_B )
+      if ( OPT__SAME_INTERFACE_B == SAME_INTERFACE_B_YES )
       {
          if ( OPT__VERBOSE  &&  MPI_Rank == 0 )
             Aux_Message( stdout, "   MHD_SameInterfaceB                       ... " );
 
          for (int lv=0; lv<NLEVEL; lv++)
-         TIMING_FUNC(   MHD_SameInterfaceB( lv ),     Timer_Main[6],   TIMER_ON   );
+         TIMING_FUNC(   MHD_SameInterfaceB( lv, amr->FluSg[lv], amr->MagSg[lv] ),
+                        Timer_Main[6],   TIMER_ON   );
 
          if ( OPT__VERBOSE  &&  MPI_Rank == 0 )
             Aux_Message( stdout, "done\n" );

src/Model_Hydro/MHD_SameInterfaceB.cpp

@@ -8,7 +8,7 @@
 //-------------------------------------------------------------------------------------------------------
 // Function    :  MHD_SameInterfaceB
 // Description :  Ensure that adjacent patches on the same level (i.e., sibling patches) have *exactly*
-//                the same B field on their shared interfaces
+//                the same B field and B energy on their shared interfaces
 //                --> Even round-off errors are the same
 //
 // Note        :  1. Applied to both real and buffer patches
@@ -17,12 +17,26 @@
 //                4. Always use the B field on the +x/+y/+z sides to overwrite that on the -x/-y/-z sides
 //                5. Mainly for debugging purposes since this consistency should already be guaranteed
 //                   even when disabling OPT__SAME_INTERFACE_B
+//                6. The following two approaches are equivalent. We currently adopt approach (1).
+//                   --> (1) Iterate over both real and buffer patches. In this case, there is no need to call
+//                           Buf_GetBufferData() afterward to exchange buffer-patch data.
+//                       (2) Iterate over real patches only. In this case, we must call Buf_GetBufferData()
+//                           afterward to exchange Flu_ParaBuf buffer-patch data for both the magnetic field and energy
+//                7. Possible optimizations:
+//                   --> (1) When calling MHD_SameInterfaceB() in EvolveLevel(), iterate over real patches only,
+//                           since Buf_GetBufferData() is always called in EvolveLevel() anyway.
+//                       (2) Related to the above, when calling Buf_GetBufferData() before MHD_SameInterfaceB()
+//                           in EvolveLevel(), we could exchange only 0 instead of Flu_ParaBuf data for the
+//                           magnetic field and skip exchanging the energy, since buffer patches do not need to be
+//                           updated inside MHD_SameInterfaceB() here.
 //
-// Parameter   :  lv : AMR level
+// Parameter   :  lv    : AMR level
+//                FluSg : Sandglass to access/store the fluid   data
+//                MagSg : Sandglass to access/store the B field data
 //
 // Return      :  Longitudinal B field on the interfaces between sibling patches
 //-------------------------------------------------------------------------------------------------------
-void MHD_SameInterfaceB( const int lv )
+void MHD_SameInterfaceB( const int lv, const int FluSg, const int MagSg )
 {
 
 // check
@@ -31,27 +45,68 @@ void MHD_SameInterfaceB( const int lv )
       Aux_Error( ERROR_INFO, "incorrect lv = %d !!\n", lv );
 #  endif
 
+// start of OpenMP parallel region
+#  pragma omp parallel
+   {
+      real ***Emag_old = NULL;
 
-   const int MagSg = amr->MagSg[lv];
+      Aux_AllocateArray3D( Emag_old, 3, PS1, PS1 );
 
 // iterate over all real and buffer patches
-#  pragma omp parallel for schedule( runtime )
+#  pragma omp for schedule( runtime )
    for (int PID=0; PID<amr->num[lv]; PID++)
    {
-//    use the B field on the +x/+y/+z sides to overwrite that on the -x/-y/-z sides
+//    always use the B field on the -x/-y/-z sides to overwrite that on the +x/+y/+z sides
 //    --> skip s=1/3/5
       for (int s=0; s<6; s+=2)
       {
+         const int d      = s/2;
          const int SibPID = amr->patch[0][lv][PID]->sibling[s];
 
+         if ( SibPID < 0 )   continue;
+
 //       some buffer patches may have magnetic == NULL --> skip them
-         if ( SibPID >= 0  &&
-              amr->patch[MagSg][lv][   PID]->magnetic != NULL  &&
-              amr->patch[MagSg][lv][SibPID]->magnetic != NULL     )
-            MHD_CopyPatchInterfaceBField( lv, PID, s, MagSg );
-      }
+         if ( amr->patch[MagSg][lv][   PID]->magnetic == NULL  ||
+              amr->patch[MagSg][lv][SibPID]->magnetic == NULL     )   continue;
+
+         int ii, jj, kk;
+
+         for (int j=0; j<PS1; j++)
+         for (int i=0; i<PS1; i++)
+         {
+            switch ( d )
+            {
+               case 0: ii = 0;  jj = i;  kk = j;  break;
+               case 1: ii = j;  jj = 0;  kk = i;  break;
+               case 2: ii = i;  jj = j;  kk = 0;  break;
+            }
+
+            Emag_old[d][j][i] = MHD_GetCellCenteredBEnergyInPatch( lv, PID, ii, jj, kk, MagSg );
+         }
+
+         MHD_CopyPatchInterfaceBField( lv, SibPID, s+1, MagSg );
+
+         for (int j=0; j<PS1; j++)
+         for (int i=0; i<PS1; i++)
+         {
+            switch ( d )
+            {
+               case 0: ii = 0;  jj = i;  kk = j;  break;
+               case 1: ii = j;  jj = 0;  kk = i;  break;
+               case 2: ii = i;  jj = j;  kk = 0;  break;
+            }
+
+            const real Emag = MHD_GetCellCenteredBEnergyInPatch( lv, PID, ii, jj, kk, MagSg );
+
+            amr->patch[FluSg][lv][PID]->fluid[ENGY][kk][jj][ii] += (Emag - Emag_old[d][j][i]);
+         }
+      } // for (int s=0; s<6; s+=2)
    } // for (int PID=0; PID<amr->num[lv]; PID++)
 
+   Aux_DeallocateArray3D( Emag_old );
+
+   } // end of OpenMP parallel region
+
 } // FUNCTION : MHD_SameInterfaceB