Add external particle fields ohms law hybrid

Python/pywarpx/HybridPICModel.py

@@ -9,3 +9,4 @@
 from .Bucket import Bucket
 
 hybridpicmodel = Bucket("hybrid_pic_model")
+external_vector_potential = Bucket("external_vector_potential")

Python/pywarpx/WarpX.py

@@ -20,7 +20,7 @@
 from .Diagnostics import diagnostics, reduced_diagnostics
 from .EB2 import eb2
 from .Geometry import geometry
-from .HybridPICModel import hybridpicmodel
+from .HybridPICModel import external_vector_potential, hybridpicmodel
 from .Interpolation import interpolation
 from .Lasers import lasers, lasers_list
 from .Particles import particles, particles_list
@@ -46,6 +46,7 @@ def create_argv_list(self, **kw):
         argv += amrex.attrlist()
         argv += geometry.attrlist()
         argv += hybridpicmodel.attrlist()
+        argv += external_vector_potential.attrlist()
         argv += boundary.attrlist()
         argv += algo.attrlist()
         argv += interpolation.attrlist()

Python/pywarpx/__init__.py

@@ -33,7 +33,7 @@
 from .Diagnostics import diagnostics, reduced_diagnostics  # noqa
 from .EB2 import eb2  # noqa
 from .Geometry import geometry  # noqa
-from .HybridPICModel import hybridpicmodel  # noqa
+from .HybridPICModel import hybridpicmodel, external_vector_potential  # noqa
 from .Interpolation import interpolation  # noqa
 from .Lasers import lasers  # noqa
 from .LoadThirdParty import load_cupy  # noqa

Python/pywarpx/fields.py

@@ -704,6 +704,87 @@ def BzFPExternalWrapper(level=0, include_ghosts=False):
     )
 
 
