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
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
|
/* Copyright 2023 The WarpX Community
*
* This file is part of WarpX.
*
* Authors: Roelof Groenewald (TAE Technologies)
*
* License: BSD-3-Clause-LBNL
*/
#ifndef WARPX_HYBRIDPICMODEL_H_
#define WARPX_HYBRIDPICMODEL_H_
#include "HybridPICModel_fwd.H"
#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H"
#include "Utils/Parser/ParserUtils.H"
#include "Utils/WarpXConst.H"
#include "Utils/WarpXProfilerWrapper.H"
#include "WarpX.H"
#include <AMReX_Array.H>
#include <AMReX_REAL.H>
/**
* \brief This class contains the parameters needed to evaluate hybrid field
* solutions (kinetic ions with fluid electrons).
*/
class HybridPICModel
{
public:
HybridPICModel (int nlevs_max); // constructor
/** Read user-defined model parameters. Called in constructor. */
void ReadParameters ();
/** Allocate hybrid-PIC specific multifabs. Called in constructor. */
void AllocateMFs (int nlevs_max);
void AllocateLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm,
const int ncomps, const amrex::IntVect& ngJ, const amrex::IntVect& ngRho,
const amrex::IntVect& jx_nodal_flag, const amrex::IntVect& jy_nodal_flag,
const amrex::IntVect& jz_nodal_flag, const amrex::IntVect& rho_nodal_flag);
/** Helper function to clear values from hybrid-PIC specific multifabs. */
void ClearLevel (int lev);
void InitData ();
/**
* \brief
* Function to calculate the total current based on Ampere's law while
* neglecting displacement current (J = curl x B). Used in the Ohm's law
* solver (kinetic-fluid hybrid model).
*
* \param[in] Bfield Magnetic field from which the current is calculated.
* \param[in] edge_lengths Length of cell edges taking embedded boundaries into account
*/
void CalculateCurrentAmpere (
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Bfield,
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths
);
void CalculateCurrentAmpere (
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
const int lev
);
/**
* \brief
* Function to update the E-field using Ohm's law (hybrid-PIC model).
*/
void HybridPICSolveE (
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Efield,
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Jfield,
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Bfield,
amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield,
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
DtType dt_type);
void HybridPICSolveE (
std::array< std::unique_ptr<amrex::MultiFab>, 3>& Efield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
std::unique_ptr<amrex::MultiFab> const& rhofield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
const int lev, DtType dt_type);
void HybridPICSolveE (
std::array< std::unique_ptr<amrex::MultiFab>, 3>& Efield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jfield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
std::unique_ptr<amrex::MultiFab> const& rhofield,
std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
const int lev, PatchType patch_type, DtType dt_type);
/**
* \brief
* Function to calculate the electron pressure at a given timestep type
* using the simulation charge density. Used in the Ohm's law solver
* (kinetic-fluid hybrid model).
*/
void CalculateElectronPressure ( DtType a_dt_type);
void CalculateElectronPressure (const int lev, DtType a_dt_type);
/**
* \brief Fill the electron pressure multifab given the kinetic particle
* charge density (and assumption of quasi-neutrality) using the user
* specified electron equation of state.
*
* \param[out] Pe scalar electron pressure MultiFab at a given level
* \param[in] rhofield scalar ion chrge density Multifab at a given level
*/
void FillElectronPressureMF (
std::unique_ptr<amrex::MultiFab> const& Pe,
amrex::MultiFab* const& rhofield );
// Declare variables to hold hybrid-PIC model parameters
/** Number of substeps to take when evolving B */
int m_substeps = 100;
/** Electron temperature in eV */
amrex::Real m_elec_temp;
/** Reference electron density */
amrex::Real m_n0_ref = 1.0;
/** Electron pressure scaling exponent */
amrex::Real m_gamma = 5.0/3.0;
/** Plasma density floor - if n < n_floor it will be set to n_floor */
amrex::Real m_n_floor = 1.0;
/** Plasma resistivity */
std::string m_eta_expression = "0.0";
std::unique_ptr<amrex::Parser> m_resistivity_parser;
amrex::ParserExecutor<1> m_eta;
// Declare multifabs specifically needed for the hybrid-PIC model
amrex::Vector< std::unique_ptr<amrex::MultiFab> > rho_fp_temp;
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_temp;
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_ampere;
amrex::Vector< std::unique_ptr<amrex::MultiFab> > electron_pressure_fp;
// Helper functions to retrieve hybrid-PIC multifabs
amrex::MultiFab * get_pointer_current_fp_ampere (int lev, int direction) const { return current_fp_ampere[lev][direction].get(); }
amrex::MultiFab * get_pointer_electron_pressure_fp (int lev) const { return electron_pressure_fp[lev].get(); }
/** Gpu Vector with index type of the Jx multifab */
amrex::GpuArray<int, 3> Jx_IndexType;
/** Gpu Vector with index type of the Jy multifab */
amrex::GpuArray<int, 3> Jy_IndexType;
/** Gpu Vector with index type of the Jz multifab */
amrex::GpuArray<int, 3> Jz_IndexType;
/** Gpu Vector with index type of the Bx multifab */
amrex::GpuArray<int, 3> Bx_IndexType;
/** Gpu Vector with index type of the By multifab */
amrex::GpuArray<int, 3> By_IndexType;
/** Gpu Vector with index type of the Bz multifab */
amrex::GpuArray<int, 3> Bz_IndexType;
/** Gpu Vector with index type of the Ex multifab */
amrex::GpuArray<int, 3> Ex_IndexType;
/** Gpu Vector with index type of the Ey multifab */
amrex::GpuArray<int, 3> Ey_IndexType;
/** Gpu Vector with index type of the Ez multifab */
amrex::GpuArray<int, 3> Ez_IndexType;
};
/**
* \brief
* This struct contains only static functions to compute the electron pressure
* using the particle density at a given point and the user provided reference
* density and temperatures.
*/
struct ElectronPressure {
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real get_pressure (amrex::Real const n0,
amrex::Real const T0,
amrex::Real const gamma,
amrex::Real const rho) {
return n0 * T0 * pow((rho/PhysConst::q_e)/n0, gamma);
}
};
#endif // WARPX_HYBRIDPICMODEL_H_
|
