Open Boundaries

CMakeLists.txt

@@ -375,6 +375,12 @@ if(BUILD_TESTING)
                  WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
         )
 
+        add_test(NAME beam_in_vacuum_open_boundary.normalized.1Rank
+                 COMMAND ${HiPACE_SOURCE_DIR}/tests/beam_in_vacuum_open_boundary.normalized.1Rank.sh
+                         $<TARGET_FILE:HiPACE> ${HiPACE_SOURCE_DIR}
+                 WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
+        )
+
     endif()
 endif()
 

docs/source/run/parameters.rst

@@ -145,6 +145,17 @@ Modeling ion motion is not yet supported by the explicit solver
     The small dst is quicker for simulations with :math:`\geq 511` transverse grid points.
     The default is set accordingly.
 
+* ``fields.extended_solve`` (`bool`) optional (default `0`)
+    Extends the area of the FFT Poisson solver to the ghost cells. This can reduce artifacts
+    originating from the boundary for long simulations.
+
+* ``fields.open_boundary`` (`bool`) optional (default `0`)
+    Uses a Taylor approximation of the Greens function to solve the Poisson equations with
+    open boundary conditions. It's recommended to use this together with
+    ``fields.extended_solve = true`` and ``geometry.is_periodic = false false false``.
+    Not implemented for the explicit Helmholtz solver.
+
+
 Predictor-corrector loop parameters
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

examples/beam_in_vacuum/analysis_open_boundary.py

