1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
|
/* Copyright 2019 Remi Lehe
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#include "FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H"
#include "FieldSolver/SpectralSolver/SpectralFieldData.H"
#include "SpectralAlgorithms/PsatdAlgorithmComoving.H"
#include "SpectralAlgorithms/PsatdAlgorithmPml.H"
#include "SpectralAlgorithms/PsatdAlgorithmFirstOrder.H"
#include "SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H"
#include "SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H"
#include "SpectralKSpace.H"
#include "SpectralSolver.H"
#include "Utils/TextMsg.H"
#include "Utils/WarpXProfilerWrapper.H"
#include <memory>
#if WARPX_USE_PSATD
SpectralSolver::SpectralSolver(
const int lev,
const amrex::BoxArray& realspace_ba,
const amrex::DistributionMapping& dm,
const int norder_x, const int norder_y,
const int norder_z, const bool nodal,
const amrex::Vector<amrex::Real>& v_galilean,
const amrex::Vector<amrex::Real>& v_comoving,
const amrex::RealVect dx, const amrex::Real dt,
const bool pml, const bool periodic_single_box,
const bool update_with_rho,
const bool fft_do_time_averaging,
const int psatd_solution_type,
const int J_in_time,
const int rho_in_time,
const bool dive_cleaning,
const bool divb_cleaning)
{
// Initialize all structures using the same distribution mapping dm
// - Initialize k space object (Contains info about the size of
// the spectral space corresponding to each box in `realspace_ba`,
// as well as the value of the corresponding k coordinates)
const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx);
m_spectral_index = SpectralFieldIndex(
update_with_rho, fft_do_time_averaging, J_in_time, rho_in_time,
dive_cleaning, divb_cleaning, pml);
// - Select the algorithm depending on the input parameters
// Initialize the corresponding coefficients over k space
if (pml) // PSATD equations in the PML region
{
algorithm = std::make_unique<PsatdAlgorithmPml>(
k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,
dt, dive_cleaning, divb_cleaning);
}
else // PSATD equations in the regular domain
{
// Comoving PSATD algorithm
if (v_comoving[0] != 0. || v_comoving[1] != 0. || v_comoving[2] != 0.)
{
algorithm = std::make_unique<PsatdAlgorithmComoving>(
k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,
v_comoving, dt, update_with_rho);
}
// Galilean PSATD algorithm (only J constant in time)
else if (v_galilean[0] != 0. || v_galilean[1] != 0. || v_galilean[2] != 0.)
{
algorithm = std::make_unique<PsatdAlgorithmJConstantInTime>(
k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,
v_galilean, dt, update_with_rho, fft_do_time_averaging,
dive_cleaning, divb_cleaning);
}
else if (psatd_solution_type == PSATDSolutionType::FirstOrder)
{
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
!fft_do_time_averaging,
"psatd.do_time_averaging=1 not supported when psatd.solution_type=first-order");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
(!dive_cleaning && !divb_cleaning) || (dive_cleaning && divb_cleaning),
"warpx.do_dive_cleaning and warpx.do_divb_cleaning must be equal when psatd.solution_type=first-order");
const bool div_cleaning = (dive_cleaning && divb_cleaning);
algorithm = std::make_unique<PsatdAlgorithmFirstOrder>(
k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,
dt, div_cleaning, J_in_time, rho_in_time);
}
else if (psatd_solution_type == PSATDSolutionType::SecondOrder)
{
if (J_in_time == JInTime::Constant)
{
algorithm = std::make_unique<PsatdAlgorithmJConstantInTime>(
k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,
v_galilean, dt, update_with_rho, fft_do_time_averaging,
dive_cleaning, divb_cleaning);
}
else if (J_in_time == JInTime::Linear)
{
algorithm = std::make_unique<PsatdAlgorithmJLinearInTime>(
k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,
dt, fft_do_time_averaging, dive_cleaning, divb_cleaning);
}
}
}
// - Initialize arrays for fields in spectral space + FFT plans
field_data = SpectralFieldData(lev, realspace_ba, k_space, dm,
m_spectral_index.n_fields, periodic_single_box);
}
void
SpectralSolver::ForwardTransform (const int lev,
const amrex::MultiFab& mf,
const int field_index,
const int i_comp)
{
WARPX_PROFILE("SpectralSolver::ForwardTransform");
field_data.ForwardTransform(lev, mf, field_index, i_comp);
}
void
SpectralSolver::BackwardTransform( const int lev,
amrex::MultiFab& mf,
const int field_index,
const amrex::IntVect& fill_guards,
const int i_comp )
{
WARPX_PROFILE("SpectralSolver::BackwardTransform");
field_data.BackwardTransform(lev, mf, field_index, fill_guards, i_comp);
}
void
SpectralSolver::pushSpectralFields(){
WARPX_PROFILE("SpectralSolver::pushSpectralFields");
// Virtual function: the actual function used here depends
// on the sub-class of `SpectralBaseAlgorithm` that was
// initialized in the constructor of `SpectralSolver`
algorithm->pushSpectralFields( field_data );
}
#endif // WARPX_USE_PSATD
|
