Staircased embedded boundaries in the YEE solver

Examples/Modules/embedded_boundary_cube/analysis_fields.py

@@ -0,0 +1,99 @@
+#! /usr/bin/env python
+
+import yt
+import os, sys
+from scipy.constants import mu_0, pi, c
+import numpy as np
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+# This is a script that analyses the simulation results from
+# the script `inputs_3d`. This simulates a TMmnp mode in a PEC cubic resonator.
+# The magnetic field in the simulation is given (in theory) by:
+# $$ B_x = \frac{-2\mu}{h^2}\, k_x k_z \sin(k_x x)\cos(k_y y)\cos(k_z z)\cos( \omega_p t)$$
+# $$ B_y = \frac{-2\mu}{h^2}\, k_y k_z \cos(k_x x)\sin(k_y y)\cos(k_z z)\cos( \omega_p t)$$
+# $$ B_z = \cos(k_x x)\cos(k_y y)\sin(k_z z)\sin( \omega_p t)$$
+# with
+# $$ h^2 = k_x^2 + k_y^2 + k_z^2$$
+# $$ k_x = \frac{m\pi}{L}$$
+# $$ k_y = \frac{n\pi}{L}$$
+# $$ k_z = \frac{p\pi}{L}$$
+
+hi = [0.8, 0.8, 0.8]
+lo = [-0.8, -0.8, -0.8]
+ncells = [48, 48, 48]
+dx = (hi[0] - lo[0])/ncells[0]
+dy = (hi[1] - lo[1])/ncells[1]
+dz = (hi[2] - lo[2])/ncells[2]
+m = 0
+n = 1
+p = 1
+Lx = 1
+Ly = 1
+Lz = 1
+h_2 = (m * pi / Lx) ** 2 + (n * pi / Ly) ** 2 + (p * pi / Lz) ** 2
+t = 1.3342563807926085e-08
+
+# Compute the analytic solution
+Bx_th = np.zeros(ncells)
+By_th = np.zeros(ncells)
+Bz_th = np.zeros(ncells)
+for i in range(ncells[0]):
+    for j in range(ncells[1]):
+        for k in range(ncells[2]):
+            x = i*dx + lo[0]
+            y = j*dy + lo[1]
+            z = k*dz + lo[2]
+
+            Bx_th[i, j, k] = -2/h_2*mu_0*(m * pi/Lx)*(p * pi/Lz) * (np.sin(m * pi/Lx * (x - Lx/2)) *
+                                                                    np.sin(n * pi/Ly * (y - Ly/2)) *
+                                                                    np.cos(p * pi/Lz * (z - Lz/2)) *
+                                                                    (-Lx/2 < x < Lx/2) *
+                                                                    (-Ly/2 < y < Ly/2) *
+                                                                    (-Lz/2 < z < Lz/2) *
+                                                                    np.cos(np.sqrt(2) *
+                                                                           np.pi / Lx * c * t))
+
+            By_th[i, j, k] = -2/h_2*mu_0*(n * pi/Ly)*(p * pi/Lz) * (np.cos(m * pi/Lx * (x - Lx/2)) *
+                                                                    np.sin(n * pi/Ly * (y - Ly/2)) *
+                                                                    np.cos(p * pi/Lz * (z - Lz/2)) *
+                                                                    (-Lx/2 < x < Lx/2) *
+                                                                    (-Ly/2 < y < Ly/2) *
+                                                                    (-Lz/2 < z < Lz/2) *
+                                                                    np.cos(np.sqrt(2) *
+                                                                    np.pi / Lx * c * t))
+
+            Bz_th[i, j, k] = mu_0*(n * pi/Ly)*(p * pi/Lz) * (np.cos(m * pi/Lx * (x - Lx/2)) *
+                                                             np.cos(n * pi/Ly * (y - Ly/2)) *
+                                                             np.sin(p * pi/Lz * (z - Lz/2)) *
+                                                             (-Lx/2 < x < Lx/2) *
+                                                             (-Ly/2 < y < Ly/2) *
+                                                             (-Lz/2 < z < Lz/2) *
+                                                             np.cos(np.sqrt(2) *
+                                                                    np.pi / Lx * c * t))
+
+
+# Open the right plot file
+filename = sys.argv[1]
+ds = yt.load(filename)
+data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+
+tol_err = 1e-5
+
+# Compute l^2 error on Bx
+Bx_sim = data['Bx'].to_ndarray()
+err_x = np.sqrt(np.sum(np.square(Bx_sim - Bx_th))*dx*dy*dz)
+assert(err_x < tol_err)
+# Compute l^2 error on By
+By_sim = data['By'].to_ndarray()
+err_y = np.sqrt(np.sum(np.square(By_sim - By_th))*dx*dy*dz)
+assert(err_y < tol_err)
+
+# Compute l^2 error on Bz
+Bz_sim = data['Bz'].to_ndarray()
+err_z = np.sqrt(np.sum(np.square(Bz_sim - Bz_th))*dx*dy*dz)
+assert(err_z < tol_err)
+
+test_name = os.path.split(os.getcwd())[1]
+
+checksumAPI.evaluate_checksum(test_name, filename)

