aboutsummaryrefslogtreecommitdiff
path: root/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
blob: 75c82319c1154d6a267fba42a0a37f7d6b310e47 (plain) (blame)
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
/* 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/PsatdAlgorithm.H"
#include "SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H"
#include "SpectralKSpace.H"
#include "SpectralSolver.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 bool do_multi_J,
                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,
                                          do_multi_J, 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 grids
    {
        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 regulard grids
    {
        // 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);
        }
        else // PSATD algorithms: standard, Galilean, averaged Galilean, multi-J
        {
            if (do_multi_J)
            {
                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);
            }
            else // standard, Galilean, averaged Galilean
            {
                algorithm = std::make_unique<PsatdAlgorithm>(
                    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);
            }
        }
    }

    // - 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