Implement div(B) Cleaning With FDTD (#1829)

* Implement div(B) Cleaning With FDTD

* Add CI Test

* Clean Up

Author: EZoni

Examples/Tests/divb_cleaning/analysis.py

@@ -0,0 +1,49 @@
+#! /usr/bin/env python
+
+# Copyright 2019
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+import sys
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import numpy as np
+import yt
+yt.funcs.mylog.setLevel(50)
+import re
+import checksumAPI
+
+# Name of the last plotfile
+fn = sys.argv[1]
+
+# Load yt data
+ds_old = yt.load('divb_cleaning_3d_plt00398')
+ds_mid = yt.load('divb_cleaning_3d_plt00399')
+ds_new = yt.load(fn) # this is the last plotfile
+
+ad_old = ds_old.covering_grid(level = 0, left_edge = ds_old.domain_left_edge, dims = ds_old.domain_dimensions)
+ad_mid = ds_mid.covering_grid(level = 0, left_edge = ds_mid.domain_left_edge, dims = ds_mid.domain_dimensions)
+ad_new = ds_new.covering_grid(level = 0, left_edge = ds_new.domain_left_edge, dims = ds_new.domain_dimensions)
+
+G_old = ad_old['boxlib', 'G'].v.squeeze()
+G_new = ad_new['boxlib', 'G'].v.squeeze()
+divB  = ad_mid['boxlib', 'divB'].v.squeeze()
+
+# Check max norm of error on div(B) = dG/dt
+# (the time interval between old and new is 2*dt)
+dt = 1.504557189e-15
+x = G_new - G_old
+y = divB * 2 * dt
+
+rel_error = np.amax(abs(x - y)) / np.amax(abs(y))
+tolerance = 1e-1
+
+assert(rel_error < tolerance)
+
+test_name = fn[:-9] # Could also be os.path.split(os.getcwd())[1]
+
+if re.search('single_precision', fn):
+    checksumAPI.evaluate_checksum(test_name, fn, rtol=1.e-3)
+else:
+    checksumAPI.evaluate_checksum(test_name, fn)

Examples/Tests/divb_cleaning/inputs_3d

@@ -0,0 +1,54 @@
+# Iterations
+max_step = 400
+
+# Domain
+amr.n_cell =  32 32 32
+amr.max_grid_size = 16
+amr.max_level = 0
+
+# Geometry
+geometry.coord_sys = 0
+geometry.is_periodic = 0 0 0
+geometry.prob_lo = -50.e-6 -50.e-6 -50.e-6
+geometry.prob_hi =  50.e-6  50.e-6  50.e-6
+
+# Shape factors
+interpolation.nox = 3
+interpolation.noy = 3
+interpolation.noz = 3
+
+# Numerics
+warpx.cfl = 0.25
+warpx.do_divb_cleaning = 1
+warpx.use_filter = 1
+
+# External magnetic field
+my_constants.qm = 1e-1
+warpx.B_ext_grid_init_style = parse_B_ext_grid_function
+warpx.Bx_external_grid_function(x,y,z) = qm * x / (x*x + y*y + z*z)
+warpx.By_external_grid_function(x,y,z) = qm * y / (x*x + y*y + z*z)
+warpx.Bz_external_grid_function(x,y,z) = qm * z / (x*x + y*y + z*z)
+
+# Particle beam
+particles.species_names = beam
+beam.charge = -q_e
+beam.mass = 1.e30
+beam.injection_style = "gaussian_beam"
+beam.x_rms = 2.e-6
+beam.y_rms = 2.e-6
+beam.z_rms = 2.e-6
+beam.x_m = 0.
+beam.y_m = 0.
+beam.z_m = 0.e-6
+beam.npart = 20000
+beam.q_tot = -1.e-20
+beam.momentum_distribution_type = "gaussian"
+beam.ux_m = 0.0
+beam.uy_m = 0.0
+beam.uz_m = 0.5
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 0, 398:400:1
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz divE divB G

Regression/Checksum/benchmarks_json/divb_cleaning_3d.json

@@ -0,0 +1,13 @@
+{
+  "lev=0": {
+    "Bx": 30110529.73244452,
+    "By": 30110529.732444517,
+    "Bz": 30110529.73244452,
+    "Ex": 137103615680453.66,
+    "Ey": 137103615680454.56,
+    "Ez": 137103615680494.5,
+    "G": 0.08248210392635753,
+    "divB": 6944252335584.074,
+    "divE": 60881293.88523142
+  }
+}

Regression/WarpX-tests.ini

@@ -980,6 +980,23 @@ runtime_params = warpx.do_dynamic_scheduling=0
 analysisRoutine = Examples/Tests/initial_plasma_profile/analysis.py
 tolerance = 1.e-14
 
+[divb_cleaning_3d]
+buildDir = .
+inputFile = Examples/Tests/divb_cleaning/inputs_3d
+runtime_params =
+dim = 3
+addToCompileString =
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Tests/divb_cleaning/analysis.py
+tolerance = 1.e-14
+
 [dive_cleaning_2d]
 buildDir = .
 inputFile = Examples/Modules/dive_cleaning/inputs_3d

Source/Diagnostics/Diagnostics.cpp

@@ -62,6 +62,13 @@ Diagnostics::BaseReadParameters ()
             "plot F only works if warpx.do_dive_cleaning = 1");
     }
 