+def AxHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="hybrid_A_fp_external",
+        idir=0,
+        level=level,
+        include_ghosts=include_ghosts,
+    )
+
+
+def AyHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="hybrid_A_fp_external",
+        idir=1,
+        level=level,
+        include_ghosts=include_ghosts,
+    )
+
+
+def AzHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="hybrid_A_fp_external",
+        idir=2,
+        level=level,
+        include_ghosts=include_ghosts,
+    )
+
+
+def ExHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="hybrid_E_fp_external",
+        idir=0,
+        level=level,
+        include_ghosts=include_ghosts,
+    )
+
+
+def EyHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="hybrid_E_fp_external",
+        idir=1,
+        level=level,
+        include_ghosts=include_ghosts,
+    )
+
+
+def EzHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="hybrid_E_fp_external",
+        idir=2,
+        level=level,
+        include_ghosts=include_ghosts,
+    )
+
+
+def BxHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="hybrid_B_fp_external",
+        idir=0,
+        level=level,
+        include_ghosts=include_ghosts,
+    )
+
+
+def ByHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="hybrid_B_fp_external",
+        idir=1,
+        level=level,
+        include_ghosts=include_ghosts,
+    )
+
+
+def BzHybridExternalWrapper(level=0, include_ghosts=False):
+    return _MultiFABWrapper(
+        mf_name="hybrid_B_fp_external",
+        idir=2,
+        level=level,
+        include_ghosts=include_ghosts,
+    )
+
+
 def JxFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(
         mf_name="current_fp", idir=0, level=level, include_ghosts=include_ghosts

Python/pywarpx/picmi.py

@@ -1782,6 +1782,9 @@ class HybridPICSolver(picmistandard.base._ClassWithInit):
 
     Jx/y/z_external_function: str
         Function of space and time specifying external (non-plasma) currents.
+
+    Ex/y/z_external_function: str
+        Function of space and time specifying external (non-plasma) E-fields.
     """
 
     def __init__(
@@ -1797,6 +1800,7 @@ def __init__(
         Jx_external_function=None,
         Jy_external_function=None,
         Jz_external_function=None,
+        A_external=None,
         **kw,
     ):
         self.grid = grid
@@ -1815,6 +1819,25 @@ def __init__(
         self.Jy_external_function = Jy_external_function
         self.Jz_external_function = Jz_external_function
 
+        self.add_external_fields = None
+
+        # It is expected that a nested dicitonary will be passed
+        # into picmi for each field that has different timings
+        # e.g.
+        # A_external = {
+        #     '<field_name1>': {
+        #         'Ax_external_function': ...,
+        #         'Ax_external_function': ...,
+        #         'Ax_external_function': ...,
+        #         'A_time_external_function': ...
+        #     },
+        #     '<field_name2>: {...}'
+        # }
+        self.A_external = A_external
+
+        if A_external is not None:
+            self.add_external_fields = True
+
         # Handle keyword arguments used in expressions
         self.user_defined_kw = {}
         for k in list(kw.keys()):
@@ -1863,6 +1886,38 @@ def solver_initialize_inputs(self):
                 self.Jz_external_function, self.mangle_dict
             ),
         )
+        pywarpx.hybridpicmodel.add_external_fields = self.add_external_fields
+        pywarpx.external_vector_potential.__setattr__(
+            "fields",
+            pywarpx.my_constants.mangle_expression(
+                list(self.A_external.keys()), self.mangle_dict
+            ),
+        )
+        for field_name, field_dict in self.A_external.items():
+            pywarpx.external_vector_potential.__setattr__(
+                f"{field_name}.Ax_external_grid_function(x,y,z)",
+                pywarpx.my_constants.mangle_expression(
+                    field_dict["Ax_external_function"], self.mangle_dict
+                ),
+            )
+            pywarpx.external_vector_potential.__setattr__(
+                f"{field_name}.Ay_external_grid_function(x,y,z)",
+                pywarpx.my_constants.mangle_expression(
+                    field_dict["Ay_external_function"], self.mangle_dict
+                ),
+            )
+            pywarpx.external_vector_potential.__setattr__(
+                f"{field_name}.Az_external_grid_function(x,y,z)",
+                pywarpx.my_constants.mangle_expression(
+                    field_dict["Az_external_function"], self.mangle_dict
+                ),
+            )
+            pywarpx.external_vector_potential.__setattr__(
+                f"{field_name}.A_time_external_function(t)",
+                pywarpx.my_constants.mangle_expression(
+                    field_dict["A_time_external_function"], self.mangle_dict
+                ),
+            )
 
 
 class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver):

Source/EmbeddedBoundary/CMakeLists.txt

@@ -7,5 +7,6 @@ foreach(D IN LISTS WarpX_DIMS)
         WarpXFaceExtensions.cpp
         WarpXFaceInfoBox.H
         Enabled.cpp
+        Covered.cpp
     )
 endforeach()

Source/EmbeddedBoundary/Covered.H

@@ -0,0 +1,116 @@
+/* Copyright 2024 S. Eric Clark (Helion Energy)
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_EB_COVERED_H_
+#define WARPX_EB_COVERED_H_
+
+
+#include "EmbeddedBoundary/Enabled.H"
+#include "Fields.H"
+
+#  include <AMReX.H>
+#  include <AMReX_Array.H>
+#  include <AMReX_Array4.H>
+#  include <AMReX_Config.H>
+#  include <AMReX_GpuControl.H>
+#  include <AMReX_GpuDevice.H>
+#  include <AMReX_GpuQualifiers.H>
+#  include <AMReX_MFIter.H>
+#  include <AMReX_MultiFab.H>
+#  include <AMReX_REAL.H>
+#  include <AMReX_SPACE.H>
+
+namespace EB {
+
+AMREX_ENUM(CoverTopology,
+    none,
+    face,
+    edge
+);
+
+class Covered
+{
+public:
+    amrex::Array4<amrex::Real> lx, ly, lz, Sx, Sy, Sz;
+
+#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
+    amrex::Dim3 lx_lo, lx_hi, lz_lo, lz_hi;
+#endif
+
+    Covered (amrex::MFIter &mfi, int lev);
+
+    // Can have topology of f or e for 'face' or 'edge'
+    [[nodiscard]] AMREX_GPU_HOST_DEVICE
+    bool isCovered (
+        [[maybe_unused]] const int idir,
+        [[maybe_unused]] const CoverTopology topology,
+        [[maybe_unused]] const int i,
+        [[maybe_unused]] const int j,
+        [[maybe_unused]] const int k) const noexcept
+    {
+#ifdef AMREX_USE_EB
+        bool is_covered = false;
+
+        switch (idir)
+        {
+        case 0:
+#ifdef WARPX_DIM_3D
+            if((topology==CoverTopology::edge and lx(i, j, k)<=0) or (topology==CoverTopology::face and Sx(i, j, k)<=0))
+            {
+                is_covered = true;
+            }
+#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
+            //In XZ and RZ Ex is associated with a x-edge, while Bx is associated with a z-edge
+            if((topology==CoverTopology::edge and lx(i, j, k)<=0) or (topology==CoverTopology::face and lz(i, j, k)<=0))
+            {
+                is_covered = true;
+            }
+#endif
+            break;
+
+        case 1:
+#ifdef WARPX_DIM_3D
+            if((topology==CoverTopology::edge and ly(i, j, k)<=0) or (topology==CoverTopology::face and Sy(i, j, k)<=0))
+            {
+                is_covered = true;
+            }
+#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
+            //In XZ and RZ Ey is associated with a mesh node, so we need to check if  the mesh node is covered
+            if( (topology==CoverTopology::edge and (lx(std::min(i  , lx_hi.x), std::min(j  , lx_hi.y), k)<=0
+                            || lx(std::max(i-1, lx_lo.x), std::min(j  , lx_hi.y), k)<=0
+                            || lz(std::min(i  , lz_hi.x), std::min(j  , lz_hi.y), k)<=0
+                            || lz(std::min(i  , lz_hi.x), std::max(j-1, lz_lo.y), k)<=0)) or
+                (topology==CoverTopology::face and Sy(i,j,k)<=0))
+            {
+                is_covered = true;
+            }
+#endif
+            break;
+
+        case 2:
+#ifdef WARPX_DIM_3D
+            if((topology==CoverTopology::edge and lz(i, j, k)<=0) or (topology==CoverTopology::face and Sz(i, j, k)<=0))
+            {
+                is_covered = true;
+            }
+#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
+            //In XZ and RZ Ez is associated with a z-edge, while Bz is associated with a x-edge
+            if((topology==CoverTopology::edge and lz(i, j, k)<=0) or (topology==CoverTopology::face and lx(i, j, k)<=0))
+            {
+                is_covered = true;
+            }
+#endif
+            break;
+        }
+        return is_covered;
+#else
+        return false;
+#endif //AMREX_USE_EB
+    }
+};
+
+} // namespace EB
+#endif

Source/EmbeddedBoundary/Covered.cpp

@@ -0,0 +1,51 @@
+/* Copyright 2024 S. Eric Clark
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include "WarpX.H"
+
+#include "EmbeddedBoundary/Covered.H"
+#include "EmbeddedBoundary/Enabled.H"
+#include "Fields.H"
+
+#  include <AMReX.H>
+#  include <AMReX_Array.H>
+#  include <AMReX_Array4.H>
+#  include <AMReX_Loop.H>
+#  include <AMReX_MFIter.H>
+#  include <AMReX_MultiFab.H>
+#  include <AMReX_REAL.H>
+#  include <AMReX_SPACE.H>
+
+using namespace amrex;
+using namespace warpx::fields;
+
+namespace EB {
+
+// Default Constructor
+Covered::Covered (amrex::MFIter &mfi, int lev)
+{
+    if (EB::enabled()) {
+        auto& warpx = WarpX::GetInstance();
+        auto edge_lengths = warpx.m_fields.get_alldirs(FieldType::edge_lengths, lev);
+        auto face_areas = warpx.m_fields.get_alldirs(FieldType::face_areas, lev);
+
+        lx = edge_lengths[0]->array(mfi);
+        ly = edge_lengths[1]->array(mfi);
+        lz = edge_lengths[2]->array(mfi);
+        Sx = face_areas[0]->array(mfi);
+        Sy = face_areas[1]->array(mfi);
+        Sz = face_areas[2]->array(mfi);
+
+#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
+        lx_lo = amrex::lbound(lx);
+        lx_hi = amrex::ubound(lx);
+        lz_lo = amrex::lbound(lz);
+        lz_hi = amrex::ubound(lz);
+#endif
+    }
+}
+
+}

Source/EmbeddedBoundary/Make.package

@@ -1,10 +1,12 @@
 CEXE_headers += Enabled.H
+CEXE_headers += Covered.H
 CEXE_headers += ParticleScraper.H
 CEXE_headers += ParticleBoundaryProcess.H
 CEXE_headers += DistanceToEB.H
 CEXE_headers += WarpXFaceInfoBox.H
 
 CEXE_sources += Enabled.cpp
+CEXE_sources += Covered.cpp
 CEXE_sources += WarpXInitEB.cpp
 CEXE_sources += WarpXFaceExtensions.cpp
 

Source/FieldSolver/FiniteDifferenceSolver/CMakeLists.txt

@@ -3,6 +3,7 @@ foreach(D IN LISTS WarpX_DIMS)
     target_sources(lib_${SD}
       PRIVATE
         ComputeDivE.cpp
+        ComputeCurlA.cpp
         EvolveB.cpp
         EvolveBPML.cpp
         EvolveE.cpp