Line data Source code
1 : //
2 : // This implements a basic heat conduction method in Alamo.
3 : // The partial differential equation to be solved is
4 : //
5 : // .. math::
6 : //
7 : // \frac{\partial T}{\partial t} = \alpha\,\Delta T
8 : //
9 : // where :math:`T` is temperature, :math:`t` is time, and :math:`alpha`
10 : // is the thermal diffusivity.
11 : // Integration is performed explicitly in time using forward Euler, and
12 : // differentiation is performed using the finite difference method.
13 : //
14 :
15 : #ifndef INTEGRATOR_HEATCONDUCTION_H // Include guards
16 : #define INTEGRATOR_HEATCONDUCTION_H //
17 :
18 : // AMReX Includes
19 : #include "AMReX_Array4.H"
20 : #include "AMReX_GpuComplex.H"
21 : #include "AMReX_MFIter.H"
22 : #include <AMReX.H>
23 : #include <AMReX_ParallelDescriptor.H>
24 : #include <AMReX_ParmParse.H>
25 : #include <AMReX_TimeIntegrator.H>
26 : #ifdef ALAMO_FFT
27 : #include <AMReX_FFT.H>
28 : #endif
29 :
30 : // Alamo Includes
31 : #include "Set/Base.H"
32 : #include "Integrator.H"
33 : #include "IO/ParmParse.H"
34 : #include "Integrator/Integrator.H"
35 : #include "BC/Constant.H"
36 : #include "BC/Expression.H"
37 : #include "IC/IC.H"
38 : #include "IC/Sphere.H"
39 : #include "IC/Constant.H"
40 : #include "IC/Expression.H"
41 : #include "Numeric/Stencil.H"
42 :
43 : namespace Integrator
44 : {
45 : class HeatConduction : virtual public Integrator
46 : {
47 : public:
48 : static constexpr const char* name = "heatconduction";
49 :
50 :
51 : // Empty constructor
52 4 : HeatConduction(int a_nghost = 2) :
53 : Integrator(),
54 4 : number_of_ghost_cells(a_nghost)
55 4 : {}
56 :
57 : // Constructor that triggers parse
58 4 : HeatConduction(IO::ParmParse& pp) : HeatConduction()
59 : {
60 4 : Parse(*this, pp);
61 4 : }
62 :
63 8 : virtual ~HeatConduction()
64 4 : {
65 4 : delete ic;
66 4 : delete bc;
67 8 : }
68 :
69 : // The Parse function initializes the HeatConduction object using
70 : // a parser, pp.
71 : // Note that this is a static function, which means it does not have
72 : // direct access to member variables. Instead, it initializes the variables
73 : // inside the argument, "value", and so all references to member items are
74 : // prefixed by "value."
75 4 : static void Parse(HeatConduction& value, IO::ParmParse& pp)
76 : {
77 : // Diffusion coefficient :math:`\alpha`.
78 : // *This is an example of a required input variable -
79 : // - program will terminate unless it is provided.*
80 8 : pp.query_required( "heat.alpha", value.alpha,
81 : Unit::ThermalDiffusivity());
82 :
83 : // Criterion for mesh refinement.
84 : // *This is an example of a default input variable.
85 : // The default value is provided here, not in the
86 : // definition of the variable.*
87 16 : pp.query_default( "heat.refinement_threshold", value.refinement_threshold, "0.01_K",
88 : Unit::Temperature());
89 :
90 : // Initial condition type.
91 8 : pp.select_default<IC::Constant,IC::Sphere,IC::Expression>( "ic",value.ic,value.geom,
92 4 : Unit::Temperature());
93 :
94 : // Select BC object for temperature
95 8 : pp.select_default<BC::Constant,BC::Expression>( "bc.temp",value.bc,1,
96 4 : Unit::Temperature());
97 :
98 : // Select between using a realspace solve or the spectral method
99 16 : pp.query_validate("method",value.method,{"realspace","spectral"});
100 :
101 : // Register the temperature and old temperature fields.
102 : // temp_mf and temp_old_mf are defined near the bottom of this Header file.
103 16 : value.RegisterNewFab( value.temp_mf, value.bc, value.number_of_components,
104 4 : value.number_of_ghost_cells, "Temp", true);
105 4 : if (value.method == "realspace")
106 : {
107 16 : value.RegisterNewFab( value.temp_old_mf, value.bc, value.number_of_components,
108 4 : value.number_of_ghost_cells, "Temp_old", false);
109 : }
110 4 : }
111 :
112 : protected:
113 :
114 : // Use the ic object to initialize the temperature field
115 16 : void Initialize(int lev)
116 : {
117 16 : ic->Initialize(lev, temp_mf);
118 16 : if (method == "realspace") ic->Initialize(lev, temp_old_mf);
119 16 : }
120 :
121 : // Integrate the heat equation
122 18000 : void Advance(int lev, Set::Scalar time, Set::Scalar dt)
123 : {
124 : // If we are solving using the spectral method, go there instead.
125 18000 : if (method == "spectral")
126 : {
127 0 : AdvanceSpectral(lev,time,dt);
128 0 : return;
129 : }
130 :
131 : // Swap the old temp fab and the new temp fab so we use
132 : // the new one.
133 18000 : std::swap(*temp_mf[lev], *temp_old_mf[lev]);
134 :
135 : // Create the time integrator object
136 18000 : amrex::TimeIntegrator timeintegrator(*temp_old_mf[lev], time);
137 :
138 : // Calculate the "RHS" of the function.
139 : // If using forward Euler, the RHS would be
140 : // T_new = T_old + dt * RHS
141 18000 : timeintegrator.set_rhs([&](amrex::MultiFab & rhs_mf, amrex::MultiFab & temp_mf, const Set::Scalar /*time*/)
142 : {
143 : // Get the cell size corresponding to this level
144 25500 : const Set::Scalar* DX = geom[lev].CellSize();
145 :
146 : // Iterate over all of the patches on this level
147 412860 : for (amrex::MFIter mfi(temp_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
148 : {
149 : // Get the box (index dimensions) for this patch
150 387360 : const amrex::Box &bx = mfi.tilebox();
151 :
152 : // Get an array-accessible handle to the data on this patch.
153 387360 : Set::Patch<const Set::Scalar> temp = temp_mf.array(mfi);
154 387360 : Set::Patch<Set::Scalar> rhs = rhs_mf.array(mfi);
155 :
156 : // Iterate over the grid on this patch
157 387360 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
158 : {
159 : // Do the physics!
160 : // Note that Numeric::Laplacian is an inlined function so there is no overhead.
161 : // You can calculate the derivatives yourself if you want.
162 93146880 : rhs(i, j, k) = alpha * Numeric::Laplacian(temp, i, j, k, 0, DX);
163 46573440 : });
164 25500 : }
165 25500 : });
166 :
167 : // We must update the boundaries here and apply boundary conditions
168 18000 : timeintegrator.set_post_stage_action([&](amrex::MultiFab & stage_mf, Set::Scalar time)
169 : {
170 7500 : bc->FillBoundary(stage_mf,0,1,time,0);
171 7500 : stage_mf.FillBoundary(true);
172 7500 : });
173 :
174 : // Do the update
175 18000 : timeintegrator.advance(*temp_old_mf[lev], *temp_mf[lev], time, dt);
176 18000 : }
177 :
178 :
179 : #ifdef ALAMO_FFT
180 : void
181 : AdvanceSpectral(int lev, Set::Scalar /*time*/, Set::Scalar dt)
182 : {
183 : Util::Assert(INFO,TEST(lev == 0), "Only single level currently supported");
184 :
185 : amrex::Box const & domain = this->geom[lev].Domain();
186 :
187 : amrex::FFT::R2C my_fft(this->geom[lev].Domain());
188 : auto const& [cba, cdm] = my_fft.getSpectralDataLayout();
189 : amrex::FabArray<amrex::BaseFab<amrex::GpuComplex<Set::Scalar> > > Temp_hat(cba, cdm, 1, 0);
190 : my_fft.forward(*temp_mf[lev], Temp_hat);
191 :
192 : const Set::Scalar* DX = geom[lev].CellSize();
193 :
194 : Set::Scalar
195 : AMREX_D_DECL(
196 : pi_Lx = 2.0 * Set::Constant::Pi / geom[lev].Domain().length(0) / DX[0],
197 : pi_Ly = 2.0 * Set::Constant::Pi / geom[lev].Domain().length(1) / DX[1],
198 : pi_Lz = 2.0 * Set::Constant::Pi / geom[lev].Domain().length(2) / DX[2]);
199 :
200 : Set::Scalar scaling = 1.0 / geom[lev].Domain().d_numPts();
201 :
202 : for (amrex::MFIter mfi(Temp_hat, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
203 : {
204 : const amrex::Box &bx = mfi.tilebox();
205 : amrex::Array4<amrex::GpuComplex<Set::Scalar>> const & T_hat = Temp_hat.array(mfi);
206 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int m, int n, int p) {
207 :
208 : AMREX_D_TERM(
209 : Set::Scalar k1 = m * pi_Lx;,
210 : Set::Scalar k2 = (n < domain.length(1)/2 ? n * pi_Ly : (n - domain.length(1)) * pi_Ly);,
211 : Set::Scalar k3 = (p < domain.length(2)/2 ? p * pi_Lz : (p - domain.length(2)) * pi_Lz););
212 :
213 : Set::Scalar lap = AMREX_D_TERM(k1 * k1, + k2 * k2, + k3*k3);
214 :
215 : Set::Scalar factor = exp( - alpha * dt * lap);
216 : T_hat(m,n,p) *= factor*scaling;
217 : });
218 : }
219 :
220 : my_fft.backward(Temp_hat, *temp_mf[lev]);
221 : }
222 : #else
223 : void
224 0 : AdvanceSpectral(int, Set::Scalar, Set::Scalar)
225 : {
226 0 : Util::Abort(INFO,"Alamo must be configured with --fft");
227 0 : }
228 : #endif
229 :
230 :
231 : // Tag cells for mesh refinement based on temperature gradient
232 1344 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/)
233 : {
234 1344 : if (method=="spectral") return;
235 :
236 : // Get cell dimensions as done above.
237 1344 : const Set::Scalar* DX = geom[lev].CellSize();
238 : // Calculate the diagonal.
239 1344 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
240 :
241 : // Iterate over the patches on this level
242 12478 : for (amrex::MFIter mfi(*temp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
243 : {
244 : // Get the box and handles as done above.
245 11134 : const amrex::Box& bx = mfi.tilebox();
246 11134 : amrex::Array4<char> const& tags = a_tags.array(mfi);
247 11134 : amrex::Array4<Set::Scalar> const& temp = (*temp_mf[lev]).array(mfi);
248 :
249 : // Iterate over the grid as done above.
250 11134 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
251 : {
252 : // Calculate the temperature gradient.
253 784036 : Set::Vector grad = Numeric::Gradient(temp, i, j, k, 0, DX);
254 :
255 : // Is the gradient * cell_size too big? If so, then
256 : // mark this cell as needing refinement.
257 784036 : if (grad.lpNorm<2>() * dr > refinement_threshold)
258 1329706 : tags(i, j, k) = amrex::TagBox::SET;
259 784036 : });
260 1344 : }
261 : }
262 :
263 : protected:
264 : Set::Field<Set::Scalar> temp_mf; // Temperature field variable (current timestep)
265 : Set::Field<Set::Scalar> temp_old_mf; // Temperature field variable (previous timestep)
266 :
267 : std::string method; // determine whether to use realspace or spectral method
268 :
269 : private:
270 :
271 : //
272 : // Definition of parameters set only at instantiation by
273 : // constructors.
274 : //
275 : const int number_of_components = 1; // Number of components
276 : const int number_of_ghost_cells = 2; // Number of ghost cells
277 :
278 : //
279 : // Definition of user-determined variables.
280 : //
281 : // Instantiate all variables to NAN if possible.
282 : // Default values may be set in Parse using query_default.
283 : //
284 :
285 : Set::Scalar alpha = NAN; // Thermal diffusivity
286 : Set::Scalar refinement_threshold = NAN ; // Criterion for cell refinement
287 :
288 : //
289 : // Definition of user-determined pointer variables.
290 : //
291 : // These should be set to nullptr. Make sure that they are deleted
292 : // in the ~HeatConduction destructor.
293 : //
294 :
295 : IC::IC<Set::Scalar>* ic = nullptr; // Object used to initialize temperature field
296 : BC::BC<Set::Scalar>* bc = nullptr; // Object used to update temp field boundary ghost cells
297 : };
298 : } // namespace Integrator
299 : #endif
|
