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 4 : HeatConduction(int a_nghost = 2) :
52 : Integrator(),
53 4 : number_of_ghost_cells(a_nghost)
54 4 : {}
55 :
56 : // Constructor that triggers parse
57 4 : HeatConduction(IO::ParmParse& pp) : HeatConduction()
58 : {
59 4 : Parse(*this, pp);
60 4 : }
61 :
62 8 : virtual ~HeatConduction()
63 4 : {
64 4 : delete ic;
65 4 : delete bc;
66 8 : }
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 4 : 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 24 : pp_query_required("heat.alpha", value.alpha);
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 20 : pp_query_default("heat.refinement_threshold", value.refinement_threshold, 0.01);
86 :
87 : // Initial condition type.
88 12 : pp.select_default<IC::Constant,IC::Sphere,IC::Expression>("ic",value.ic,value.geom);
89 :
90 : // Select BC object for temperature
91 12 : pp.select_default<BC::Constant,BC::Expression>("bc.temp",value.bc,1);
92 :
93 : // Select between using a realspace solve or the spectral method
94 28 : pp.query_validate("method",value.method,{"realspace","spectral"});
95 :
96 : // Register the temperature and old temperature fields.
97 : // temp_mf and temp_old_mf are defined near the bottom of this Header file.
98 12 : value.RegisterNewFab(value.temp_mf, value.bc, value.number_of_components, value.number_of_ghost_cells, "Temp", true);
99 4 : if (value.method == "realspace")
100 : {
101 12 : value.RegisterNewFab(value.temp_old_mf, value.bc, value.number_of_components, value.number_of_ghost_cells, "Temp_old", false);
102 : }
103 4 : }
104 :
105 : protected:
106 :
107 : // Use the ic object to initialize the temperature field
108 15 : void Initialize(int lev)
109 : {
110 15 : ic->Initialize(lev, temp_mf);
111 15 : if (method == "realspace") ic->Initialize(lev, temp_old_mf);
112 15 : }
113 :
114 : // Integrate the heat equation
115 57560 : void Advance(int lev, Set::Scalar time, Set::Scalar dt)
116 : {
117 : // If we are solving using the spectral method, go there instead.
118 57560 : if (method == "spectral")
119 : {
120 0 : AdvanceSpectral(lev,time,dt);
121 0 : return;
122 : }
123 :
124 : // Swap the old temp fab and the new temp fab so we use
125 : // the new one.
126 57560 : std::swap(*temp_mf[lev], *temp_old_mf[lev]);
127 :
128 : // Get the cell size corresponding to this level
129 57560 : const Set::Scalar* DX = geom[lev].CellSize();
130 :
131 : // Iterate over all of the patches on this level
132 223148 : for (amrex::MFIter mfi(*temp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
133 : {
134 : // Get the box (index dimensions) for this patch
135 165588 : const amrex::Box &bx = mfi.tilebox();
136 :
137 : // Get an array-accessible handle to the data on this patch.
138 165588 : Set::Patch<const Set::Scalar> temp_old = temp_old_mf.Patch(lev,mfi);
139 165588 : Set::Patch<Set::Scalar> temp = temp_mf.Patch(lev,mfi);
140 :
141 : // Iterate over the grid on this patch
142 165588 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
143 : {
144 : // Do the physics!
145 : // Note that Numeric::Laplacian is an inlined function so there is no overhead.
146 : // You can calculate the derivatives yourself if you want.
147 658453920 : temp(i, j, k) = temp_old(i, j, k) + dt * alpha * Numeric::Laplacian(temp_old, i, j, k, 0, DX);
148 219484640 : });
149 57560 : }
150 : }
151 :
152 :
153 : #ifdef ALAMO_FFT
154 : void
155 : AdvanceSpectral(int lev, Set::Scalar /*time*/, Set::Scalar dt)
156 : {
157 : Util::Assert(INFO,TEST(lev == 0), "Only single level currently supported");
158 :
159 : amrex::Box const & domain = this->geom[lev].Domain();
160 :
161 : amrex::FFT::R2C my_fft(this->geom[lev].Domain());
162 : auto const& [cba, cdm] = my_fft.getSpectralDataLayout();
163 : amrex::FabArray<amrex::BaseFab<amrex::GpuComplex<Set::Scalar> > > Temp_hat(cba, cdm, 1, 0);
164 : my_fft.forward(*temp_mf[lev], Temp_hat);
165 :
166 : const Set::Scalar* DX = geom[lev].CellSize();
167 :
168 : Set::Scalar
169 : AMREX_D_DECL(
170 : pi_Lx = 2.0 * Set::Constant::Pi / geom[lev].Domain().length(0) / DX[0],
171 : pi_Ly = 2.0 * Set::Constant::Pi / geom[lev].Domain().length(1) / DX[1],
172 : pi_Lz = 2.0 * Set::Constant::Pi / geom[lev].Domain().length(2) / DX[2]);
173 :
174 : Set::Scalar scaling = 1.0 / geom[lev].Domain().d_numPts();
175 :
176 : for (amrex::MFIter mfi(Temp_hat, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
177 : {
178 : const amrex::Box &bx = mfi.tilebox();
179 : amrex::Array4<amrex::GpuComplex<Set::Scalar>> const & T_hat = Temp_hat.array(mfi);
180 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int m, int n, int p) {
181 :
182 : AMREX_D_TERM(
183 : Set::Scalar k1 = m * pi_Lx;,
184 : Set::Scalar k2 = (n < domain.length(1)/2 ? n * pi_Ly : (n - domain.length(1)) * pi_Ly);,
185 : Set::Scalar k3 = (p < domain.length(2)/2 ? p * pi_Lz : (p - domain.length(2)) * pi_Lz););
186 :
187 : Set::Scalar lap = AMREX_D_TERM(k1 * k1, + k2 * k2, + k3*k3);
188 :
189 : Set::Scalar factor = exp( - alpha * dt * lap);
190 : T_hat(m,n,p) *= factor*scaling;
191 : });
192 : }
193 :
194 : my_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 583 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/)
207 : {
208 583 : if (method=="spectral") return;
209 :
210 : // Get cell dimensions as done above.
211 583 : const Set::Scalar* DX = geom[lev].CellSize();
212 : // Calculate the diagonal.
213 583 : 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 5371 : for (amrex::MFIter mfi(*temp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
217 : {
218 : // Get the box and handles as done above.
219 4788 : const amrex::Box& bx = mfi.tilebox();
220 4788 : amrex::Array4<char> const& tags = a_tags.array(mfi);
221 4788 : amrex::Array4<Set::Scalar> const& temp = (*temp_mf[lev]).array(mfi);
222 :
223 : // Iterate over the grid as done above.
224 4788 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
225 : {
226 : // Calculate the temperature gradient.
227 408872 : 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 408872 : if (grad.lpNorm<2>() * dr > refinement_threshold)
232 601732 : tags(i, j, k) = amrex::TagBox::SET;
233 408872 : });
234 583 : }
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
|
