Curl2 matrix BC masks

.azure-pipelines.yml

@@ -44,6 +44,9 @@ jobs:
       # Cartesian 2D
       cartesian_2d:
         WARPX_CMAKE_FLAGS: -DWarpX_DIMS=2 -DWarpX_FFT=ON -DWarpX_PYTHON=ON
+      # Cartesian 2D with PETSc
+      cartesian_petsc_2d:
+        WARPX_CMAKE_FLAGS: -DWarpX_DIMS=2 -DWarpX_FFT=OFF -DWarpX_PYTHON=ON -DAMReX_PETSC=YES
       # Cartesian 3D
       cartesian_3d:
         WARPX_CMAKE_FLAGS: -DWarpX_DIMS=3 -DWarpX_FFT=ON -DWarpX_PYTHON=ON
@@ -108,6 +111,18 @@ jobs:
       #    -DopenPMD_USE_PYTHON=OFF -DBUILD_TESTING=OFF -DBUILD_EXAMPLES=OFF -DBUILD_CLI_TOOLS=OFF
       #  #python3 -m pip install --upgrade openpmd-api
       #fi
+      if [ "${AMReX_PETSC:-YES}" == "TRUE" ]; then
+         export warpx_dir=${PWD}
+         export PETSC_DIR=${warpx_dir}/petsc
+         export PETSC_ARCH=arch-opt
+         # clone and configure PETSc
+         cd ${warpx_dir}
+         git clone -b release https://gitlab.com/petsc/petsc.git petsc
+         cd ${PETSC_DIR}
+         ./configure --with-fortran-bindings=no --with-x=no --with-cc=mpicc --with-fc=mpif90 --with-cxx=mpicxx --with-debugging=1 --download-make --download-cmake --download-superlu --download-superlu_dist
+         make -j 4
+         export WARPX_CMAKE_FLAGS="${WARPX_CMAKE_FLAGS} -DPETSC_DIR=${PETSC_DIR} -DPETSC_ARCH=${PETSC_ARCH}"
+      fi
       if [ "${WARPX_RZ_FFT:-FALSE}" == "TRUE" ]; then
         # BLAS++
         cmake-easyinstall --prefix=/usr/local           \

Examples/Tests/implicit/CMakeLists.txt

@@ -31,6 +31,18 @@ add_warpx_test(
     OFF  # dependency
 )
 
