2D EM solver with EB

Examples/Modules/embedded_boundary_cube/analysis_fields_2d.py

@@ -0,0 +1,60 @@
+#! /usr/bin/env python
+
+import yt
+import os
+import 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_y = \mu \cos(k_x x)\cos(k_z z)\cos( \omega_p t)$$
+# with
+# $$ k_x = \frac{m\pi}{L}$$
+# $$ k_y = \frac{n\pi}{L}$$
+# $$ k_z = \frac{p\pi}{L}$$
+
+hi = [0.8, 0.8]
+lo = [-0.8, -0.8]
+ncells = [32, 32, 1]
+dx = (hi[0] - lo[0]) / ncells[0]
+dz = (hi[1] - lo[1]) / ncells[1]
+m = 0
+n = 1
+Lx = 1
+Lz = 1
+
+# 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)
+
+t = ds.current_time.to_value()
+
+# Compute the analytic solution
+By_th = np.zeros(ncells)
+for i in range(ncells[0]):
+    for j in range(ncells[1]):
+        x = (i+0.5) * dx + lo[0]
+        z = (j+0.5) * dz + lo[1]
+
+        By_th[i, j, 0] = mu_0 * (np.cos(m * pi / Lx * (x - Lx / 2)) *
+                                 np.cos(n * pi / Lz * (z - Lz / 2)) *
+                                 (-Lx / 2 <= x < Lx / 2) *
+                                 (-Lz / 2 <= z < Lz / 2) *
+                                 np.cos(np.pi / Lx * c * t))
+
+rel_tol_err = 1e-3
+
+# Compute relative l^2 error on By
+By_sim = data['By'].to_ndarray()
+rel_err_y = np.sqrt(np.sum(np.square(By_sim - By_th)) / np.sum(np.square(By_th)))
+assert (rel_err_y < rel_tol_err)
+
+test_name = os.path.split(os.getcwd())[1]
+
+checksumAPI.evaluate_checksum(test_name, filename)

Examples/Modules/embedded_boundary_cube/inputs_2d

@@ -0,0 +1,36 @@
+stop_time = 1.3342563807926085e-08
+amr.n_cell = 32 32
+amr.max_grid_size = 128
+amr.max_level = 0
+
+geometry.coord_sys   = 0
+geometry.prob_lo     = -0.8 -0.8
+geometry.prob_hi     =  0.8  0.8
+warpx.cfl = 1
+
+boundary.field_lo = pec pec
+boundary.field_hi = pec pec
+
+my_constants.xmin = -0.5
+my_constants.zmin = -0.5
+my_constants.xmax = 0.5
+my_constants.zmax = 0.5
+# Alternatively one could use parser to build EB
+# Note that for amrex EB implicit function, >0 is covered, =0 is boundary and <0 is regular.
+warpx.eb_implicit_function = "max(max(x+xmin,-(x+xmax)), max(z+zmin,-(z+zmax)))"
+
+warpx.B_ext_grid_init_style = parse_B_ext_grid_function
+
+my_constants.m = 0
+my_constants.p = 1
+my_constants.Lx = 1
+my_constants.Lz = 1
+
+warpx.Bz_external_grid_function(x,y,z) = 0
+warpx.Bx_external_grid_function(x,y,z) = 0
+warpx.By_external_grid_function(x,y,z) = cos(m * pi / Lx * (x - Lx / 2)) * cos(p * pi / Lz * (z - Lz / 2))*mu0
+
+diagnostics.diags_names = diag1
+diag1.intervals = 1
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz

Examples/Modules/embedded_boundary_rotated_cube/analysis_fields_2d.py

@@ -0,0 +1,64 @@
+#! /usr/bin/env python
+
+import yt
+import os
+import 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_y = \mu \cos(k_x x)\cos(k_z z)\cos( \omega_p t)$$
+# with
+# $$ k_x = \frac{m\pi}{L}$$
+# $$ k_y = \frac{n\pi}{L}$$
+# $$ k_z = \frac{p\pi}{L}$$
+
+hi = [0.8, 0.8]
+lo = [-0.8, -0.8]
+ncells = [32, 32]
+dx = (hi[0] - lo[0]) / ncells[0]
+dz = (hi[1] - lo[1]) / ncells[1]
+m = 0
+n = 1
+Lx = 1.06
+Lz = 1.06
+
+# 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)
+my_grid = ds.index.grids[0]
+
+By_sim = my_grid['By'].squeeze().v
+
+t = ds.current_time.to_value()
+
+theta = np.pi/8
+
+# Compute the analytic solution
+By_th = np.zeros(ncells)
+for i in range(ncells[0]):
+    for j in range(ncells[1]):
+        x = i * dx + lo[0]
+        z = j * dz + lo[1]
+        xr = x*np.cos(-theta) + z*np.sin(-theta)
+        zr = -x*np.sin(-theta) + z*np.cos(-theta)
+
+        By_th[i, j] = mu_0 * (np.cos(m * pi / Lx * (xr - Lx / 2)) *
+                              np.cos(n * pi / Lz * (zr - Lz / 2)) *
+                              np.cos(np.pi / Lx * c * t))*(By_sim[i, j] != 0)
+
+rel_tol_err = 1e-1
+
+# Compute relative l^2 error on By
+rel_err_y = np.sqrt(np.sum(np.square(By_sim - By_th)) / np.sum(np.square(By_th)))
+assert (rel_err_y < rel_tol_err)
+
+test_name = os.path.split(os.getcwd())[1]
+
+checksumAPI.evaluate_checksum(test_name, filename)

