Fixing some bugs that lead to non-convergence. Relaxing tolerance whi… (

#5446)

…ch is currently hard coded, and adding assert to force EB being enabled
in order to avoid a seg fault in AMReX when computing the gradient
solution.

This PR partially addresses
#5444. This PR adds
semi-coarsening in 3D and then adds an assert to keep the magnetostatic
solver from being run without an EB. This is required since in AMReX
MLMG->getGradSolution will segfault when not using an EB.

It should also be noted that
#5175 will use a different scheme
around the embedded boundaries to compute gradients, and will likely
mitigate these issues.

A work around in RZ to use the outer edge is to enable the embedded
boundary and set the boundary radius larger than the outer grid radius.
This works like it would without an embedded boundary and can be used
until either the refactor or the bugfix in AMReX for getGradSolution.

---------

Signed-off-by: S. Eric Clark <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp

@@ -65,7 +65,16 @@ WarpX::ComputeMagnetostaticField()
     // Fields have been reset in Electrostatic solver for this time step, these fields
     // are added into the B fields after electrostatic solve
 
-    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(this->max_level == 0, "Magnetostatic solver not implemented with mesh refinement.");
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(this->max_level == 0,
+        "Magnetostatic solver not implemented with mesh refinement.");
+
+#if defined(AMREX_USE_EB)
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(EB::enabled(),
+        "Magnetostatic Solver currently requires an embedded boundary to be installed for "
+        "compatibility with AMReX when compiling with EB support. "
+        "Current workaround is to install an EB outside of domain or recompile with EB support off."
+        "Workaround for https://github.com/AMReX-Codes/amrex/issues/4223");
+#endif
 
     AddMagnetostaticFieldLabFrame();
 }
@@ -128,7 +137,13 @@ WarpX::AddMagnetostaticFieldLabFrame()
     // const amrex::Real magnetostatic_absolute_tolerance = self_fields_absolute_tolerance*PhysConst::c;
     // temporary fix!!!
     const amrex::Real magnetostatic_absolute_tolerance = 0.0;
-    const amrex::Real self_fields_required_precision = 1e-12;
+    amrex::Real self_fields_required_precision;
+    if constexpr (std::is_same<Real, float>::value) {
+        self_fields_required_precision = 1e-5;
+    }
+    else {
+        self_fields_required_precision = 1e-11;
+    }
     const int self_fields_max_iters = 200;
     const int self_fields_verbosity = 2;
 
@@ -187,11 +202,16 @@ WarpX::computeVectorPotential (ablastr::fields::MultiLevelVectorField const& cur
     });
 
 #if defined(AMREX_USE_EB)
-    amrex::Vector<amrex::EBFArrayBoxFactory const *> factories;
-    for (int lev = 0; lev <= finest_level; ++lev) {
-        factories.push_back(&WarpX::fieldEBFactory(lev));
+    std::optional<amrex::Vector<amrex::EBFArrayBoxFactory const *> > eb_farray_box_factory;
+    auto &warpx = WarpX::GetInstance();
+
+    if (EB::enabled()) {
+        amrex::Vector<amrex::EBFArrayBoxFactory const *> factories;
+        for (int lev = 0; lev <= finest_level; ++lev) {
+            factories.push_back(&warpx.fieldEBFactory(lev));
+        }
+        eb_farray_box_factory = factories;
     }
-    const std::optional<amrex::Vector<amrex::EBFArrayBoxFactory const *> > eb_farray_box_factory({factories});
 #else
     const std::optional<amrex::Vector<amrex::FArrayBoxFactory const *> > eb_farray_box_factory;
 #endif

Source/ablastr/fields/VectorPoissonSolver.H

@@ -1,4 +1,4 @@
-/* Copyright 2022 S. Eric Clark, LLNL
+/* Copyright 2022-2024 S. Eric Clark (Helion Energy, formerly LLNL)
  *
  * This file is part of WarpX.
  *
@@ -137,10 +137,41 @@ computeVectorPotential   (  amrex::Vector<amrex::Array<amrex::MultiFab*, 3> > co
         );
     }
 
-    const amrex::LPInfo& info = amrex::LPInfo();
-
     // Loop over dimensions of A to solve each component individually
     for (int lev=0; lev<=finest_level; lev++) {
+        amrex::LPInfo info;
+
+#ifdef WARPX_DIM_RZ
+        constexpr bool is_rz = true;
+#else
+        constexpr bool is_rz = false;
+#endif
+
+        amrex::Array<amrex::Real,AMREX_SPACEDIM> const dx
+            {AMREX_D_DECL(geom[lev].CellSize(0),
+                          geom[lev].CellSize(1),
+                          geom[lev].CellSize(2))};
+
+
+        if (!eb_enabled && !is_rz) {
+            // Determine whether to use semi-coarsening
+            int max_semicoarsening_level = 0;
+            int semicoarsening_direction = -1;
+            const auto min_dir = static_cast<int>(std::distance(dx.begin(),
+                                                                std::min_element(dx.begin(), dx.end())));
+            const auto max_dir = static_cast<int>(std::distance(dx.begin(),
+                                                                std::max_element(dx.begin(), dx.end())));
+            if (dx[max_dir] > dx[min_dir]) {
+                semicoarsening_direction = max_dir;
+                max_semicoarsening_level = static_cast<int>(std::log2(dx[max_dir] / dx[min_dir]));
+            }
+            if (max_semicoarsening_level > 0) {
+                info.setSemicoarsening(true);
+                info.setMaxSemicoarseningLevel(max_semicoarsening_level);
+                info.setSemicoarseningDirection(semicoarsening_direction);
+            }
+        }
+
         amrex::MLEBNodeFDLaplacian linopx, linopy, linopz;
         if (eb_enabled) {
 #ifdef AMREX_USE_EB