@@ -0,0 +1,179 @@
+#! /usr/bin/env python3
+
+# Copyright 2020-2022
+#
+# This file is part of HiPACE++.
+#
+# Authors: AlexanderSinn, MaxThevenet, Severin Diederichs
+# License: BSD-3-Clause-LBNL
+
+
+# This script compares the transverse field By with the theoretical value, plots both
+# the simulation result and the theory on the same plot, and asserts that the
+# difference is small.
+#
+# To use it, run the simulation and execute this script with
+# > ../../build/bin/hipace inputs_SI
+# > python analysis.py
+
+import numpy as np
+import scipy.constants as scc
+import argparse
+import sys
+from openpmd_viewer import OpenPMDTimeSeries
+import matplotlib
+import matplotlib.pyplot as plt
+
+parser = argparse.ArgumentParser(description='Script to analyze the correctness of the beam in vacuum')
+parser.add_argument('--normalized-units',
+                    dest='norm_units',
+                    action='store_true',
+                    default=False,
+                    help='Run the analysis in normalized units')
+parser.add_argument('--do-plot',
+                    dest='do_plot',
+                    action='store_true',
+                    default=False,
+                    help='Plot figures and save them to file')
+parser.add_argument('--output-dir',
+                    dest='output_dir',
+                    default='diags/hdf5',
+                    help='Path to the directory containing output files')
+args = parser.parse_args()
+
+ts = OpenPMDTimeSeries(args.output_dir)
+
+if args.norm_units:
+    c = 1.
+    jz0 = -1.
+    rho0 = -1.
+    mu_0 = 1.
+    eps_0 = 1.
+    R = 1.
+else:
+    # Density of the can beam
+    dens = 2.8239587008591567e23 # at this density, 1/kp = 10um, allowing for an easy comparison with normalized units
+    # Define array for transverse coordinate and theory for By and Bx
+    jz0 = - scc.e * scc.c * dens
+    rho0 = - scc.e * dens
+    c = scc.c
+    mu_0 = scc.mu_0
+    eps_0 = scc.epsilon_0
+    # Radius of the can beam
+    R = 10.e-6
+
+x_beam_mid = 2
+y_beam_mid = -1
+x_domain_len = 8
+y_domain_len = 8
+
+# Load HiPACE++ data for By in SI units
+Bx_sim, Bx_meta = ts.get_field(field='Bx', iteration=0, slice_across=['x','z'], slice_relative_position=[2*x_beam_mid/x_domain_len,0])
+By_sim, By_meta = ts.get_field(field='By', iteration=0, slice_across=['y','z'], slice_relative_position=[2*y_beam_mid/y_domain_len,0])
+jz_sim = ts.get_field(field='jz_beam', iteration=0, slice_across=['y','z'], slice_relative_position=[2*y_beam_mid/y_domain_len,0])[0]
+rho_sim = ts.get_field(field='rho', iteration=0, slice_across=['y','z'], slice_relative_position=[2*y_beam_mid/y_domain_len,0])[0]
+Ex_sim = ts.get_field(field='ExmBy', iteration=0, slice_across=['y','z'], slice_relative_position=[2*y_beam_mid/y_domain_len,0])[0] + c*By_sim
+Ey_sim = ts.get_field(field='EypBx', iteration=0, slice_across=['x','z'], slice_relative_position=[2*x_beam_mid/x_domain_len,0])[0] - c*Bx_sim
+y = Bx_meta.y
+x = By_meta.x
+
+By_th = mu_0 * jz0 * (x-x_beam_mid) / 2.
+By_th[abs(x-x_beam_mid)>=R] = mu_0 * jz0 * R**2/(2*(x[abs(x-x_beam_mid)>R]-x_beam_mid))
+Ex_th = rho0 / eps_0 * (x-x_beam_mid) / 2.
+Ex_th[abs(x-x_beam_mid)>=R] = rho0 / eps_0 * R**2/(2*(x[abs(x-x_beam_mid)>R]-x_beam_mid))
+
+Bx_th = -mu_0 * jz0 * (y-y_beam_mid) / 2.
+Bx_th[abs(y-y_beam_mid)>=R] = -mu_0 * jz0 * R**2/(2*(y[abs(y-y_beam_mid)>R]-y_beam_mid))
+Ey_th = rho0 / eps_0 * (y-y_beam_mid) / 2.
+Ey_th[abs(y-y_beam_mid)>=R] = rho0 / eps_0 * R**2/(2*(y[abs(y-y_beam_mid)>R]-y_beam_mid))
+
+jz_th = np.ones_like(x) * jz0
+jz_th[abs(x-x_beam_mid)>=R] = 0.
+rho_th = np.ones_like(x) * rho0
+rho_th[abs(x-x_beam_mid)>=R] = 0.
+
+# Plot simulation result and theory
+if args.do_plot:
+    matplotlib.rcParams.update({'font.size': 14})
+    plt.figure(figsize=(12,4))
+
+    if not args.norm_units:
+        plt.subplot(131)
+        plt.plot(1.e6*y, Bx_sim, '+-', label='HiPACE++')
+        plt.plot(1.e6*y, Bx_th, 'k--', label='theory')
+        plt.grid()
+        plt.legend()
+        plt.xlim(-50., 50.)
+        plt.xlabel('y (um)')
+        plt.ylabel('Bx (T)')
+
+        plt.subplot(132)
+        plt.plot(1.e6*x, By_sim, '+-', label='HiPACE++')
+        plt.plot(1.e6*x, By_th, 'k--', label='theory')
+        plt.grid()
+        plt.legend()
+        plt.xlim(-50., 50.)
+        plt.xlabel('x (um)')
+        plt.ylabel('By (T)')
+
+        plt.subplot(133)
+        plt.plot(1.e6*x, jz_sim, '+-', label='HiPACE++')
+        plt.plot(1.e6*x, jz_th, 'k--', label='theory')
+        plt.grid()
+        plt.legend()
+        plt.xlim(-50., 50.)
+        plt.xlabel('x (um)')
+        plt.ylabel('jz (A/m2)')
+    else:
+        plt.subplot(131)
+        plt.plot(y, Bx_sim, '+-', label='HiPACE++')
+        plt.plot(y, Bx_th, 'k--', label='theory')
+        plt.grid()
+        plt.legend()
+        plt.xlim(-5., 5.)
+        plt.xlabel('kp y')
+        plt.ylabel('c Bx / E0')
+
+        plt.subplot(132)
+        plt.plot(x, By_sim, '+-', label='HiPACE++')
+        plt.plot(x, By_th, 'k--', label='theory')
+        plt.grid()
+        plt.legend()
+        plt.xlim(-5., 5.)
+        plt.xlabel('kp x')
+        plt.ylabel('c By / E0')
+
+        plt.subplot(133)
+        plt.plot(x, jz_sim, '+-', label='HiPACE++')
+        plt.plot(x, jz_th, 'k--', label='theory')
+        plt.grid()
+        plt.legend()
+        plt.xlim(-5., 5.)
+        plt.xlabel('kp x')
+        plt.ylabel('jz /IA')
+
+    plt.tight_layout()
+
+    plt.savefig("beam_in_vacuum.png", bbox_inches="tight")
+
+# Assert that the simulation result is close enough to theory
+error_jz = np.sum((jz_sim-jz_th)**2) / np.sum((jz_th)**2)
+print("total relative error jz: " + str(error_jz) + " (tolerance = 0.1)")
+
+error_Bx = np.sum((Bx_sim-Bx_th)**2) / np.sum((Bx_th)**2)
+print("total relative error Bx: " + str(error_Bx) + " (tolerance = 0.005)")
+
+error_By = np.sum((By_sim-By_th)**2) / np.sum((By_th)**2)
+print("total relative error By: " + str(error_By) + " (tolerance = 0.015)")
+
+error_Ex = np.sum((Ex_sim-Ex_th)**2) / np.sum((Ex_th)**2)
+print("total relative error Ex: " + str(error_Ex) + " (tolerance = 0.015)")
+
+error_Ey = np.sum((Ey_sim-Ey_th)**2) / np.sum((Ey_th)**2)
+print("total relative error Ey: " + str(error_Ey) + " (tolerance = 0.005)")
+
+assert(error_jz < .1)
+assert(error_Bx < .005)
+assert(error_By < .015)
+assert(error_Ex < .015)
+assert(error_Ey < .005)

