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
|
/* Copyright 2023 The WarpX Community
*
* This file is part of WarpX.
*
* Authors: Roelof Groenewald (TAE Technologies)
*
* License: BSD-3-Clause-LBNL
*/
#include "Evolve/WarpXDtType.H"
#include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H"
#include "Particles/MultiParticleContainer.H"
#include "Utils/TextMsg.H"
#include "Utils/WarpXProfilerWrapper.H"
#include "WarpX.H"
using namespace amrex;
void WarpX::HybridPICEvolveFields ()
{
WARPX_PROFILE("WarpX::HybridPICEvolveFields()");
// The below deposition is hard coded for a single level simulation
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
finest_level == 0,
"Ohm's law E-solve only works with a single level.");
// The particles have now been pushed to their t_{n+1} positions.
// Perform charge deposition in component 0 of rho_fp at t_{n+1}.
mypc->DepositCharge(rho_fp, 0._rt);
// Perform current deposition at t_{n+1/2}.
mypc->DepositCurrent(current_fp, dt[0], -0.5_rt * dt[0]);
// Synchronize J and rho:
// filter (if used), exchange guard cells, interpolate across MR levels
// and apply boundary conditions
SyncCurrentAndRho();
// SyncCurrent does not include a call to FillBoundary, but it is needed
// for the hybrid-PIC solver since current values are interpolated to
// a nodal grid
for (int lev = 0; lev <= finest_level; ++lev)
for (int idim = 0; idim < 3; ++idim)
current_fp[lev][idim]->FillBoundary(Geom(lev).periodicity());
// Get requested number of substeps to use
int sub_steps = m_hybrid_pic_model->m_substeps / 2;
// Reference hybrid-PIC multifabs
auto& rho_fp_temp = m_hybrid_pic_model->rho_fp_temp;
auto& current_fp_temp = m_hybrid_pic_model->current_fp_temp;
// During the above deposition the charge and current density were updated
// so that, at this time, we have rho^{n} in rho_fp_temp, rho{n+1} in the
// 0'th index of `rho_fp`, J_i^{n-1/2} in `current_fp_temp` and J_i^{n+1/2}
// in `current_fp`.
// TODO: To speed up the algorithm insert Runge-Kutta integration logic
// for B update instead of the substep update used here - can test with
// small timestep using this simpler implementation
// Note: E^{n} is recalculated with the accurate J_i^{n} since at the end
// of the last step we had to "guess" it. It also needs to be
// recalculated to include the resistivity before evolving B.
// J_i^{n} is calculated as the average of J_i^{n-1/2} and J_i^{n+1/2}.
for (int lev = 0; lev <= finest_level; ++lev)
{
for (int idim = 0; idim < 3; ++idim) {
// Perform a linear combination of values in the 0'th index (1 comp)
// of J_i^{n-1/2} and J_i^{n+1/2} (with 0.5 prefactors), writing
// the result into the 0'th index of `current_fp_temp[lev][idim]`
MultiFab::LinComb(
*current_fp_temp[lev][idim],
0.5_rt, *current_fp_temp[lev][idim], 0,
0.5_rt, *current_fp[lev][idim], 0,
0, 1, current_fp_temp[lev][idim]->nGrowVect()
);
}
}
// Calculate the electron pressure at t=n using rho^n
m_hybrid_pic_model->CalculateElectronPressure(DtType::FirstHalf);
// Push the B field from t=n to t=n+1/2 using the current and density
// at t=n, while updating the E field along with B using the electron
// momentum equation
for (int sub_step = 0; sub_step < sub_steps; sub_step++)
{
m_hybrid_pic_model->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths);
m_hybrid_pic_model->HybridPICSolveE(
Efield_fp, current_fp, Bfield_fp, rho_fp, m_edge_lengths,
DtType::FirstHalf
);
FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
EvolveB(0.5 / sub_steps * dt[0], DtType::FirstHalf);
FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
}
// Average rho^{n} and rho^{n+1} to get rho^{n+1/2} in rho_fp_temp
for (int lev = 0; lev <= finest_level; ++lev)
{
// Perform a linear combination of values in the 0'th index (1 comp)
// of rho^{n} and rho^{n+1} (with 0.5 prefactors), writing
// the result into the 0'th index of `rho_fp_temp[lev]`
MultiFab::LinComb(
*rho_fp_temp[lev], 0.5_rt, *rho_fp_temp[lev], 0,
0.5_rt, *rho_fp[lev], 0, 0, 1, rho_fp_temp[lev]->nGrowVect()
);
}
// Calculate the electron pressure at t=n+1/2
m_hybrid_pic_model->CalculateElectronPressure(DtType::SecondHalf);
// Now push the B field from t=n+1/2 to t=n+1 using the n+1/2 quantities
for (int sub_step = 0; sub_step < sub_steps; sub_step++)
{
m_hybrid_pic_model->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths);
m_hybrid_pic_model->HybridPICSolveE(
Efield_fp, current_fp, Bfield_fp, rho_fp, m_edge_lengths,
DtType::SecondHalf
);
FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
EvolveB(0.5 / sub_steps * dt[0], DtType::SecondHalf);
FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
}
// Extrapolate the ion current density to t=n+1 using
// J_i^{n+1} = 1/2 * J_i^{n-1/2} + 3/2 * J_i^{n+1/2}, and recalling that
// now current_fp_temp = J_i^{n} = 1/2 * (J_i^{n-1/2} + J_i^{n+1/2})
for (int lev = 0; lev <= finest_level; ++lev)
{
for (int idim = 0; idim < 3; ++idim) {
// Perform a linear combination of values in the 0'th index (1 comp)
// of J_i^{n-1/2} and J_i^{n+1/2} (with -1.0 and 2.0 prefactors),
// writing the result into the 0'th index of `current_fp_temp[lev][idim]`
MultiFab::LinComb(
*current_fp_temp[lev][idim],
-1._rt, *current_fp_temp[lev][idim], 0,
2._rt, *current_fp[lev][idim], 0,
0, 1, current_fp_temp[lev][idim]->nGrowVect()
);
}
}
// Calculate the electron pressure at t=n+1
m_hybrid_pic_model->CalculateElectronPressure(DtType::Full);
// Update the E field to t=n+1 using the extrapolated J_i^n+1 value
m_hybrid_pic_model->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths);
m_hybrid_pic_model->HybridPICSolveE(
Efield_fp, current_fp, Bfield_fp, rho_fp, m_edge_lengths,
DtType::Full
);
FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
// Copy the rho^{n+1} values to rho_fp_temp and the J_i^{n+1/2} values to
// current_fp_temp since at the next step those values will be needed as
// rho^{n} and J_i^{n-1/2}.
for (int lev = 0; lev <= finest_level; ++lev)
{
// copy 1 component value starting at index 0 to index 0
MultiFab::Copy(*rho_fp_temp[lev], *rho_fp[lev],
0, 0, 1, rho_fp_temp[lev]->nGrowVect());
for (int idim = 0; idim < 3; ++idim) {
MultiFab::Copy(*current_fp_temp[lev][idim], *current_fp[lev][idim],
0, 0, 1, current_fp_temp[lev][idim]->nGrowVect());
}
}
}
|
