From 1fd412a51b662ced98fc952b5317aaa8f2ecd342 Mon Sep 17 00:00:00 2001
From: "S. Eric Clark" <25495882+clarkse@users.noreply.github.com>
Date: Fri, 8 Nov 2024 13:17:15 -0800
Subject: [PATCH] Fixing some bugs that lead to non-convergence. Relaxing
 tolerance which 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.

Signed-off-by: S. Eric Clark <25495882+clarkse@users.noreply.github.com>
---
 .../MagnetostaticSolver.cpp                   | 27 +++++++++++---
 Source/ablastr/fields/VectorPoissonSolver.H   | 37 +++++++++++++++++--
 2 files changed, 55 insertions(+), 9 deletions(-)

diff --git a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp
index 5c28ff1f3c7..d02432e13f3 100644
--- a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp
+++ b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp
@@ -65,7 +65,11 @@ 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.");
+
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(EB::enabled(),
+        "Magnetostatic Solver currently requires an embedded boundary.");
 
     AddMagnetostaticFieldLabFrame();
 }
@@ -128,7 +132,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 +197,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
diff --git a/Source/ablastr/fields/VectorPoissonSolver.H b/Source/ablastr/fields/VectorPoissonSolver.H
index f6dd2a99cf1..d91b5b9a7ce 100644
--- a/Source/ablastr/fields/VectorPoissonSolver.H
+++ b/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