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