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
|
#include "Interpolate.H"
#include <AMReX_FillPatchUtil_F.H>
namespace Interpolate
{
using namespace amrex;
std::unique_ptr<MultiFab>
getInterpolatedScalar(
const MultiFab& F_cp, const MultiFab& F_fp,
const DistributionMapping& dm, const int r_ratio,
const Real* /*dx*/, const int ngrow )
{
// Prepare the structure that will contain the returned fields
std::unique_ptr<MultiFab> interpolated_F;
interpolated_F.reset( new MultiFab(F_fp.boxArray(), dm, 1, ngrow) );
interpolated_F->setVal(0.);
// Loop through the boxes and interpolate the values from the _cp data
#ifdef _OPENMP
#pragma omp parallel
#endif
{
FArrayBox ffab; // Temporary array ; contains interpolated fields
for (MFIter mfi(*interpolated_F); mfi.isValid(); ++mfi)
{
Box finebx = mfi.fabbox();
finebx.coarsen(r_ratio).refine(r_ratio); // so that finebx is coarsenable
const FArrayBox& cfab = (F_cp)[mfi];
ffab.resize(finebx);
// - Fully nodal
if ( F_fp.is_nodal() ){
IntVect refinement_vector{AMREX_D_DECL(r_ratio, r_ratio, r_ratio)};
node_bilinear_interp.interp(cfab, 0, ffab, 0, 1,
finebx, refinement_vector, {}, {}, {}, 0, 0, RunOn::Cpu);
} else {
amrex::Abort("Unknown field staggering.");
}
// Add temporary array to the returned structure
const Box& bx = (*interpolated_F)[mfi].box();
(*interpolated_F)[mfi].plus<RunOn::Host>(ffab, bx, bx, 0, 0, 1);
}
}
return interpolated_F;
}
std::array<std::unique_ptr<MultiFab>, 3>
getInterpolatedVector(
const MultiFab* Fx_cp,
const MultiFab* Fy_cp,
const MultiFab* Fz_cp,
const MultiFab* Fx_fp,
const MultiFab* Fy_fp,
const MultiFab* Fz_fp,
const DistributionMapping& dm, const int r_ratio,
const Real* dx, const int ngrow )
{
// Prepare the structure that will contain the returned fields
std::array<std::unique_ptr<MultiFab>, 3> interpolated_F;
interpolated_F[0].reset( new MultiFab(Fx_fp->boxArray(), dm, 1, ngrow) );
interpolated_F[1].reset( new MultiFab(Fy_fp->boxArray(), dm, 1, ngrow) );
interpolated_F[2].reset( new MultiFab(Fz_fp->boxArray(), dm, 1, ngrow) );
for (int i=0; i<3; i++) interpolated_F[i]->setVal(0.);
// Loop through the boxes and interpolate the values from the _cp data
const int use_limiter = 0;
#ifdef _OPENMP
#pragma omp parallel
#endif
{
std::array<FArrayBox,3> ffab; // Temporary array ; contains interpolated fields
for (MFIter mfi(*interpolated_F[0]); mfi.isValid(); ++mfi)
{
Box ccbx = mfi.fabbox();
ccbx.enclosedCells();
ccbx.coarsen(r_ratio).refine(r_ratio); // so that ccbx is coarsenable
const FArrayBox& cxfab = (*Fx_cp)[mfi];
const FArrayBox& cyfab = (*Fy_cp)[mfi];
const FArrayBox& czfab = (*Fz_cp)[mfi];
ffab[0].resize(amrex::convert(ccbx,(*Fx_fp)[mfi].box().type()));
ffab[1].resize(amrex::convert(ccbx,(*Fy_fp)[mfi].box().type()));
ffab[2].resize(amrex::convert(ccbx,(*Fz_fp)[mfi].box().type()));
// - Face centered, in the same way as B on a Yee grid
if ( (*Fx_fp)[mfi].box().type() == IntVect{AMREX_D_DECL(1,0,0)} ){
#if (AMREX_SPACEDIM == 3)
amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(),
BL_TO_FORTRAN_ANYD(ffab[0]),
BL_TO_FORTRAN_ANYD(ffab[1]),
BL_TO_FORTRAN_ANYD(ffab[2]),
BL_TO_FORTRAN_ANYD(cxfab),
BL_TO_FORTRAN_ANYD(cyfab),
BL_TO_FORTRAN_ANYD(czfab),
dx, &r_ratio, &use_limiter);
#else
amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(),
BL_TO_FORTRAN_ANYD(ffab[0]),
BL_TO_FORTRAN_ANYD(ffab[2]),
BL_TO_FORTRAN_ANYD(cxfab),
BL_TO_FORTRAN_ANYD(czfab),
dx, &r_ratio, &use_limiter);
amrex_interp_cc_bfield(ccbx.loVect(), ccbx.hiVect(),
BL_TO_FORTRAN_ANYD(ffab[1]),
BL_TO_FORTRAN_ANYD(cyfab),
&r_ratio, &use_limiter);
#endif
// - Edge centered, in the same way as E on a Yee grid
} else if ( (*Fx_fp)[mfi].box().type() == IntVect{AMREX_D_DECL(0,1,1)} ){
#if (AMREX_SPACEDIM == 3)
amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(),
BL_TO_FORTRAN_ANYD(ffab[0]),
BL_TO_FORTRAN_ANYD(ffab[1]),
BL_TO_FORTRAN_ANYD(ffab[2]),
BL_TO_FORTRAN_ANYD(cxfab),
BL_TO_FORTRAN_ANYD(cyfab),
BL_TO_FORTRAN_ANYD(czfab),
&r_ratio, &use_limiter);
#else
amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(),
BL_TO_FORTRAN_ANYD(ffab[0]),
BL_TO_FORTRAN_ANYD(ffab[2]),
BL_TO_FORTRAN_ANYD(cxfab),
BL_TO_FORTRAN_ANYD(czfab),
&r_ratio,&use_limiter);
amrex_interp_nd_efield(ccbx.loVect(), ccbx.hiVect(),
BL_TO_FORTRAN_ANYD(ffab[1]),
BL_TO_FORTRAN_ANYD(cyfab),
&r_ratio);
#endif
} else {
amrex::Abort("Unknown field staggering.");
}
// Add temporary array to the returned structure
for (int i = 0; i < 3; ++i) {
const Box& bx = (*interpolated_F[i])[mfi].box();
(*interpolated_F[i])[mfi].plus<RunOn::Host>(ffab[i], bx, bx, 0, 0, 1);
}
}
}
return interpolated_F;
}
}
|