src/Hipace.cpp

@@ -573,7 +573,7 @@ Hipace::SolveOneSlice (int islice_coarse, const int ibox,
                                  ijz == ijx+4 && ijz_beam == ijx+5 && irho == ijx+6 );
             amrex::MultiFab j_slice(m_fields.getSlices(lev, WhichSlice::This),
                                     amrex::make_alias, Comps[WhichSlice::This]["jx"], 7);
-            j_slice.FillBoundary(Geom(lev).periodicity());
+            if (!m_fields.m_extended_solve) j_slice.FillBoundary(Geom(lev).periodicity());
 
             m_fields.SolvePoissonExmByAndEypBx(Geom(), m_comm_xy, lev, islice);
 
@@ -583,7 +583,7 @@ Hipace::SolveOneSlice (int islice_coarse, const int ibox,
                                              WhichSlice::This);
             m_fields.AddBeamCurrents(lev, WhichSlice::This);
 
-            j_slice.FillBoundary(Geom(lev).periodicity());
+            if (!m_fields.m_extended_solve) j_slice.FillBoundary(Geom(lev).periodicity());
 
             m_fields.SolvePoissonEz(Geom(), lev, islice);
             m_fields.SolvePoissonBz(Geom(), lev, islice);
@@ -912,10 +912,13 @@ Hipace::PredictorCorrectorLoopToSolveBxBy (const int islice_local, const int lev
 
     /* Guess Bx and By */
     m_fields.InitialBfieldGuess(relative_Bfield_error, m_predcorr_B_error_tolerance, lev);
-    amrex::ParallelContext::push(m_comm_xy);
-     // exchange ExmBy EypBx Ez Bx By Bz
-    m_fields.getSlices(lev, WhichSlice::This).FillBoundary(Geom(lev).periodicity());
-    amrex::ParallelContext::pop();
+
+    if (!m_fields.m_extended_solve) {
+        amrex::ParallelContext::push(m_comm_xy);
+        // exchange ExmBy EypBx Ez Bx By Bz
+        m_fields.getSlices(lev, WhichSlice::This).FillBoundary(Geom(lev).periodicity());
+        amrex::ParallelContext::pop();
+    }
 
     /* creating temporary Bx and By arrays for the current and previous iteration */
     amrex::MultiFab Bx_iter(m_fields.getSlices(lev, WhichSlice::This).boxArray(),
@@ -976,12 +979,14 @@ Hipace::PredictorCorrectorLoopToSolveBxBy (const int islice_local, const int lev
                                          ibox, m_do_beam_jx_jy_deposition, WhichSlice::Next);
         m_fields.AddBeamCurrents(lev, WhichSlice::Next);
 
-        amrex::ParallelContext::push(m_comm_xy);
-        // need to exchange jx jy jx_beam jy_beam
-        amrex::MultiFab j_slice_next(m_fields.getSlices(lev, WhichSlice::Next),
-                                     amrex::make_alias, Comps[WhichSlice::Next]["jx"], 4);
-        j_slice_next.FillBoundary(Geom(lev).periodicity());
-        amrex::ParallelContext::pop();
+        if (!m_fields.m_extended_solve) {
+            amrex::ParallelContext::push(m_comm_xy);
+            // need to exchange jx jy jx_beam jy_beam
+            amrex::MultiFab j_slice_next(m_fields.getSlices(lev, WhichSlice::Next),
+                                         amrex::make_alias, Comps[WhichSlice::Next]["jx"], 4);
+            j_slice_next.FillBoundary(Geom(lev).periodicity());
+            amrex::ParallelContext::pop();
+        }
 
         /* Calculate Bx and By */
         m_fields.SolvePoissonBx(Bx_iter, Geom(), lev, islice);
@@ -1010,10 +1015,12 @@ Hipace::PredictorCorrectorLoopToSolveBxBy (const int islice_local, const int lev
         jx_beam_next.setVal(0., m_fields.m_slices_nguards);
         jy_beam_next.setVal(0., m_fields.m_slices_nguards);
 
-        amrex::ParallelContext::push(m_comm_xy);
-         // exchange Bx By
-        m_fields.getSlices(lev, WhichSlice::This).FillBoundary(Geom(lev).periodicity());
-        amrex::ParallelContext::pop();
+        if (!m_fields.m_extended_solve) {
+            amrex::ParallelContext::push(m_comm_xy);
+            // exchange Bx By
+            m_fields.getSlices(lev, WhichSlice::This).FillBoundary(Geom(lev).periodicity());
+            amrex::ParallelContext::pop();
+        }
 
         /* Update force terms using the calculated Bx and By */
         m_multi_plasma.AdvanceParticles(m_fields, m_laser, geom[lev],

src/fields/Fields.H

@@ -307,6 +307,9 @@ public:
     static amrex::IntVect m_slices_nguards;
     /** Number of guard cells for poisson solver MultiFab */
     static amrex::IntVect m_poisson_nguards;
+    /** If the poisson solver should include the guard cells */
+    bool m_extended_solve = false;
+
 private:
     /** Vector over levels, array of 4 slices required to compute current slice */
     amrex::Vector<std::array<amrex::MultiFab, m_nslices>> m_slices;
@@ -316,6 +319,12 @@ private:
     amrex::Vector<amrex::FArrayBox> m_tmp_densities;
     /** Stores temporary values for z interpolation in Fields::Copy */
     amrex::Gpu::DeviceVector<amrex::Real> m_rel_z_vec;
+    /** Number of guard cells where ExmBy and EypBx are calculated */
+    amrex::IntVect m_exmby_eypbx_nguard;
+    /** Number of guard cells where sources for poisson equation are included */
+    amrex::IntVect m_source_nguard;
+    /** If lev_0 should be solved with open boundary conditions */
+    bool m_open_boundary = false;
 };
 
 #endif