Heat Equation Tutorial¶
This tutorial walks through the two-dimensional heat-equation example in
examples/heat_equation. It shows how a FleCSI mesh and fields become
flecsolve vectors, how the discrete Laplacian is exposed as an operator,
and how explicit and implicit time integrators use that operator.
The complete, buildable source for this tutorial is in the repository under
examples/heat_equation. The snippets below are included directly from
that source so that the tutorial stays aligned with the example.
Problem¶
The example solves the heat equation
with homogeneous Dirichlet boundary conditions and a hot square as the initial condition:
The mesh is distributed across flecsi colors. The command-line mesh
extents specify the number of points in the x and y
directions.
Build the example¶
Configure the project with examples enabled and build it:
cmake -S . -B build -DFLECSOLVE_BUILD_EXAMPLES=ON
cmake --build build
The executables and their configuration files are written to
build/examples/heat_equation.
Run the example¶
Run the explicit RK23 version on four MPI processes:
cd build/examples/heat_equation
mpirun -np 4 ./heat-explicit 100 100 -d 1.5 -o true
Run the implicit BDF version with the same mesh and diffusivity:
mpirun -np 4 ./heat-implicit 100 100 -d 1.5 -o true
The positional arguments are the x and y mesh extents. The -d
option sets the diffusivity alpha and -o true writes every accepted
time step. The explicit.cfg and implicit.cfg files control the time
integrator and linear-solver settings.
The final output is written by the finalize control point. With
-o true, intermediate files are also written using names such as
timestep-0-0.dat. The files contain x, y, and u columns and
can be visualized with the utilities in examples/heat_equation/util.
Mesh and vectors¶
The example specializes FleCSI’s narray topology as a two-dimensional
mesh. Two fields hold the current and next solution values. The control
policy turns those fields into flecsolve topology views after the mesh has
been allocated:
struct control_policy : flecsi::run::control_base {
using control_points_enum = cp;
using control = flecsi::run::control<control_policy>;
using control_points =
list<point<cp::initialize>, point<cp::advance>, point<cp::finalize>>;
double diffusivity;
heat::mesh::ptr m;
auto & mesh() { return *m; }
using vec = decltype(flecsolve::vec::make(ud[0](*m)));
vec & u() { return u_.value(); }
vec & unew() { return unew_.value(); }
void initialize_vectors() {
u_.emplace(flecsolve::vec::make(ud[0](mesh())));
unew_.emplace(flecsolve::vec::make(ud[1](mesh())));
}
template<class T>
void save_geometry(const T & g, std::vector<std::size_t> extents) {
dx = std::abs(g[0][1] - g[0][0]) / (extents[0] - 1);
dy = std::abs(g[1][1] - g[1][0]) / (extents[1] - 1);
}
The topology and field definitions are in
examples/heat_equation/mesh.hh. A topology view lets the time integrator
operate on FleCSI fields through the standard flecsolve vector interface.
Initial and boundary conditions¶
The initialize control point first allocates the mesh and then launches
the ics task. That task sets the value to 50 inside the square
[4, 6] x [4, 6] and to zero elsewhere. The implementation is in
examples/heat_equation/heat.cc.
The discrete Laplacian applies the zero Dirichlet boundary condition before computing interior values. The boundary checks account for MPI subdomains, so only processes that own a global boundary write those boundary values.
Discrete operator¶
The laplace task uses the mesh spacing and a centered finite-difference
stencil. It computes alpha * Laplacian(u) into the output vector:
inline void laplace(mesh::accessor<ro> m,
const double c,
field<double>::accessor<wo, na> unewa,
field<double>::accessor<rw, ro> ua) {
auto unew = m.mdcolex<mesh::vertices>(unewa);
auto u = m.mdcolex<mesh::vertices>(ua);
const auto dx_over_dy = m.xdelta() / m.ydelta();
const auto dy_over_dx = m.ydelta() / m.xdelta();
const auto dxdy = m.dxdy();
// boundary conditions (Dirichlet)
auto xverts = m.vertices<mesh::x_axis, mesh::all>();
const auto is = *xverts.begin();
const auto ie = *(xverts.end() - 1);
if (m.is_boundary<mesh::x_axis, mesh::boundary::low>(is)) {
for (auto j : m.vertices<mesh::y_axis, mesh::extended>()) {
u(is, j) = 0.;
}
}
if (m.is_boundary<mesh::x_axis, mesh::boundary::high>(ie)) {
for (auto j : m.vertices<mesh::y_axis, mesh::extended>()) {
u(ie, j) = 0.;
}
}
auto yverts = m.vertices<mesh::y_axis, mesh::all>();
const auto js = *yverts.begin();
const auto je = *(yverts.end() - 1);
if (m.is_boundary<mesh::y_axis, mesh::boundary::low>(js)) {
for (auto i : m.vertices<mesh::x_axis, mesh::extended>()) {
u(i, js) = 0.;
}
}
if (m.is_boundary<mesh::y_axis, mesh::boundary::high>(je)) {
for (auto i : m.vertices<mesh::x_axis, mesh::extended>()) {
u(i, je) = 0.;
}
}
for (auto j : m.vertices<mesh::y_axis>()) {
for (auto i : m.vertices<mesh::x_axis>()) {
unew(i, j) =
c * (1. / dxdy) *
(dy_over_dx * (u(i + 1, j) - 2 * u(i, j) + u(i - 1, j)) +
dx_over_dy * (u(i, j + 1) - 2 * u(i, j) + u(i, j - 1)));
}
}
}
The task is wrapped in heat_op, which implements the flecsolve operator
interface. Its apply method launches the task with the topology and
field references obtained from the input and output vectors:
struct heat_params {
double diffusivity;
};
struct heat_op : flecsolve::op::base<heat_params> {
explicit heat_op(double d) : flecsolve::op::base<heat_params>(d) {}
template<class Domain, class Range>
void apply(const Domain & x, Range & y) const {
flecsi::execute<task::laplace>(
y.data.topo(), params.diffusivity, y.data.ref(), x.data.ref());
}
};
Explicit integration¶
The explicit driver constructs an RK23 integrator using the heat operator,
the settings read from explicit.cfg, and work vectors created from the
solution vector:
void time_integration(control_policy & cp) {
flog(info) << "Time integration" << std::endl;
flecsi::flog::flush();
auto & u = cp.u();
auto & unew = cp.unew();
using namespace flecsolve;
using namespace flecsolve::time_integrator;
op::core<heat_op> F(cp.diffusivity);
rk23::integrator ti(rk23::parameters(
read_config("explicit.cfg", rk23::options("time-integrator")),
op::ref(F),
rk23::make_work(u)));
auto output = [&]() {
if (output_steps.value()) {
std::string fname{"timestep" +
std::to_string(ti.get_current_step())};
flecsi::execute<task::output, flecsi::mpi>(
flecsi::exec::on, cp.mesh(), u.data.ref(), fname.c_str());
}
};
output();
auto dt = ti.get_current_dt();
while (ti.get_current_time() < ti.get_final_time()) {
ti.advance(dt, u, unew);
auto good_solution = ti.check_solution();
if (good_solution) {
flog(info) << "Step " << ti.get_current_step() << " advanced " << dt
<< "s to time " << ti.get_current_time() << "s"
<< std::endl;
ti.update();
output();
std::swap(u, unew);
}
dt = ti.get_next_dt(good_solution);
}
}
}
At each iteration, advance proposes a new solution. The driver checks
the result, updates the integrator state, swaps the current and next vectors,
and asks the integrator for the next time step.
Implicit integration¶
The implicit driver uses BDF and a Krylov solver. An operator_adapter
turns the heat operator F = alpha Laplacian into the operator required by
the implicit method, such as I - gamma F. The solver factory selects the
linear solver from implicit.cfg:
namespace heat {
void time_integration(control_policy & cp) {
flog(info) << "Time integration" << std::endl;
auto & u = cp.u();
auto & unew = cp.unew();
using namespace flecsolve;
using namespace flecsolve::time_integrator;
auto [ti_settings, slv_settings] =
read_config("implicit.cfg",
bdf::options("time-integrator"),
krylov_factory::options("linear-solver"));
auto F = op::make_shared<operator_adapter<heat_op>>(cp.diffusivity);
bdf::integrator ti(
bdf::parameters(ti_settings,
F,
bdf::make_work(u),
krylov_factory::make_shared(slv_settings, u, F)));
auto output = [&]() {
if (output_steps.value()) {
std::string fname{"timestep" +
std::to_string(ti.get_current_step())};
flecsi::execute<task::output, flecsi::mpi>(
flecsi::exec::on, cp.mesh(), u.data.ref(), fname.c_str());
}
};
output();
auto dt = ti.get_current_dt();
bool first_step = true;
while (ti.get_current_time() < ti.get_final_time()) {
ti.advance(dt, first_step, u, unew);
auto good_solution = ti.check_solution();
if (good_solution) {
flog(info) << "Step " << ti.get_current_step() << " advanced " << dt
<< "s to time " << ti.get_current_time() << "s"
<< std::endl;
ti.update();
std::swap(u, unew);
first_step = false;
output();
}
dt = ti.get_next_dt(good_solution);
}
}
Next steps¶
To experiment with the example, try changing the following:
Set
diffusivitywith-dto change the rate of diffusion.Change
initial-dt,max-dt, or the final time in the configuration files.Increase the mesh extents to study spatial resolution and parallel scaling.
Modify
task::icsinheat.ccto use a different initial condition.Modify
task::laplaceinheat.hhto experiment with another stencil or boundary condition.