From 51e866e9a3f450d7767df97b0a2011cf65857c7a Mon Sep 17 00:00:00 2001
From: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
Date: Thu, 27 Feb 2020 11:17:56 -0800
Subject: Implement div(E) diagnostics for spectral case (#720)

* Implement div(E) diagnostics for spectral case.

* split travis tests in bigger matrix

* split more TravisCI tests, add electrostatic, use defaults values

* typo

* Move computation of div(E) to base class SpectralBaseAlgorithm.

* need to split psatd too

* consistent variable names and use function to avoid duplication

* fix typo for qed tests

* typo

* also need to update run_tests.sg

* Update copyright tags.

* change matrix

* Add test of div(E) vs rho/epsilon_0 in PML test.

* SpectralFieldIndex: reuse memory slot for Bx when computing divE.

Co-authored-by: MaxThevenet <mthevenet@lbl.gov>
---
 .../SpectralSolver/SpectralAlgorithms/Make.package |  1 +
 .../SpectralAlgorithms/SpectralBaseAlgorithm.H     |  9 ++-
 .../SpectralAlgorithms/SpectralBaseAlgorithm.cpp   | 69 ++++++++++++++++++++++
 3 files changed, 78 insertions(+), 1 deletion(-)
 create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp

(limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms')

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
index 82763f501..69758e7a9 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
@@ -1,4 +1,5 @@
 CEXE_headers += SpectralBaseAlgorithm.H
+CEXE_sources += SpectralBaseAlgorithm.cpp
 CEXE_headers += PsatdAlgorithm.H
 CEXE_sources += PsatdAlgorithm.cpp
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
index 2487eae78..743d6d1c7 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
@@ -1,4 +1,4 @@
-/* Copyright 2019 Remi Lehe
+/* Copyright 2019 Remi Lehe, Edoardo Zoni
  *
  * This file is part of WarpX.
  *
@@ -31,6 +31,13 @@ class SpectralBaseAlgorithm
         // calls the subclass's destructor.
         virtual ~SpectralBaseAlgorithm() {};
 
+        /**
+         * \brief Compute spectral divergence of E
+         */
+        void ComputeSpectralDivE ( SpectralFieldData& field_data,
+                                   const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
+                                   amrex::MultiFab& divE );
+
     protected: // Meant to be used in the subclasses
 
         using SpectralRealCoefficients = \
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp
new file mode 100644
index 000000000..c917391f9
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp
@@ -0,0 +1,69 @@
+/* Copyright 2019 Edoardo Zoni
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include <SpectralBaseAlgorithm.H>
+#include <cmath>
+
+using namespace amrex;
+
+/**
+ * \brief Compute spectral divergence of E
+ */
+void
+SpectralBaseAlgorithm::ComputeSpectralDivE (
+    SpectralFieldData& field_data,
+    const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
+    amrex::MultiFab& divE )
+{
+    using Idx = SpectralFieldIndex;
+
+    // Forward Fourier transform of E
+    field_data.ForwardTransform( *Efield[0], Idx::Ex, 0 );
+    field_data.ForwardTransform( *Efield[1], Idx::Ey, 0 );
+    field_data.ForwardTransform( *Efield[2], Idx::Ez, 0 );
+
+    // Loop over boxes
+    for (MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
+
+        const Box& bx = field_data.fields[mfi].box();
+
+        // Extract arrays for the fields to be updated
+        Array4<Complex> fields = field_data.fields[mfi].array();
+        // Extract pointers for the k vectors
+        const Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr();
+#if (AMREX_SPACEDIM==3)
+        const Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr();
+#endif
+        const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr();
+
+        // Loop over indices within one box
+        ParallelFor(bx,
+        [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            using Idx = SpectralFieldIndex;
+            // Shortcuts for the components of E
+            const Complex Ex = fields(i,j,k,Idx::Ex);
+            const Complex Ey = fields(i,j,k,Idx::Ey);
+            const Complex Ez = fields(i,j,k,Idx::Ez);
+            // k vector values
+            const Real kx = modified_kx_arr[i];
+#if (AMREX_SPACEDIM==3)
+            const Real ky = modified_ky_arr[j];
+            const Real kz = modified_kz_arr[k];
+#else
+            constexpr Real ky = 0;
+            const Real kz = modified_kz_arr[j];
+#endif
+            const Complex I = Complex{0,1};
+
+            // div(E) in Fourier space
+            fields(i,j,k,Idx::divE) = I*(kx*Ex+ky*Ey+kz*Ez);
+        });
+    }
+
+    // Backward Fourier transform
+    field_data.BackwardTransform( divE, Idx::divE, 0 );
+};
-- 
cgit v1.2.3