Examples/Modules/embedded_boundary_cube/inputs_3d

@@ -0,0 +1,37 @@
+stop_time = 1.3342563807926085e-08
+amr.n_cell = 48 48 48
+amr.max_grid_size = 128
+amr.max_level = 0
+
+geometry.coord_sys   = 0
+geometry.is_periodic = 0 0 0
+geometry.prob_lo     = -0.8 -0.8 -0.8
+geometry.prob_hi     =  0.8  0.8  0.8
+warpx.do_pml = 0
+warpx.const_dt = 1e-6
+warpx.cfl = 1
+
+eb2.geom_type = box
+eb2.box_lo = -0.5 -0.5 -0.5
+eb2.box_hi = 0.5 0.5 0.5
+eb2.box_has_fluid_inside = true
+
+warpx.B_ext_grid_init_style = parse_B_ext_grid_function
+my_constants.m = 0
+my_constants.n = 1
+my_constants.p = 1
+my_constants.pi = 3.141592653589793
+my_constants.Lx = 1
+my_constants.Ly = 1
+my_constants.Lz = 1
+my_constants.h_2 = (m * pi / Lx) ** 2 + (n * pi / Ly) ** 2 + (p * pi / Lz) ** 2
+my_constants.mu_0 = 1.25663706212e-06
+
+warpx.By_external_grid_function(x,y,z) = -2/h_2 * (n * pi / Ly) * (p * pi / Lz) * cos(m * pi / Lx * (x - Lx / 2)) * sin(n * pi / Ly * (y - Ly / 2)) * cos(p * pi / Lz * (z - Lz / 2))*mu_0*(x>-Lx/2)*(x<Lx/2)*(y>-Ly/2)*(y<Ly/2)*(z>-Lz/2)*(z<Lz/2)
+warpx.Bx_external_grid_function(x,y,z) = -2/h_2 * (m * pi / Lx) * (p * pi / Lz) * sin(m * pi / Lx * (x - Lx / 2)) * cos(n * pi / Ly * (y - Ly / 2)) * cos(p * pi / Lz * (z - Lz / 2))*mu_0*(x>-Lx/2)*(x<Lx/2)*(y>-Ly/2)*(y<Ly/2)*(z>-Lz/2)*(z<Lz/2)
+warpx.Bz_external_grid_function(x,y,z) = cos(m * pi / Lx * (x - Lx / 2)) * cos(n * pi / Ly * (y - Ly / 2)) * sin(p * pi / Lz * (z - Lz / 2))*mu_0*(x>-0.5)*(x<0.5)*(y>-0.5)*(y<0.5)*(z>-0.5)*(z<0.5)
+
+diagnostics.diags_names = diag1
+diag1.intervals = 1
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz

Regression/Checksum/benchmarks_json/embedded_boundary_cube.json

@@ -0,0 +1,11 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 0.0066283741197868335,
+    "Bz": 0.006628374119786833,
+    "Ex": 5102618.47115243,
+    "Ey": 0.0,
+    "Ez": 0.0
+  }
+}
+

Regression/WarpX-tests.ini

@@ -2002,3 +2002,19 @@ compileTest = 0
 doVis = 0
 compareParticles = 0
 analysisRoutine = Examples/Modules/resampling/analysis_leveling_thinning.py
+
+[embedded_boundary_cube]
+buildDir = .
+inputFile = Examples/Modules/embedded_boundary_cube/inputs_3d
+runtime_params =
+dim = 3
+addToCompileString = USE_EB=TRUE
+restartTest = 0
+useMPI = 1
+numprocs = 4
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0s
+analysisRoutine = Examples/Modules/embedded_boundary_cube/analysis_fields.py

Source/EmbeddedBoundary/WarpXInitEB.cpp

@@ -15,8 +15,162 @@ WarpX::InitEB ()
 
     amrex::ParmParse pp_eb2("eb2");
     if (!pp_eb2.contains("geom_type")) {
-        pp_eb2.add("geom_type", "all_regular"); // use all_regular by default
+        std::string geom_type = "all_regular";
+        pp_eb2.add("geom_type", geom_type); // use all_regular by default
     }
     amrex::EB2::Build(Geom(maxLevel()), maxLevel(), maxLevel());
