ECT conformal solver (#1923)

* adding base implementation of the conformal solver

* adding some preprocessor directives

* qualifying the isnan's correctly

* getting rid of some unused fields

* computing S_stab on the fly

* using empty in length check

* Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp

* replaced a looponcpu with a parallelfor

* trying to not pass lending/borrowing info by reference in evolveB

* passing data using dataPtr

* converting inds into a device vector

* simplifying some +=

* extracting the inds BaseFabs as DataPtr

* Revert "extracting the inds BaseFabs as DataPtr"

This reverts commit bc709d2fcaa347c119514de79a3f57169a05af65.

* implementing new way ogf handling cell extensions (gpu compatible)

* fixing some white spaces

* removed a typo

* made a function static

* tidying up

* tidying up

* getting rid of the rho vector

* moving the rho update to evolveE

* refectoring the cell extension

* small chenge in input parameters

* updating WarpX.H

* bug fix

* fixing another bug

* rotating the cube

* updated gitignore

* using the conformal soler in the cube test

* fixing the signature of a function

* trying to fix the setVal problem

* Revert "trying to fix the setVal problem"

This reverts commit c7d0e5e3709b730ff570183b2a6df5f587ca4640.

* trying to fix the setVal problem

* cleaning some comments

* removing an abort in device code

* Including geometric information in the external field initializer

* cleaned up the test

* moving the geometric data initialization before the external fields initialization

* changing way of saving info about intruded cells

* fixed some white spaces

* Using S_mod instead of S_red and S_enl

* substituting a std::vector with amrex::Array1D

* bug fix in the uint8_t -> coords maps

* Condensed ect cell info into one single MultiFab

* fixing some includes

* fixing some more includes

* fixed a typo

* making some functions available in gpu code

* using tilebox instead of validbox for the geometry initialization

* communicating the geometric info in the guard cells

* fixing the guard cells initialization for ect

* fixing an unused parameter

* fixing an unused parameter

* ignored some unused variables when EB is not supported

* cleaning up

* cleaning up

* ignored some unused variables when EB is not supported

* evolving rho just for ect

* handling some more unused variables

* removing some white spaces

* adding the rotated cube test

* removed some white spaces

* change cells into faces

* small fix in unused parameters

* fixed a comment

* for safety for now just initialize the geometric data when lev==maxLevel

* adding some preprocessor directives to isolate EB code

* updating Makefile stuff

* Improving the rotated cube test

* bug fix in mesh refinement

* fixing boundary conditions in rotated cube test

* deactivating MPI in rotated cube test

* activating mpi in test

* improving outputs of cells extension

* improved the documentation

* Update Source/EmbeddedBoundary/WarpXInitEB.cpp

Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>

* Update Source/EmbeddedBoundary/WarpXInitEB.cpp

Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>

* Removed some obsolete isnan's

* undone change to gitignore

* adding some missing AMREX_GPU_DEVICEs

* Changing some amrex::Real(0) into 0.

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Making some variables const

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Adding amrex:: to some ignore_unused

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Enhancing readability of some parallelfor's

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Removed some unecessary includes

* Fixing some tags

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Fixing the function CountExtFaces

* added a comment for Rhofield

* Fixed a typo in CountExtCells

* trying to remove accesses to private members in ComputeEdgeLengths

* making some functions public

* undoing some useless changes

* undoing some useless changes

* adding some AMREX_GPU_DEVICEs

* adding some unused variables

* adding some AMREX_GPU_DEVICEs to fix warnings

* fixing relative error

* Making several variables const

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Fixed a comment

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* changing some double to amrex::Real

* removing commented lines

* changing some double to amrex::Real

* removing commented lines

* inizialing nelems_x,y,z to avoid a warning

* Fixing number of cells

* Adding missing analysis routine in rotated cube test

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Making some variables constexpr

Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>

* Improving some reduce operations

Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>

* reading time directly from the output

* fixed a few data types

* fixing another int

* Replacing or->||, and->&&, not->!

* Adding license info

* Adding a white space

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* removed some unused imports

* Moving the doxygen documentation

* Including AMReX_LayoutData.H

* removing some useless parameters

* adding ect solver to the doc

* Removing some useless reductions

* Small change for consistency

* Handling the resizing of arrays

* Handling correctly the resing of arrays

* Fixing some whitespaces

* Fixing an indentation

* Improving a comment

* Removed a useless initialization

* Renamed Rhofield to ECTRhofield

* Renamed Rhofield to ECTRhofield

* Added some if conditions to isolate some EB related code

* Improved some comments

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

* Made some functions not member of WarpX anymore

* Isolated some EB-only code

* Isolated some EB-only code

* Using the square brackets to access vectors

* Ignoring some unused variables

* Bug fix in initialization

* Removed a redundant if

* Using enumeration for connectivity and improved names

* Added a comment

* added a few comments

* Removed some useless template parameters

* Printing some info only in verbose mode

* Revert "Removed some useless template parameters"

This reverts commit 2c527980e6872c0212767c27f00e2b53ecbcfd0a.

* Fixed a bug with the neighbours enumeration

* Initializing geom data alsoo in the ghost cells

Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>
Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

3 files changed, 195 insertions, 6 deletions

diff --git a/Examples/Modules/embedded_boundary_cube/analysis_fields.py b/Examples/Modules/embedded_boundary_cube/analysis_fields.py
index 58dad1bb7..e4b461edd 100755
--- a/Examples/Modules/embedded_boundary_cube/analysis_fields.py
+++ b/Examples/Modules/embedded_boundary_cube/analysis_fields.py
@@ -32,7 +32,13 @@ Lx = 1
 Ly = 1
 Lz = 1
 h_2 = (m * pi / Lx) ** 2 + (n * pi / Ly) ** 2 + (p * pi / Lz) ** 2
-t = 1.3342563807926085e-08
+
+# 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
 Bx_th = np.zeros(ncells)
@@ -65,11 +71,6 @@ for i in range(ncells[0]):
                                    (-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)
-
 rel_tol_err = 1e-1
 
 # Compute relative l^2 error on By
diff --git a/Examples/Modules/embedded_boundary_rotated_cube/analysis_fields.py b/Examples/Modules/embedded_boundary_rotated_cube/analysis_fields.py
new file mode 100755
index 000000000..7765c5d69
--- /dev/null
+++ b/Examples/Modules/embedded_boundary_rotated_cube/analysis_fields.py
@@ -0,0 +1,116 @@
+#! /usr/bin/env python
+
+# Copyright 2021 Lorenzo Giacomel
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+
+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 rotated by pi/8.
+# 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]
+m = 0
+n = 1
+p = 1
+Lx = 1.06
+Ly = 1.06
+Lz = 1.06
+h_2 = (m * pi / Lx) ** 2 + (n * pi / Ly) ** 2 + (p * pi / Lz) ** 2
+theta = np.pi/8
+
+# 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()
+
+rel_tol_err = 1e-2
+my_grid = ds.index.grids[0]
+
+By_sim = my_grid['raw', 'By_fp'].squeeze().v
+Bz_sim = my_grid['raw', 'Bz_fp'].squeeze().v
+
+ncells = np.array(np.shape(By_sim[:, :, :, 0]))
+dx = (hi[0] - lo[0])/ncells[0]
+dy = (hi[1] - lo[1])/ncells[1]
+dz = (hi[2] - lo[2])/ncells[2]
+
+# 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]):
+            x0 = (i+0.5)*dx + lo[0]
+            y0 = j*dy + lo[1]
+            z0 = (k+0.5)*dz + lo[2]
+
+            x = x0
+            y = y0*np.cos(-theta)-z0*np.sin(-theta)
+            z = y0*np.sin(-theta)+z0*np.cos(-theta)
+            By = -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)) *
+                                                        np.cos(np.sqrt(2) *
+                                                        np.pi / Lx * c * t))
+
+            Bz = mu_0*(np.cos(m * pi/Lx * (x - Lx/2)) *
+                       np.cos(n * pi/Ly * (y - Ly/2)) *
+                       np.sin(p * pi/Lz * (z - Lz/2)) *
+                       np.cos(np.sqrt(2) * np.pi / Lx * c * t))
+
+            By_th[i, j, k] = (By*np.cos(theta) - Bz*np.sin(theta))*(By_sim[i, j, k, 0] != 0)
+
+            x0 = (i+0.5)*dx + lo[0]
+            y0 = (j+0.5)*dy + lo[1]
+            z0 = k*dz + lo[2]
+
+            x = x0
+            y = y0*np.cos(-theta)-z0*np.sin(-theta)
+            z = y0*np.sin(-theta)+z0*np.cos(-theta)
+
+            By = -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)) *
+                                                        np.cos(np.sqrt(2) *
+                                                               np.pi / Lx * c * t))
+
+            Bz = mu_0*(np.cos(m * pi/Lx * (x - Lx/2)) *
+                       np.cos(n * pi/Ly * (y - Ly/2)) *
+                       np.sin(p * pi/Lz * (z - Lz/2)) *
+                       np.cos(np.sqrt(2) * np.pi / Lx * c * t))
+
+            Bz_th[i, j, k] = (By*np.sin(theta) + Bz*np.cos(theta))*(Bz_sim[i, j, k, 0] != 0)
+
+
+# Compute relative l^2 error on By
+rel_err_y = np.sqrt( np.sum(np.square(By_sim[:, :, :, 0] - By_th)) / np.sum(np.square(By_th)))
+assert(rel_err_y < rel_tol_err)
+# Compute relative l^2 error on Bz
+rel_err_z = np.sqrt( np.sum(np.square(Bz_sim[:, :, :, 0] - Bz_th)) / np.sum(np.square(Bz_th)))
+assert(rel_err_z < rel_tol_err)
+
+test_name = os.path.split(os.getcwd())[1]
+
+checksumAPI.evaluate_checksum(test_name, filename)
diff --git a/Examples/Modules/embedded_boundary_rotated_cube/inputs_3d b/Examples/Modules/embedded_boundary_rotated_cube/inputs_3d
new file mode 100644
index 000000000..7a6d51cf7
--- /dev/null
+++ b/Examples/Modules/embedded_boundary_rotated_cube/inputs_3d
@@ -0,0 +1,72 @@
+stop_time = 5.303669113650618e-09
+amr.n_cell = 32 32 32
+amr.max_grid_size = 128
+amr.max_level = 0
+
+geometry.coord_sys   = 0
+geometry.prob_lo     = -0.8 -0.8 -0.8
+geometry.prob_hi     =  0.8  0.8  0.8
+warpx.cfl = 1
+
+boundary.field_lo = pec pec pec
+boundary.field_hi = pec pec pec
+
+algo.maxwell_solver = ect
+
+my_constants.xmin = -0.53
+my_constants.ymin = -0.53
+my_constants.zmin = -0.53
+my_constants.xmax = 0.53
+my_constants.ymax = 0.53
+my_constants.zmax = 0.53
+my_constants.pi = 3.141592653589793
+my_constants.theta = pi/8
+
+warpx.eb_implicit_function = "max(max(max(x+xmin,-(x+xmax)), max(y*cos(-theta)-z*sin(-theta)+ymin,-(y*cos(-theta)-z*sin(-theta)+ymax))), max(y*sin(-theta)+z*cos(-theta)+zmin,-(y*sin(-theta)+z*cos(-theta)+zmax)))"
+
+my_constants.m = 0
+my_constants.n = 1
+my_constants.p = 1
+my_constants.Lx = 1.06
+my_constants.Ly = 1.06
+my_constants.Lz = 1.06
+my_constants.x_cent = 0.
+my_constants.y_cent = 0.
+my_constants.z_cent = 0.
+my_constants.h_2 = (m * pi / Lx) ** 2 + (n * pi / Ly) ** 2 + (p * pi / Lz) ** 2
+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) = -2/h_2 * mu_0 * (m * pi / Lx) * (p * pi / Lz) *
+                                          sin(m * pi / Lx * (x - Lx / 2 - x_cent)) *
+                                          cos(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) *
+                                          cos(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent))
+
+warpx.By_external_grid_function(x,y,z) = -2/h_2 * mu_0 * (n * pi / Ly) * (p * pi / Lz) *
+                                          cos(m * pi / Lx * (x - Lx / 2 - x_cent)) *
+                                          sin(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) *
+                                          cos(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) *
+                                          cos(theta) -
+                                          mu_0 *
+                                          cos(m * pi / Lx * (x - Lx / 2 - x_cent)) *
+                                          cos(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) *
+                                          sin(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) *
+                                          sin(theta)
+
+warpx.Bz_external_grid_function(x,y,z) = mu_0 *
+                                         cos(m * pi / Lx * (x - Lx / 2 - x_cent)) *
+                                         cos(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) *
+                                         sin(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) *
+                                         cos(theta) -
+                                         2/h_2 * mu_0 * (n * pi / Ly) * (p * pi / Lz) *
+                                         cos(m * pi / Lx * (x - Lx / 2 - x_cent)) *
+                                         sin(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) *
+                                         cos(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) *
+                                         sin(theta)
+
+
+diagnostics.diags_names = diag1
+diag1.intervals = 1000
+diag1.diag_type = Full
+diag1.plot_raw_fields = 1
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz