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