+
+#endif
+}
+
+/**
+ * \brief Compute the length of the mesh edges. Here the length is a value in [0, 1].
+ *        An edge of length 0 is fully covered.
+ */
+void
+WarpX::ComputeEdgeLengths () {
+#ifdef AMREX_USE_EB
+    BL_PROFILE("ComputeEdgeLengths");
+
+    auto const eb_fact = fieldEBFactory(maxLevel());
+
+    auto const &flags = eb_fact.getMultiEBCellFlagFab();
+    auto const &edge_centroid = eb_fact.getEdgeCent();
+    for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi){
+        amrex::Box const &box = mfi.validbox();
+        amrex::FabType fab_type = flags[mfi].getType(box);
+        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim){
+            auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi);
+            if (fab_type == amrex::FabType::regular) {
+                // every cell in box is all regular
+                amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_lengths_dim).ixType()),
+                                 [=](int i, int j, int k) {
+                    edge_lengths_dim(i, j, k) = 1.;
+                });
+            } else if (fab_type == amrex::FabType::covered) {
+                // every cell in box is all covered
+                amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_lengths_dim).ixType()),
+                                 [=](int i, int j, int k) {
+                    edge_lengths_dim(i, j, k) = 1.;
+                });
+            } else {
+                auto const &edge_cent = edge_centroid[idim]->const_array(mfi);
+                amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_cent).ixType()),
+                                [=](int i, int j, int k) {
+                    if (edge_cent(i, j, k) == amrex::Real(-1.0)) {
+                        // This edge is all covered
+                        edge_lengths_dim(i, j, k) = 0.;
+                    } else if (edge_cent(i, j, k) == amrex::Real(1.0)) {
+                        // This edge is all open
+                        edge_lengths_dim(i, j, k) = 1.;
+                    } else {
+                        // This edge is cut.
+                        edge_lengths_dim(i, j, k) = 1 - abs(amrex::Real(2.0)*edge_cent(i, j, k));
+                    }
+                });
+            }
+        }
+    }
+#endif
+}
+
+/**
+ * \brief Compute the ara of the mesh faces. Here the area is a value in [0, 1].
+ *        An edge of area 0 is fully covered.
+ */
+void
+WarpX::ComputeFaceAreas () {
+#ifdef AMREX_USE_EB
+    BL_PROFILE("ComputeFaceAreas");
+
+    auto const eb_fact = fieldEBFactory(maxLevel());
+    auto const &flags = eb_fact.getMultiEBCellFlagFab();
+    auto const &area_frac = eb_fact.getAreaFrac();
+
+    for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) {
+        amrex::Box const &box = mfi.validbox();
+        amrex::FabType fab_type = flags[mfi].getType(box);
+        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+            auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi);
+            if (fab_type == amrex::FabType::regular) {
+                // every cell in box is all regular
+                amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face_areas_dim).ixType()),
+                                [=](int i, int j, int k) {
+                    face_areas_dim(i, j, k) = amrex::Real(1.);
+                });
+            } else if (fab_type == amrex::FabType::covered) {
+                // every cell in box is all covered
+                amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face_areas_dim).ixType()),
+                                 [=](int i, int j, int k) {
+                    face_areas_dim(i, j, k) = amrex::Real(0.);
+                });
+            } else {
+                auto const &face = area_frac[idim]->const_array(mfi);
+                amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face).ixType()),
+                                [=](int i, int j, int k) {
+                    face_areas_dim(i, j, k) = face(i, j, k);
+                });
+            }
+        }
+    }
+#endif
+}
+
+/**
+ * \brief Scale the edges lengths by the mesh width to obtain the real lengths.
+ */
+void
+WarpX::ScaleEdges () {
+#ifdef AMREX_USE_EB
+    BL_PROFILE("ScaleEdges");
+
+    auto const &cell_size = CellSize(maxLevel());
+    auto const eb_fact = fieldEBFactory(maxLevel());
+    auto const &flags = eb_fact.getMultiEBCellFlagFab();
+
+    for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) {
+        amrex::Box const &box = mfi.validbox();
+        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+            auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi);
+            amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_lengths_dim).ixType()),
+                             [=](int i, int j, int k) {
+                edge_lengths_dim(i, j, k) *= cell_size[idim];
+            });
+        }
+    }
+#endif
+}
+
+/**
+ * \brief Scale the edges areas by the mesh width to obtain the real areas.
+ */
+void
+WarpX::ScaleAreas() {
+#ifdef AMREX_USE_EB
+    BL_PROFILE("ScaleAreas");
+
+    auto const& cell_size = CellSize(maxLevel());
+    amrex::Real full_area;
+
+    auto const eb_fact = fieldEBFactory(maxLevel());
+    auto const &flags = eb_fact.getMultiEBCellFlagFab();
+
+    for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) {
+        amrex::Box const &box = mfi.validbox();
+        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+            if (idim == 0) {
+                full_area = cell_size[1]*cell_size[2];
+            } else if (idim == 1) {
+                full_area = cell_size[0]*cell_size[2];
+            } else {
+                full_area = cell_size[0]*cell_size[1];
+            }
+            auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi);
+            amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face_areas_dim).ixType()),
+                             [=](int i, int j, int k) {
+                                face_areas_dim(i, j, k) *= full_area;
+            });
+        }
+    }
 #endif
 }

Source/Evolve/WarpXEvolve.cpp

@@ -44,6 +44,7 @@ WarpX::Evolve (int numsteps)
     bool early_params_checked = false; // check typos in inputs after step 1
 
     Real walltime, walltime_start = amrex::second();
+
     for (int step = istep[0]; step < numsteps_max && cur_time < stop_time; ++step)
     {
         Real walltime_beg_step = amrex::second();