+    // G can be written to file only if WarpX::do_divb_cleaning = 1
+    if (WarpXUtilStr::is_in(m_varnames, "G"))
+    {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+            warpx.do_divb_cleaning, "G can be written to file only if warpx.do_divb_cleaning = 1");
+    }
+
     // If user requests to plot proc_number for a serial run,
     // delete proc_number from fields_to_plot
     if (amrex::ParallelDescriptor::NProcs() == 1){

Source/Diagnostics/FullDiagnostics.cpp

@@ -425,6 +425,8 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
             i++;
         } else if ( m_varnames[comp] == "F" ){
             m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_F_fp(lev), lev, m_crse_ratio);
+        } else if ( m_varnames[comp] == "G" ){
+            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_G_fp(lev), lev, m_crse_ratio);
         } else if ( m_varnames[comp] == "phi" ){
             m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_phi_fp(lev), lev, m_crse_ratio);
         } else if ( m_varnames[comp] == "part_per_cell" ){

Source/Evolve/WarpXEvolve.cpp

@@ -359,7 +359,9 @@ WarpX::OneStep_nosub (Real cur_time)
         }
     } else {
         EvolveF(0.5_rt * dt[0], DtType::FirstHalf);
+        EvolveG(0.5_rt * dt[0], DtType::FirstHalf);
         FillBoundaryF(guard_cells.ng_FieldSolverF);
+        FillBoundaryG(guard_cells.ng_FieldSolverG);
         EvolveB(0.5_rt * dt[0]); // We now have B^{n+1/2}
 
         if (do_silver_mueller) ApplySilverMuellerBoundary( dt[0] );
@@ -377,6 +379,7 @@ WarpX::OneStep_nosub (Real cur_time)
 
         FillBoundaryE(guard_cells.ng_FieldSolver);
         EvolveF(0.5_rt * dt[0], DtType::SecondHalf);
+        EvolveG(0.5_rt * dt[0], DtType::SecondHalf);
         EvolveB(0.5_rt * dt[0]); // We now have B^{n+1}
 
         // Synchronize E and B fields on nodal points

Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt

@@ -7,6 +7,7 @@ target_sources(WarpX
     EvolveEPML.cpp
     EvolveF.cpp
     EvolveFPML.cpp
+    EvolveG.cpp
     FiniteDifferenceSolver.cpp
     MacroscopicEvolveE.cpp
     ApplySilverMuellerBoundary.cpp

Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp

@@ -25,6 +25,7 @@ using namespace amrex;
 void FiniteDifferenceSolver::EvolveB (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+    std::unique_ptr<amrex::MultiFab> const& Gfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
     int lev, amrex::Real const dt ) {
 
@@ -38,21 +39,20 @@ void FiniteDifferenceSolver::EvolveB (
 #else
     if (m_do_nodal) {
 
-        EvolveBCartesian <CartesianNodalAlgorithm> ( Bfield, Efield, face_areas, lev, dt );
+        EvolveBCartesian <CartesianNodalAlgorithm> ( Bfield, Efield, Gfield, face_areas, lev, dt );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
 
-        EvolveBCartesian <CartesianYeeAlgorithm> ( Bfield, Efield, face_areas, lev, dt );
+        EvolveBCartesian <CartesianYeeAlgorithm> ( Bfield, Efield, Gfield, face_areas, lev, dt );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
 
-        EvolveBCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, face_areas, lev, dt );
+        EvolveBCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, Gfield, face_areas, lev, dt );
 
 #endif
     } else {
         amrex::Abort("EvolveB: Unknown algorithm");
     }
-
 }
 
 
@@ -62,6 +62,7 @@ template<typename T_Algo>
 void FiniteDifferenceSolver::EvolveBCartesian (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+    std::unique_ptr<amrex::MultiFab> const& Gfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas,
     int lev, amrex::Real const dt ) {
 
@@ -71,6 +72,8 @@ void FiniteDifferenceSolver::EvolveBCartesian (
 
     amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
 
+    Real constexpr c2 = PhysConst::c * PhysConst::c;
+
     // Loop through the grids, and over the tiles within each grid
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
@@ -138,9 +141,32 @@ void FiniteDifferenceSolver::EvolveBCartesian (
                 Bz(i, j, k) += dt * T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k)
                              - dt * T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k);
             }
-
         );
 