Examples/Modules/embedded_boundary_rotated_cube/inputs_2d

@@ -0,0 +1,44 @@
+stop_time = 8.019424744948937e-09
+amr.n_cell = 32 32
+amr.max_grid_size = 128
+amr.max_level = 0
+
+geometry.coord_sys   = 0
+geometry.prob_lo     = -0.8 -0.8
+geometry.prob_hi     =  0.8  0.8
+warpx.cfl = 1
+
+boundary.field_lo = pec pec
+boundary.field_hi = pec pec
+
+algo.maxwell_solver = ect
+
+my_constants.xmin = -0.53
+my_constants.zmin = -0.53
+my_constants.xmax = 0.53
+my_constants.zmax = 0.53
+my_constants.pi = 3.141592653589793
+my_constants.theta = pi/8
+
+warpx.eb_implicit_function = "xr=x*cos(-theta)+z*sin(-theta); zr=-x*sin(-theta)+z*cos(-theta); max(max(xr+xmin,-(xr+xmax)), max(zr+zmin,-(zr+zmax)))"
+
+my_constants.m = 0
+my_constants.p = 1
+my_constants.Lx = 1.06
+my_constants.Lz = 1.06
+my_constants.x_cent = 0.
+my_constants.z_cent = 0.
+my_constants.mu_0 = 1.25663706212e-06
+
+warpx.B_ext_grid_init_style = parse_B_ext_grid_function
+
+warpx.Bx_external_grid_function(x,y,z) = 0
+warpx.By_external_grid_function(x,y,z) = mu_0 *
+                                         cos(m * pi / Lx * (x*cos(-theta)+z*sin(-theta) - Lx / 2 - x_cent)) *
+                                         cos(p * pi / Lz * (-x*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent))
+warpx.Bz_external_grid_function(x,y,z) = 0
+
+diagnostics.diags_names = diag1
+diag1.intervals = 1000
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz

Regression/Checksum/benchmarks_json/embedded_boundary_cube_2d.json

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

Regression/Checksum/benchmarks_json/embedded_boundary_rotated_cube_2d.json

@@ -0,0 +1,10 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 0.0002722100736459592,
+    "Bz": 0.0,
+    "Ex": 78275.78138124556,
+    "Ey": 0.0,
+    "Ez": 31681.635842293996
+  }
+}

Regression/WarpX-tests.ini

@@ -2310,6 +2310,22 @@ doVis = 0
 compareParticles = 0
 analysisRoutine = Examples/Modules/embedded_boundary_cube/analysis_fields.py
 
+[embedded_boundary_cube_2d]
+buildDir = .
+inputFile = Examples/Modules/embedded_boundary_cube/inputs_2d
+runtime_params =
+dim = 2
+addToCompileString = USE_EB=TRUE
+restartTest = 0
+useMPI = 1
+numprocs = 1
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Modules/embedded_boundary_cube/analysis_fields_2d.py
+
 [embedded_boundary_rotated_cube]
 buildDir = .
 inputFile = Examples/Modules/embedded_boundary_rotated_cube/inputs_3d
@@ -2326,6 +2342,22 @@ doVis = 0
 compareParticles = 0
 analysisRoutine = Examples/Modules/embedded_boundary_rotated_cube/analysis_fields.py
 
+[embedded_boundary_rotated_cube_2d]
+buildDir = .
+inputFile = Examples/Modules/embedded_boundary_rotated_cube/inputs_2d
+runtime_params =
+dim = 2
+addToCompileString = USE_EB=TRUE
+restartTest = 0
+useMPI = 1
+numprocs = 1
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Modules/embedded_boundary_rotated_cube/analysis_fields_2d.py
+
 [dirichletbc]
 buildDir = .
 inputFile = Examples/Tests/ElectrostaticDirichletBC/inputs_2d