1D3V Cartesian Support (#2307)

* Build System: Add 1D Geometry

* test PR

* test PR

* 1D cartesian yee algorithm

* fixed typo

* Fixes for PML

* 1D support related multiple changes

* Fix compilation

* change 1D to 1D_Z

* 1D Field Gather and typo fix

* 1D Charge Deposition

* Particle Pusher

* multiple changes related to 1D

* 1D diagnostics and initialization

* PlasmaInjector and PEC fixes for 1D

* clean-up delete diags file

* mobility 1D laser particle and bilinear filter

* deleted diags files

* update laser particle weight formula

* delete diags files

* Azure: Add 1D Cartesian Runner

* 1D fixes for FieldProbe

* Update Docs/source/developers/dimensionality.rst

Co-authored-by: Remi Lehe <[email protected]>

* 1d laser injection and langmuir test input files

* 1d tests

* clean up : delete print statements

* analyse simulation result for laser injection and Langmuir tests

* EOL

* delete input files for which there are no automated tests

* delete input files for which there are no automated tests

* add ignore_unused to remove warnings

* remove space

* Fix compilation issues

* fix error : macro name must be an identifier

* Small bug fix

* cleanup Python script for analysis

* bug fix

* bug fix

* Update ParticleProbe: Check 1D in-domain

* Update Source/Make.WarpX

* Update .azure-pipelines.yml

* Add USE_OPENPMD=FALSE to .azure-pipeline.yml

* resolve conflict

* resolve conflict

* fix typo

* Correct out-of-bound access

* Fix Particle BC in WarpXParticleContainer and correct path to checksumAPI in python analysis scripts

* EOL

* Fix bug : accessing out of bound index of cell in 1D

* remove 1d test for cartesian3d

* Fix CI check

* Slight style change

* Address review comments

* Fix GPU compilation Filter.cpp

* Fix CI

* Fix Indentation

* Address review comments

* More consistent ifdef for dimension bigger than 1

* Update Examples/Tests/Langmuir/analysis_langmuir_multi_1d.py

Co-authored-by: Axel Huebl <[email protected]>

* Update GNUmakefile

Co-authored-by: Axel Huebl <[email protected]>

* Update Regression/prepare_file_ci.py

Co-authored-by: Axel Huebl <[email protected]>

* Update Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H

Co-authored-by: Axel Huebl <[email protected]>

* Update Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H

Co-authored-by: Axel Huebl <[email protected]>

* Update Source/Filter/Filter.cpp

Co-authored-by: Axel Huebl <[email protected]>

* Update Source/Filter/Filter.cpp

Co-authored-by: Axel Huebl <[email protected]>

* Update Source/Filter/Filter.cpp

Co-authored-by: Axel Huebl <[email protected]>

* Update Source/Filter/Filter.cpp

Co-authored-by: Axel Huebl <[email protected]>

* Update Source/Initialization/PlasmaInjector.cpp

Co-authored-by: Axel Huebl <[email protected]>

* Update Source/Initialization/PlasmaInjector.cpp

Co-authored-by: Axel Huebl <[email protected]>

* add comment inline to explain twice push_back

* Add amrex::Abort for NCIGodfreyFilter

Co-authored-by: Axel Huebl <[email protected]>
Co-authored-by: Prabhat Kumar <[email protected]>
Co-authored-by: Remi Lehe <[email protected]>
Co-authored-by: Remi Lehe <[email protected]>

Author: 4 people

.azure-pipelines.yml

@@ -20,6 +20,10 @@ jobs:
 
   strategy:
     matrix:
+      cartesian1d:
+        WARPX_CI_REGULAR_CARTESIAN_1D: 'TRUE'
+        WARPX_CI_PSATD: 'FALSE'
+        WARPX_CI_OPENPMD: 'FALSE'
       cartesian2d:
         WARPX_CI_REGULAR_CARTESIAN_2D: 'TRUE'
       cartesian3d:

.github/workflows/source/test_ci_matrix.sh

@@ -13,8 +13,10 @@ grep "\[" ci-tests.ini > ci_all_tests.txt
 export WARPX_CI_PSATD=TRUE
 
 # Concatenate the names of all elements in CI matrix into another test file
-WARPX_CI_REGULAR_CARTESIAN_2D=TRUE python prepare_file_ci.py
+WARPX_CI_REGULAR_CARTESIAN_1D=TRUE python prepare_file_ci.py
 grep "\[" ci-tests.ini >  ci_matrix_elements.txt
+WARPX_CI_REGULAR_CARTESIAN_2D=TRUE python prepare_file_ci.py
+grep "\[" ci-tests.ini >>  ci_matrix_elements.txt
 WARPX_CI_REGULAR_CARTESIAN_3D=TRUE python prepare_file_ci.py
 grep "\[" ci-tests.ini >>  ci_matrix_elements.txt
 WARPX_CI_PYTHON_MAIN=TRUE       python prepare_file_ci.py

CMakeLists.txt

@@ -60,8 +60,8 @@ option(WarpX_SENSEI        "SENSEI in situ diagnostics"                 OFF)
 option(WarpX_QED           "QED support (requires PICSAR)"                    ON)
 option(WarpX_QED_TABLE_GEN "QED table generation (requires PICSAR and Boost)" OFF)
 
-set(WarpX_DIMS_VALUES 2 3 RZ)
-set(WarpX_DIMS 3 CACHE STRING "Simulation dimensionality (2/3/RZ)")
+set(WarpX_DIMS_VALUES 1 2 3 RZ)
+set(WarpX_DIMS 3 CACHE STRING "Simulation dimensionality (1/2/3/RZ)")
 set_property(CACHE WarpX_DIMS PROPERTY STRINGS ${WarpX_DIMS_VALUES})
 if(NOT WarpX_DIMS IN_LIST WarpX_DIMS_VALUES)
     message(FATAL_ERROR "WarpX_DIMS (${WarpX_DIMS}) must be one of ${WarpX_DIMS_VALUES}")
@@ -263,6 +263,8 @@ if(WarpX_DIMS STREQUAL 3)
     target_compile_definitions(WarpX PUBLIC WARPX_DIM_3D)
 elseif(WarpX_DIMS STREQUAL 2)
     target_compile_definitions(WarpX PUBLIC WARPX_DIM_XZ)
+elseif(WarpX_DIMS STREQUAL 1)
+    target_compile_definitions(WarpX PUBLIC WARPX_DIM_1D_Z)
 elseif(WarpX_DIMS STREQUAL RZ)
     target_compile_definitions(WarpX PUBLIC WARPX_DIM_RZ)
 endif()
@@ -331,12 +333,10 @@ install(TARGETS ${WarpX_INSTALL_TARGET_NAMES}
 # simplified library alias
 # this is currently expected by Python bindings
 if(WarpX_LIB)
-    if(WarpX_DIMS STREQUAL 3)
-        set(lib_suffix "3d")
-    elseif(WarpX_DIMS STREQUAL 2)
-        set(lib_suffix "2d")
-    elseif(WarpX_DIMS STREQUAL RZ)
+    if(WarpX_DIMS STREQUAL RZ)
         set(lib_suffix "rz")
+    else()
+        set(lib_suffix "${WarpX_DIMS}d")
     endif()
     if(WIN32)
         set(mod_ext "dll")

Docs/source/developers/dimensionality.rst

@@ -13,6 +13,7 @@ Dimensions  CMake Option
 ==========  ===================
 **3D3V**    ``WarpX_DIMS=3``
 **2D3V**    ``WarpX_DIMS=2``
+**1D3V**    ``WarpX_DIMS=1``
 **RZ**      ``WarpX_DIMS=RZ``
 ==========  ===================
 
@@ -23,35 +24,37 @@ Defines
 
 Depending on the build variant of WarpX, the following preprocessor macros will be set:
 
-==================  ===========  ===========  ===========
-Macro               3D3V         2D3V         RZ
-==================  ===========  ===========  ===========
-``AMREX_SPACEDIM``  ``3``        ``2``        ``2``
-``WARPX_DIM_3D``    **defined**  *undefined*  *undefined*
-``WARPX_DIM_XZ``    *undefined*  **defined**  *undefined*
-``WARPX_DIM_RZ``    *undefined*  *undefined*  **defined**
-==================  ===========  ===========  ===========
+==================  ===========  ===========  ===========  ===========
+Macro               3D3V         2D3V         1D3V         RZ
+==================  ===========  ===========  ===========  ===========
+``AMREX_SPACEDIM``  ``3``        ``2``        ``1``        ``2``
+``WARPX_DIM_3D``    **defined**  *undefined*  *undefined*  *undefined*
+``WARPX_DIM_1D_Z``  *undefined*  *undefined*  **defined**  *undefined*
+``WARPX_DIM_XZ``    *undefined*  **defined**  *undefined*  *undefined*
+``WARPX_DIM_RZ``    *undefined*  *undefined*  *undefined*  **defined**
+==================  ===========  ===========  ===========  ===========
 
 At the same time, the following conventions will apply:
 
-====================  ===========  ===========  ===========
-**Convention**        **3D3V**     **2D3V**     **RZ**
---------------------  -----------  -----------  -----------
+====================  ===========  ===========  ===========  ===========
+**Convention**        **3D3V**     **2D3V**     **1D3V**     **RZ**
+--------------------  -----------  -----------  -----------  -----------
 *Fields*
------------------------------------------------------------
-AMReX Box dimensions  ``3``         ``2``       ``2``
-WarpX axis labels     ``x, y, z``   ``x, z``    ``x, z``
---------------------  -----------  -----------  -----------
+------------------------------------------------------------------------
+AMReX Box dimensions  ``3``         ``2``       ``1``        ``2``
+WarpX axis labels     ``x, y, z``   ``x, z``    ``z``        ``x, z``
+--------------------  -----------  -----------  -----------  -----------
 *Particles*
------------------------------------------------------------
-AMReX AoS ``.pos()``  ``0, 1, 2``  ``0, 1``     ``0, 1``
-WarpX position names  ``x, y, z``  ``x, z``     ``r, z``
-extra SoA attribute                             ``theta``
-====================  ===========  ===========  ===========
+------------------------------------------------------------------------
+AMReX AoS ``.pos()``  ``0, 1, 2``  ``0, 1``     ``0``        ``0, 1``
+WarpX position names  ``x, y, z``  ``x, z``     ``z``        ``r, z``
+extra SoA attribute                                          ``theta``
+====================  ===========  ===========  ===========  ===========
 
 Please see the following sections for particle AoS and SoA details.
 
 Conventions
 -----------
 
 In 2D3V, we assume that the position of a particle in ``y`` is equal to ``0``.
+In 1D3V, we assume that the position of a particle in ``x`` and ``y`` is equal to ``0``.

Examples/Modules/laser_injection/analysis_1d.py

@@ -0,0 +1,175 @@
+#! /usr/bin/env python
+
+# Copyright 2021 Prabhat Kumar, Remi Lehe
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+# This file is part of the WarpX automated test suite. Its purpose is to test the
+# injection of a Gaussian laser pulse from an antenna in a 1D simulation.
+# The test calculates the envelope of each component of the laser pulse at the end of
+# the simulation and it compares it with theory. It also checks that the
+# central frequency of the Fourier transform is the expected one.
+
+import yt
+import sys
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.signal import hilbert
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+# Maximum acceptable error for this test
+relative_error_threshold = 0.05
+
+# A small number
+small_num = 1.0e-8
+
+# Physical parameters
+um = 1.e-6
+fs = 1.e-15
+c = 299792458
+
+# Parameters of the gaussian beam
+wavelength = 1.*um
+w0 = 5.*um
+tt = 10.*fs
+t_c = 24.*fs
+E_max = 4e12
+
+# laser direction
+dir_vector = np.array([0,0,1.0])
+dir_vector /= np.linalg.norm(dir_vector)
+
+
+# polarization vector
+pol_vector = np.array([1.0,1.0,0.0])
+pol_vector /= np.linalg.norm(pol_vector)
+
+# Calculates the envelope of a Gaussian beam
+def gauss_env(T,Z):
+    '''Function to compute the theory for the envelope
+    '''
+    inv_tau2 = 1./tt/tt
+    exp_arg = - inv_tau2 / c/c * (Z-T*c)*(Z-T*c)
+    return E_max * np.real(np.exp(exp_arg))
+
+# Checks envelope and central frequency for a given laser component
+def check_component(data, component, t_env_theory, coeff,Z,dz):
+    print("*** Checking " + component + " ***")
+    field = data['boxlib', component].v.squeeze()
+    env = abs(hilbert(field))
+
+    env_theory = t_env_theory*np.abs(coeff)
+
+    # Plot results
+    fig = plt.figure(figsize=(12,6))
+
+    ax1 = fig.add_subplot(221)
+    ax1.set_title('PIC field')
+    ax1.plot(Z,field)
+
+    ax2 = fig.add_subplot(222)
+    ax2.set_title('PIC envelope')
+    ax2.plot(Z,env)
+
+    ax3 = fig.add_subplot(223)
+    ax3.set_title('Theory envelope')
+    ax3.plot(Z,env_theory, label="theory")
+    ax3.plot(Z,env, label="simulation")
+    ax3.legend(loc="upper right")
+
+    ax4 = fig.add_subplot(224)
+    ax4.set_title('Difference')
+    ax4.plot(Z,env-env_theory)
+
+    plt.tight_layout()
+    plt.savefig("plt_" + component + ".png", bbox_inches='tight')
+
+    if(np.abs(coeff) < small_num):
+        is_field_zero = np.sum(np.abs(env)) < small_num
+        if is_field_zero :
+            print("[OK] Field component expected to be 0 is ~ 0")
+        else :
+            print("[FAIL] Field component expected to be 0 is NOT ~ 0")
+        assert(is_field_zero)
+        print("******\n")
+        return
+
+    fft_field = np.fft.fft(field)
+
+    freq_cols = np.fft.fftfreq(fft_field.shape[0],dz/c)
+
+    pos_max = np.unravel_index(np.abs(fft_field).argmax(), fft_field.shape)
+
+    freq = np.abs(freq_cols[pos_max[0]])
+    exp_freq = c/wavelength
+
+    relative_error_freq = np.abs(freq-exp_freq)/exp_freq
+    is_freq_ok = relative_error_freq < relative_error_threshold
+    if is_freq_ok :
+        print("[OK] Relative error frequency: {:6.3f} %".format(relative_error_freq*100))
+    else :
+        print("[FAIL] Relative error frequency: {:6.3f} %".format(relative_error_freq*100))
+    assert(is_freq_ok)
+
+    print("******\n")
+
+    relative_error_env = np.sum(np.abs(env-env_theory)) / np.sum(np.abs(env_theory))
+    is_env_ok = relative_error_env < relative_error_threshold
+    if is_env_ok :
+        print("[OK] Relative error envelope: {:6.3f} %".format(relative_error_env*100))
+    else :
+        print("[FAIL] Relative error envelope: {:6.3f} %".format(relative_error_env*100))
+    assert(is_env_ok)
+
+def check_laser(filename):
+    ds = yt.load(filename)
+
+    # yt 4.0+ has rounding issues with our domain data:
+    # RuntimeError: yt attempted to read outside the boundaries
+    # of a non-periodic domain along dimension 0.
+    if 'force_periodicity' in dir(ds): ds.force_periodicity()
+
+    z = np.linspace(
+        ds.domain_left_edge[0].v,
+        ds.domain_right_edge[0].v,
+        ds.domain_dimensions[0])
+
+    dz = (ds.domain_right_edge[0].v-ds.domain_left_edge[0].v)/(ds.domain_dimensions[0]-1)
+
+    # Compute the theory for envelope
+    env_theory = gauss_env(+t_c-ds.current_time.to_value(),z)+gauss_env(-t_c+ds.current_time.to_value(),z)
+
+    # Read laser field in PIC simulation, and compute envelope
+    all_data_level_0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+
+    b_vector = np.cross(dir_vector, pol_vector)
+
+    components = ["Ex", "Ey", "Ez", "Bx", "By", "Bz"]
+    coeffs = [
+        pol_vector[0],
+        pol_vector[1],
+        pol_vector[2],
+        b_vector[0],
+        b_vector[1],
+        b_vector[2]]
+
+    field_facts = [1, 1, 1, 1/c, 1/c, 1/c]
+
+    for comp, coeff, field_fact in zip(components, coeffs, field_facts):
+        check_component(all_data_level_0, comp, field_fact*env_theory, coeff, z, dz)
+
+def main():
+    filename_end = sys.argv[1]
+
+    check_laser(filename_end)
+
+    test_name = filename_end[:-9] # Could also be os.path.split(os.getcwd())[1]
+    checksumAPI.evaluate_checksum(test_name, filename_end)
+
+if __name__ == "__main__":
+    main()

Examples/Modules/laser_injection/inputs_1d_rt

@@ -0,0 +1,62 @@
+# Maximum number of time steps
+max_step = 240
+
+# number of grid points
+amr.n_cell =  352
+
+# Maximum allowable size of each subdomain in the problem domain;
+#    this is used to decompose the domain for parallel calculations.
+amr.max_grid_size = 32
+
+# Maximum level in hierarchy (for now must be 0, i.e., one level in total)
+amr.max_level = 0
+
+# Geometry
+geometry.coord_sys   = 0                     # 0: Cartesian
+geometry.prob_lo     = -15.e-6    # physical domain
+geometry.prob_hi     =  15.e-6
+
+boundary.field_lo = pec
+boundary.field_hi = pec
+
+warpx.serialize_ics = 1
+
+# Verbosity
+warpx.verbose = 1
+
+# Algorithms
+algo.current_deposition = esirkepov
+warpx.use_filter = 0
+
+# CFL
+warpx.cfl = 0.9
+
+# Order of particle shape factors
+algo.particle_shape = 1
+
+# Laser
+lasers.names        = laser1
+laser1.profile      = Gaussian
+laser1.position     = 0.e-6  0.e-6  0.e-6 # This point is on the laser plane
+laser1.direction    = 0.  0.  1.           # The plane normal direction
+laser1.polarization = 1. 1.  0.           # The main polarization vector
+laser1.e_max        = 4.e12                # Maximum amplitude of the laser field (in V/m)
+laser1.wavelength   = 1.0e-6               # The wavelength of the laser (in meters)
+laser1.profile_waist = 5.e-6               # The waist of the laser (in meters)
+laser1.profile_duration = 10.e-15          # The duration of the laser (in seconds)
+laser1.profile_t_peak = 24.e-15            # The time at which the laser reaches its peak (in seconds)
+laser1.profile_focal_distance = 13.109e-6  # Focal distance from the antenna (in meters)
+                                           #     With this focal distance the laser is at focus
+                                           #     at the end of the simulation.
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 20
+diag1.diag_type = Full
+
+# Moving window
+warpx.do_moving_window = 1
+warpx.moving_window_dir = z
+warpx.moving_window_v = 1.0 # in units of the speed of light
+warpx.start_moving_window_step = 20
+warpx.end_moving_window_step = 200