+if(AMReX_PETSC)
+    add_warpx_test(
+        test_2d_curl_curl_petsc_pc  # name
+        2  # dims
+        2  # nprocs
+        inputs_test_2d_curl_curl_petsc_pc  # inputs
+        "analysis_petsc_matrix.py diags/diag1000200"  # analysis
+        "analysis_default_regression.py --path diags/diag1000200"  # checksum
+        OFF  # dependency
+    )
+endif()
+
 add_warpx_test(
     test_2d_theta_implicit_planar_pinch  # name
     2  # dims

Examples/Tests/implicit/analysis_petsc_matrix.py

@@ -0,0 +1,26 @@
+#!/usr/bin/env python3
+
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+#
+# This is a script that analyses the simulation results from
+# and script `inputs_test_2d_curl_curl_petsc_pc`.
+# This simulates a time-varing EM wave injected from an insulator boundary
+# using the theta-implicit solver with PETSc's LU preconditer.
+# Since LU is an exact solver, if the preconditioner matrix is constructed
+# correctly, then here should be 1 Newton and 1 GMRES iteration per time step.
+
+import numpy as np
+
+newton_solver = np.loadtxt("diags/reduced_files/newton_solver.txt", skiprows=1)
+num_steps = newton_solver[-1, 0]
+total_newton_iters = newton_solver[-1, 3]
+total_gmres_iters = newton_solver[-1, 7]
+
+# check that there is 1 Newton and 1 GMRES iteration per step
+print(f"total steps: {num_steps}")
+print(f"total gmres iters: {total_gmres_iters}")
+print(f"total newton iters: {total_newton_iters}")
+assert total_gmres_iters == num_steps
+assert total_newton_iters == num_steps

Examples/Tests/implicit/inputs_test_2d_curl_curl_petsc_pc

@@ -0,0 +1,162 @@
+#################################
+########## CONSTANTS ############
+#################################
+
+###   characteristic plasma values
+my_constants.n0 = 1.e23    # plasma density [1/m^3]
+my_constants.T0 = 1.       # plasma temperature [eV]
+my_constants.R0 = 1.5e-2   # plasma radius [m]
+my_constants.I0 = 200.e3   # pinch current [Amps]
+
+my_constants.AD = 2.01410  # atomic number deuterium
+my_constants.mD = AD*m_u - m_e   # deuterium mass [kg]
+
+###   derived characteristic pinch values
+my_constants.B0 = mu0*I0/(2*pi*R0)  # B at r=R0 [Tesla]
+my_constants.PB0 = B0^2/(2.0*mu0)   # B pressure at r=R0
+my_constants.rho0 = AD*m_u*n0       # mass density [kg/m^3]
+my_constants.u0 = sqrt(PB0/rho0)    # implosion speed [m/s]
+my_constants.t0 = R0/u0             # implosion time [m/s]
+
+###   simulation parameters
+my_constants.tsim = 220.e-9       # sim time = 1.35*t0 [s]
+my_constants.rise_time = t0/20.0  # current rise time [s]
+my_constants.dt = 1.0e-3*1.0e-9   # time step [s]
+my_constants.Nppc_z = 10          # particles per cell per dim
+
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+stop_time = tsim
+max_step = 200
+amr.n_cell =  224 8
+amr.max_grid_size = 8
+amr.max_level = 0
+geometry.dims = 2
+geometry.prob_lo = 0.0 0.0
+geometry.prob_hi = 1.54e-2 0.11e-2
+
+#################################
+####### Boundary condition ######
+#################################
+boundary.field_lo = pmc periodic
+boundary.field_hi = pec_insulator periodic
+boundary.particle_lo = reflecting periodic
+boundary.particle_hi = absorbing periodic
+
+insulator.area_x_hi(y,z) = (-1.0 <= z and z <= 1.0)
+insulator.By_x_hi(y,z,t) = min(t/rise_time,1)*B0
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.serialize_initial_conditions = 1
+warpx.verbose = 1
+warpx.limit_verbose_step = false
+warpx.const_dt = dt
+warpx.use_filter = 0
+
+amrex.fpe_trap_invalid = 1
+amrex.fpe_trap_overflow = 1
+amrex.fpe_trap_zero = 1
+
+algo.maxwell_solver = Yee
+algo.evolve_scheme = "theta_implicit_em"
+
+implicit_evolve.theta = 0.5
+implicit_evolve.max_particle_iterations = 21
+implicit_evolve.particle_tolerance = 1.0e-12
+
+implicit_evolve.nonlinear_solver = "newton"
+newton.verbose = true
+newton.linear_solver= petsc_ksp
+newton.max_iterations = 100
+newton.relative_tolerance = 1.0e-6
+newton.absolute_tolerance = 0.0
+newton.require_convergence = false
+newton.diagnostic_file = "diags/reduced_files/newton_solver.txt"
+newton.diagnostic_interval = 10
+
+gmres.verbose_int = 2
+gmres.max_iterations = 1000
+gmres.relative_tolerance = 1.0e-3
+gmres.absolute_tolerance = 0.0
+
+implicit_evolve.use_mass_matrices_jacobian = true
+implicit_evolve.skip_particle_picard_init = true
+implicit_evolve.use_mass_matrices_pc = true
+jacobian.pc_type = "pc_petsc"
+pc_petsc.verbose = true
+pc_petsc.type = "lu"
+
+particles.max_grid_crossings = 2
+algo.particle_pusher = "higuera"
+algo.particle_shape = 2
+algo.current_deposition = "villasenor"
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = electrons deuterium
+
+electrons.charge = -q_e
+electrons.mass = m_e
+electrons.injection_style = "NUniformPerCell"
+electrons.num_particles_per_cell_each_dim = Nppc_z Nppc_z
+electrons.profile = constant
+electrons.density = n0
+electrons.xmax = 0.0
+electrons.momentum_distribution_type = "gaussian"
+electrons.ux_th = sqrt(T0*q_e/m_e)/clight
+electrons.uy_th = sqrt(T0*q_e/m_e)/clight
+electrons.uz_th = sqrt(T0*q_e/m_e)/clight
+
+deuterium.charge = q_e
+deuterium.mass = mD
+deuterium.injection_style = "NUniformPerCell"
+deuterium.num_particles_per_cell_each_dim = Nppc_z Nppc_z
+deuterium.profile = constant
+deuterium.density = n0
+deuterium.xmax = 0.0
+deuterium.momentum_distribution_type = "gaussian"
+deuterium.ux_th = sqrt(T0*q_e/mD)/clight
+deuterium.uy_th = sqrt(T0*q_e/mD)/clight
+deuterium.uz_th = sqrt(T0*q_e/mD)/clight
+
+###   Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 20
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho divE
+diag1.electrons.variables = x z w ux uy uz
+diag1.deuterium.variables = x z w ux uy uz
+
+###   reduced diagnostics
+warpx.reduced_diags_names = particle_energy field_energy poynting_flux
+
+particle_energy.type = ParticleEnergy
+particle_energy.intervals = 10
+particle_energy.precision = 18
+particle_energy.path = diags/reduced_files/
+
+field_energy.type = FieldEnergy
+field_energy.intervals = 10
+field_energy.precision = 18
+field_energy.path = diags/reduced_files/
+
+poynting_flux.type = FieldPoyntingFlux
+poynting_flux.intervals = 10
+poynting_flux.precision = 18
+poynting_flux.path = diags/reduced_files/
+
+#################################
+############ COLLISION ##########
+#################################
+collisions.collision_names = collision1 collision2 collision3
+
+collision1.species = electrons electrons
+collision2.species = deuterium deuterium
+collision3.species = electrons deuterium
+collision1.CoulombLog = 10.0
+collision2.CoulombLog = 10.0
+collision3.CoulombLog = 10.0

Regression/Checksum/benchmarks_json/test_2d_curl_curl_petsc_pc.json

@@ -0,0 +1,31 @@
+{
+  "lev=0": {
+    "Bx": 9.069162586111297,
+    "By": 19.361810144285307,
+    "Bz": 12.107034717007863,
+    "Ex": 4193780781.9712296,
+    "Ey": 3007724460.40948,
+    "Ez": 4187550905.2162943,
+    "divE": 46810632867690.65,
+    "jx": 0.0,
+    "jy": 0.0,
+    "jz": 0.0,
+    "rho": 0.0
+  },
+  "electrons": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_position_x": 0.0,
+    "particle_position_y": 0.0,
+    "particle_weight": 0.0
+  },
+  "deuterium": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_position_x": 0.0,
+    "particle_position_y": 0.0,
+    "particle_weight": 0.0
+  }
+}

Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H

@@ -101,6 +101,15 @@ public:
         else { return nullptr; }
     }
 
+    /**
+     * \brief Return reference to the multifabs for curl curl mask
+     */
+    [[nodiscard]] virtual const amrex::MultiFab* GetCurl2BCmask (const int lev, const int dir) const
+    {
+        amrex::ignore_unused(lev,dir);
+        return nullptr;
+    }
+
     // This parameter is used for the time-step fraction in the PC for implicit
     // treatment of light waves in the curl-curl MLMG solver.
     // This function should return zero if light waves are not treated implicitly
@@ -242,7 +251,7 @@ protected:
     /**
      * \brief Convert from WarpX FieldBoundaryType to amrex::LinOpBCType
      */
-    [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> convertFieldBCToLinOpBC ( const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>& ) const;
+    [[nodiscard]] amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> convertFieldBCToLinOpBC ( const amrex::Array<FieldBoundaryType,AMREX_SPACEDIM>&, const int bdry_side ) const;
 
 };
 

Source/FieldSolver/ImplicitSolvers/ImplicitSolver.cpp

@@ -52,15 +52,15 @@ const Array<FieldBoundaryType,AMREX_SPACEDIM>& ImplicitSolver::GetFieldBoundaryH
 
 Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::GetLinOpBCLo () const
 {
-    return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryLo());
+    return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryLo(),0);
 }
 
 Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::GetLinOpBCHi () const
 {
-    return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryHi());
+    return convertFieldBCToLinOpBC(m_WarpX->GetFieldBoundaryHi(),1);
 }
 
-Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::convertFieldBCToLinOpBC (const Array<FieldBoundaryType,AMREX_SPACEDIM>& a_fbc) const
+Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::convertFieldBCToLinOpBC (const Array<FieldBoundaryType,AMREX_SPACEDIM>& a_fbc, const int bdry_side) const
 {
     Array<LinOpBCType, AMREX_SPACEDIM> lbc;
     for (auto& bc : lbc) { bc = LinOpBCType::interior; }
@@ -82,10 +82,15 @@ Array<LinOpBCType,AMREX_SPACEDIM> ImplicitSolver::convertFieldBCToLinOpBC (const
             // Also for FieldBoundaryType::PMC
             lbc[i] = LinOpBCType::symmetry;
         } else if (a_fbc[i] == FieldBoundaryType::PECInsulator) {
-            ablastr::warn_manager::WMRecordWarning("Implicit solver",
-                "With PECInsulator, in the Curl-Curl preconditioner Neumann boundary will be used since the full boundary is not yet implemented.",
-                ablastr::warn_manager::WarnPriority::medium);
-            lbc[i] = LinOpBCType::symmetry;
+            const int voltage_driven = m_WarpX->GetPECInsulator_IsESet(i,bdry_side);
+            if (voltage_driven) { // Dirichlet for E
+                lbc[i] = LinOpBCType::Dirichlet;
+            } else { // Dirichlet for B
+                ablastr::warn_manager::WMRecordWarning("Implicit solver with current-driven PECInsulator",
+                    "in the Curl-Curl preconditioner. Neumann boundary will be used since the full boundary is not yet implemented.",
+                    ablastr::warn_manager::WarnPriority::medium);
+                lbc[i] = LinOpBCType::symmetry;
+            }
         } else if (a_fbc[i] == FieldBoundaryType::None) {
             WARPX_ABORT_WITH_MESSAGE("LinOpBCType not set for this FieldBoundaryType");
         } else if (a_fbc[i] == FieldBoundaryType::Open) {
@@ -530,7 +535,8 @@ void ImplicitSolver::InitializeMassMatrices ()
 
     // check that PC is being used by nonlinear solver
     if (m_use_mass_matrices_pc) {
-        if (!m_nlsolver->UsePreconditioner()) {
+        const PreconditionerType pc_type = m_nlsolver->GetPreconditionerType();
+        if (pc_type == PreconditionerType::none) {
             m_use_mass_matrices_pc = false;
         }
     }
@@ -721,7 +727,7 @@ void ImplicitSolver::InitializeMassMatrices ()
     // Set the pointer to mass matrix MultiFab
     if (m_use_mass_matrices_pc) {
         for (int lev = 0; lev < m_num_amr_levels; ++lev) {
-            m_mmpc_mfarrvec.push_back(m_WarpX->m_fields.get_alldirs(FieldType::MassMatrices_PC, 0));
+            m_mmpc_mfarrvec.push_back(m_WarpX->m_fields.get_alldirs(FieldType::MassMatrices_PC, lev));
         }
     }
 
@@ -835,16 +841,21 @@ void ImplicitSolver::SetMassMatricesForPC ( const amrex::Real a_theta_dt )
     using ablastr::fields::Direction;
     using warpx::fields::FieldType;
 
-    // Scale mass matrices used by preconditioner by c^2*mu0*theta*dt and add 1 to diagonal terms
+    // Scale mass matrices used by preconditioner by c^2*mu0*theta*dt.
+    // Add one to diagonal terms when using the curl_curl_mlmg pc_type.
+    // The pc_type petsc already has the one from the curl curl operator
     // Note: This should be done after Sync/communication has been called
 
     const amrex::Real pc_factor = PhysConst::c * PhysConst::c * PhysConst::mu0 * a_theta_dt;
+    const PreconditionerType pc_type = m_nlsolver->GetPreconditionerType();
     const int diag_comp = 0;
     for (int lev = 0; lev < m_num_amr_levels; ++lev) {
         for (int dir = 0 ; dir < 3 ; dir++) {
             amrex::MultiFab* MM_PC = m_WarpX->m_fields.get(FieldType::MassMatrices_PC, Direction{dir}, lev);
             MM_PC->mult(pc_factor, 0, MM_PC->nComp());
-            MM_PC->plus(1.0_rt, diag_comp, 1, 0);
+            if (pc_type == PreconditionerType::pc_curl_curl_mlmg) {
+                MM_PC->plus(1.0_rt, diag_comp, 1, 0);
+            }
         }
     }
 

Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.H

@@ -89,6 +89,16 @@ public:
     // This function should return zero if light waves are not treated implicitly
     [[nodiscard]] virtual amrex::Real  GetThetaForPC () const override { return m_theta; }
 
+    /**
+     * \brief Return reference to the multifabs for curl curl mask in virtual dirs
+     */
+    [[nodiscard]] virtual const amrex::MultiFab* GetCurl2BCmask (const int lev, const int field_dir) const override;
+
+    /**
+     * \brief Inititialize the BC masks for the curl curl operator
+     */
+    void InitializeCurlCurlBCMasks ();
+
 private:
 
     /**