+        // div(B) cleaning correction for errors in magnetic Gauss law (div(B) = 0)
+        if (Gfield)
+        {
+            // Extract field data for this grid/tile
+            Array4<Real> G = Gfield->array(mfi);
+
+            // Loop over cells and update G
+            amrex::ParallelFor(tbx, tby, tbz,
+
+                [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                {
+                    Bx(i,j,k) += c2 * dt * T_Algo::DownwardDx(G, coefs_x, n_coefs_x, i, j, k);
+                },
+                [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                {
+                    By(i,j,k) += c2 * dt * T_Algo::DownwardDy(G, coefs_y, n_coefs_y, i, j, k);
+                },
+                [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                {
+                    Bz(i,j,k) += c2 * dt * T_Algo::DownwardDz(G, coefs_z, n_coefs_z, i, j, k);
+                }
+            );
+        }
+
         if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
         {
             amrex::Gpu::synchronize();

Source/FieldSolver/FiniteDifferenceSolver/EvolveG.cpp

@@ -0,0 +1,95 @@
+/* Copyright 2020
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "Utils/WarpXAlgorithmSelection.H"
+#include "FiniteDifferenceSolver.H"
+#ifdef WARPX_DIM_RZ
+#   include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H"
+#else
+#   include "FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H"
+#   include "FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H"
+#   include "FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H"
+#endif
+#include "Utils/WarpXConst.H"
+#include "WarpX.H"
+#include <AMReX_Gpu.H>
+
+using namespace amrex;
+
+void FiniteDifferenceSolver::EvolveG (
+    std::unique_ptr<amrex::MultiFab>& Gfield,
+    std::array<std::unique_ptr<amrex::MultiFab>,3> const& Bfield,
+    amrex::Real const dt)
+{
+#ifdef WARPX_DIM_RZ
+    // TODO Implement G update equation in RZ geometry
+    amrex::ignore_unused(Gfield, Bfield, dt);
+#else
+    // Select algorithm
+    if (m_do_nodal)
+    {
+        EvolveGCartesian<CartesianNodalAlgorithm>(Gfield, Bfield, dt);
+    }
+    else if (m_fdtd_algo == MaxwellSolverAlgo::Yee)
+    {
+        EvolveGCartesian<CartesianYeeAlgorithm>(Gfield, Bfield, dt);
+    }
+    else if (m_fdtd_algo == MaxwellSolverAlgo::CKC)
+    {
+        EvolveGCartesian<CartesianCKCAlgorithm>(Gfield, Bfield, dt);
+    }
+    else
+    {
+        amrex::Abort("EvolveG: unknown FDTD algorithm");
+    }
+#endif
+}
+
+#ifndef WARPX_DIM_RZ
+
+template<typename T_Algo>
+void FiniteDifferenceSolver::EvolveGCartesian (
+    std::unique_ptr<amrex::MultiFab>& Gfield,
+    std::array<std::unique_ptr<amrex::MultiFab>,3> const& Bfield,
+    amrex::Real const dt)
+{
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+
+    // Loop over grids and over tiles within each grid
+    for (amrex::MFIter mfi(*Gfield, TilingIfNotGPU()); mfi.isValid(); ++mfi)
+    {
+        // Extract field data for this grid/tile
+        amrex::Array4<amrex::Real> const& G = Gfield->array(mfi);
+        amrex::Array4<amrex::Real> const& Bx = Bfield[0]->array(mfi);
+        amrex::Array4<amrex::Real> const& By = Bfield[1]->array(mfi);
+        amrex::Array4<amrex::Real> const& Bz = Bfield[2]->array(mfi);
+
+        // Extract stencil coefficients
+        amrex::Real const* const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr();
+        amrex::Real const* const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr();
+        amrex::Real const* const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr();
+
+        const int n_coefs_x = m_stencil_coefs_x.size();
+        const int n_coefs_y = m_stencil_coefs_y.size();
+        const int n_coefs_z = m_stencil_coefs_z.size();
+
+        // Extract tilebox to loop over
+        amrex::Box const& tf = mfi.tilebox(Gfield->ixType().toIntVect());
+
+        // Loop over cells and update G
+        amrex::ParallelFor(tf, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+        {
+            G(i,j,k) += dt * (T_Algo::UpwardDx(Bx, coefs_x, n_coefs_x, i, j, k)
+                            + T_Algo::UpwardDy(By, coefs_y, n_coefs_y, i, j, k)
+                            + T_Algo::UpwardDz(Bz, coefs_z, n_coefs_z, i, j, k));
+        });
+    }
+}
+
+#endif