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 : // Standard library includes
19 : #include <iostream>
20 : #include <fstream>
21 : #include <iomanip>
22 :
23 : // AMReX Includes
24 : #include "AMReX.H"
25 : #include "AMReX_ParallelDescriptor.H"
26 : #include "AMReX_ParmParse.H"
27 :
28 : // Alamo Includes
29 : #include "IO/ParmParse.H"
30 : #include "Integrator/Integrator.H"
31 : #include "BC/Constant.H"
32 : #include "BC/Expression.H"
33 : #include "IC/IC.H"
34 : #include "IC/Sphere.H"
35 : #include "IC/Constant.H"
36 : #include "IC/Expression.H"
37 : #include "Numeric/Stencil.H"
38 :
39 : namespace Integrator
40 : {
41 : class HeatConduction : virtual public Integrator
42 : {
43 : public:
44 : // Empty constructor
45 3 : HeatConduction(int a_nghost = 2) :
46 : Integrator(),
47 3 : number_of_ghost_cells(a_nghost)
48 3 : {}
49 :
50 : // Constructor that triggers parse
51 3 : HeatConduction(IO::ParmParse& pp) : HeatConduction()
52 : {
53 3 : Parse(*this, pp);
54 3 : }
55 :
56 6 : virtual ~HeatConduction()
57 3 : {
58 3 : delete ic;
59 3 : delete bc;
60 6 : }
61 :
62 : // The Parse function initializes the HeatConduction object using
63 : // a parser, pp.
64 : // Note that this is a static function, which means it does not have
65 : // direct access to member variables. Instead, it initializes the variables
66 : // inside the argument, "value", and so all references to member items are
67 : // prefixed by "value."
68 3 : static void Parse(HeatConduction& value, IO::ParmParse& pp)
69 : {
70 : // Diffusion coefficient :math:`\alpha`.
71 : // *This is an example of a required input variable -
72 : // - program will terminate unless it is provided.*
73 3 : pp_query_required("heat.alpha", value.alpha);
74 :
75 : // Criterion for mesh refinement.
76 : // *This is an example of a default input variable.
77 : // The default value is provided here, not in the
78 : // definition of the variable.*
79 3 : pp_query_default("heat.refinement_threshold", value.refinement_threshold, 0.01);
80 :
81 6 : std::string type;
82 :
83 : // Initial condition type.
84 : // *This is an example of type validation:
85 : // the input variable must be one of the three provided values,
86 : // and will error if not. The default selection is the first
87 : // argument.*
88 3 : pp_query_validate("ic.type", type, {"constant", "sphere","expression"});
89 3 : if (type == "sphere") value.ic = new IC::Sphere(value.geom, pp, "ic.sphere");
90 0 : else if (type == "constant") value.ic = new IC::Constant(value.geom, pp, "ic.constant");
91 0 : else if (type == "expression") value.ic = new IC::Expression(value.geom, pp, "ic.expression");
92 0 : else Util::Abort(INFO, "Invalid ic.type ", type);
93 :
94 3 : std::string bc_type;
95 : // Select BC object for temperature
96 3 : pp_query_validate("bc.temp.type",bc_type,{"constant","expression"});
97 3 : if (bc_type == "expression") value.bc = new BC::Expression(1, pp, "bc.temp.expression");
98 3 : else if (bc_type == "constant") value.bc = new BC::Constant(1, pp, "bc.temp");
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 3 : value.RegisterNewFab(value.temp_mf, value.bc, value.number_of_components, value.number_of_ghost_cells, "Temp", true);
103 3 : value.RegisterNewFab(value.temp_old_mf, value.bc, value.number_of_components, value.number_of_ghost_cells, "Temp_old", false);
104 3 : }
105 :
106 : protected:
107 :
108 : // Use the ic object to initialize the temperature field
109 14 : void Initialize(int lev)
110 : {
111 14 : ic->Initialize(lev, temp_old_mf);
112 14 : }
113 :
114 : // Integrate the heat equation
115 7560 : void Advance(int lev, Set::Scalar /*time*/, Set::Scalar dt)
116 : {
117 : // Swap the old temp fab and the new temp fab so we use
118 : // the new one.
119 7560 : std::swap(*temp_mf[lev], *temp_old_mf[lev]);
120 :
121 : // Get the cell size corresponding to this level
122 7560 : const Set::Scalar* DX = geom[lev].CellSize();
123 :
124 : // Iterate over all of the patches on this level
125 123148 : for (amrex::MFIter mfi(*temp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
126 : {
127 : // Get the box (index dimensions) for this patch
128 115588 : const amrex::Box& bx = mfi.tilebox();
129 :
130 : // Get an array-accessible handle to the data on this patch.
131 : // amrex::Array4<const Set::Scalar> const& temp_old = (*temp_old_mf[lev]).array(mfi);
132 : // amrex::Array4<Set::Scalar> const& temp = temp_mf.Patch(lev,mfi);//(*temp_mf[lev]).array(mfi);
133 115588 : Set::Patch<const Set::Scalar> temp_old = temp_old_mf.Patch(lev,mfi);
134 115588 : Set::Patch<Set::Scalar> temp = temp_mf.Patch(lev,mfi);
135 :
136 : // Iterate over the grid on this patch
137 115588 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
138 : {
139 : // Do the physics!
140 : // Note that Numeric::Laplacian is an inlined function so there is no overhead.
141 : // You can calculate the derivatives yourself if you want.
142 44053880 : temp(i, j, k) = temp_old(i, j, k) + dt * alpha * Numeric::Laplacian(temp_old, i, j, k, 0, DX);
143 14684628 : });
144 : }
145 7560 : }
146 :
147 : // Tag cells for mesh refinement based on temperature gradient
148 583 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/)
149 : {
150 : // Get cell dimensions as done above.
151 583 : const Set::Scalar* DX = geom[lev].CellSize();
152 : // Calculate the diagonal.
153 583 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
154 :
155 : // Iterate over the patches on this level
156 5371 : for (amrex::MFIter mfi(*temp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
157 : {
158 : // Get the box and handles as done above.
159 4788 : const amrex::Box& bx = mfi.tilebox();
160 4788 : amrex::Array4<char> const& tags = a_tags.array(mfi);
161 4788 : amrex::Array4<Set::Scalar> const& temp = (*temp_mf[lev]).array(mfi);
162 :
163 : // Iterate over the grid as done above.
164 4788 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
165 : {
166 : // Calculate the temperature gradient.
167 408872 : Set::Vector grad = Numeric::Gradient(temp, i, j, k, 0, DX);
168 :
169 : // Is the gradient * cell_size too big? If so, then
170 : // mark this cell as needing refinement.
171 408872 : if (grad.lpNorm<2>() * dr > refinement_threshold)
172 601732 : tags(i, j, k) = amrex::TagBox::SET;
173 408872 : });
174 : }
175 583 : }
176 :
177 : protected:
178 : Set::Field<Set::Scalar> temp_mf; // Temperature field variable (current timestep)
179 : Set::Field<Set::Scalar> temp_old_mf; // Temperature field variable (previous timestep)
180 :
181 : private:
182 :
183 : //
184 : // Definition of parameters set only at instantiation by
185 : // constructors.
186 : //
187 : const int number_of_components = 1; // Number of components
188 : const int number_of_ghost_cells = 2; // Number of ghost cells
189 :
190 : //
191 : // Definition of user-determined variables.
192 : //
193 : // Instantiate all variables to NAN if possible.
194 : // Default values may be set in Parse using query_default.
195 : //
196 :
197 : Set::Scalar alpha = NAN; // Thermal diffusivity
198 : Set::Scalar refinement_threshold = NAN ; // Criterion for cell refinement
199 :
200 : //
201 : // Definition of user-determined pointer variables.
202 : //
203 : // These should be set to nullptr. Make sure that they are deleted
204 : // in the ~HeatConduction destructor.
205 : //
206 :
207 : IC::IC* ic = nullptr; // Object used to initialize temperature field
208 : BC::BC<Set::Scalar>* bc = nullptr; // Object used to update temp field boundary ghost cells
209 : };
210 : } // namespace Integrator
211 : #endif
|
