Heat Conduction Simulation Tutorial This tutorial walks through implementing heat conduction simulations, starting with basic vanilla C implementations and progressing to more powerful Basilisk C versions. We build from 1D steady-state problems to 2D, annular, and axisymmetric geometries, demonstrating how Basilisk’s domain-specific language dramatically simplifies code while adding powerful features like adaptive mesh refinement and embedded boundaries. The tutorial includes eight progressive simulation cases, showing how the same physics can be implemented with increasingly sophisticated numerical methods, and concludes with custom post-processing tools for data extraction and visualization. Basilisk ultimately reduces hundreds of lines of vanilla C code to concise, readable implementations that handle complex geometries with minimal effort.
Introduction
This guide walks you through implementing heat conduction simulations using both vanilla C and the Basilisk framework. We’ll progress from simple steady-state problems to more complex transient simulations in multiple dimensions, highlighting the advantages of Basilisk for computational fluid dynamics (CFD) and heat transfer applications.
Learning Objectives
Implement basic heat conduction solvers in vanilla C
Transition to equivalent implementations in Basilisk C
Understand the advantages of Basilisk’s domain-specific language
Develop simulations with increasing complexity (1D → 2D → axisymmetric)
Appreciate how Basilisk simplifies numerical methods implementation
Learn to analyze and visualize simulation results with post-processing tools
Throughout this tutorial, we’ll follow a logical progression across eight core simulation cases:
Steady-state heat conduction in vanilla C (1-conduction-simple.c)
Transient heat conduction in vanilla C (1-conduction-transient.c)
Steady-state heat conduction in Basilisk (1-conduction-simple-basilisk.c)
Transient heat conduction in Basilisk (fill-in exercise) (1-conduction-transient-basilisk.c)
Enhanced transient solution with diffusion module (1-conduction-transient-basilisk-withHeaders.c)
2D heat conduction (1-conduction-2D.c)
Heat conduction in an annulus (1-conduction-2D-annulus.c)
Axisymmetric heat conduction (1-conduction-Axi.c)
After completing the simulations, we’ll also cover how to post-process and visualize the results using custom tools.
Let’s implement the transient heat equation in Basilisk. Here’s a template with some parts for you to fill in:
#include "grid/cartesian1D.h"#include "run.h"// Declare scalar field for temperaturescalar T[];// Simulation parameters#define EPS 0.1 // Width of initial temperature peak#define tmax 1.0#define tsnap 0.1int main() { // Domain setup L0 = 10.0; // Domain length X0 = -L0/2; // Left boundary N = 200; // Number of cells // Set the timestep based on stability criterion (CFL condition) // dt = dx^2/2 for explicit scheme DT = (L0/N)*(L0/N)/2; // Create output directory char comm[80]; sprintf (comm, "mkdir -p intermediate"); system(comm); // Run simulation run();}/ * Initialize temperature field * * Sets up a "Dirac delta" approximated by a thin rectangle * centered at x=0 with total integral = 1. */event init (t = 0) { foreach() T[] = (fabs(x) < EPS) ? 1.0/EPS/2.0 : 0.0; }/ * Time integration using explicit finite volume method */event integration (i++) { // Get timestep for this iteration double dt = dtnext(DT); // Compute fluxes at the faces of the cells // q_i = -(T_i - T_{i-1})/Delta scalar q[]; foreach() q[] = -(T[] - T[-1])/Delta; // Update temperature // dT/dt = -(q_{i+1} - q_i)/Delta // this enforces the central difference scheme scalar dT[]; foreach() dT[] = -(q[1] - q[])/Delta; // Explicit Euler step foreach() T[] += dt*dT[];}/ * Save snapshots at regular intervals */event writingFiles (t += tsnap; t < tmax+tsnap) { char filename[100]; sprintf(filename, "intermediate/snapshot-%5.4f.csv", t); FILE *fp = fopen(filename, "w"); foreach() { fprintf(fp, "%g,%g\n", x, T[]); } fclose(fp);}/ * Save final results and comparison with analytical solution */event end (t = end) { char filename[100]; sprintf(filename, "conduction-transient.csv"); FILE *fp = fopen(filename, "w"); foreach() { fprintf(fp, "%g,%g\n", x, T[]); } fclose(fp);}
Stability Condition (K value) The stability condition for the 1D heat equation with explicit Euler time-stepping requires: So K = 2 is the appropriate value to ensure numerical stability.
Understanding The Heat Equation Discretization The heat equation can be discretized in two steps:
Compute heat fluxes at cell interfaces:
Update temperature based on flux divergence:
In Basilisk’s staggered grid, q[] is at cell faces while T[] is at cell centers.
Part 4: Enhanced Transient Solution with Diffusion Module
Now let’s look at an enhanced implementation using Basilisk’s built-in diffusion module.
#include "grid/multigrid1D.h" /* Multigrid solver is required by diffusion.h */#include "run.h"#include "diffusion.h"// Declare scalar field for temperaturescalar T[];// Boundary conditions// The diffusion solver will use homogeneous Neumann conditions by defaultT[left] = neumann(0.);T[right] = neumann(0.);// Simulation parameters#define EPS 0.1 // Width of initial temperature peak#define tmax 1.0#define tsnap 0.1int main() { // Domain setup L0 = 10.0; // Domain length X0 = -L0/2; // Left boundary N = 10000; // Number of cells // We can use a larger timestep with the implicit solver // compared to the explicit version which requires dt = dx^2/2 DT = 0.01; // Create output directory char comm[80]; sprintf (comm, "mkdir -p intermediate"); system(comm); // Run simulation run();}/ * Initialize temperature field */event init (t = 0) { foreach() T[] = (fabs(x) < EPS) ? 1.0/EPS/2.0 : 0.0; }/ * Time integration using implicit diffusion solver */event integration (i++) { // Get timestep for this iteration double dt = dtnext(DT); // Use the diffusion() function from diffusion.h to solve the equation diffusion(T, dt);}
Advantages of Using diffusion.h
Implicit time stepping allows for much larger timesteps
Unconditional stability removes the CFL restriction
Higher accuracy due to careful implementation of boundary conditions
Multigrid acceleration for faster convergence
Adaptive mesh refinement compatible
Handles complex geometries through embedded boundaries
Part 5: Advanced Heat Conduction Problems
5.1 2D Heat Conduction (1-conduction-2D.c)
Now, let’s extend our simulations to two dimensions.
Problem Statement Solve the 2D transient heat equation: With Dirichlet boundary conditions:
Customizing Visualizations You can adjust the visualization parameters to highlight different aspects of your results:
Change the colormap using the cmap parameter in ax.imshow()
Adjust the value range with vmin and vmax
Modify the plot layout or add additional subplots
Include contour lines with ax.contour()
Conclusion
Throughout this tutorial, we’ve progressed from simple 1D steady-state heat conduction problems in vanilla C to complex geometries and dimensions using Basilisk C. The key takeaways include:
Basilisk significantly reduces the complexity and verbosity of numerical simulation code
The event-based programming model provides a clean separation of concerns
Basilisk’s domain-specific language makes numerical methods more intuitive
Built-in modules like diffusion.h and embed.h provide powerful, optimized solvers
Extending to higher dimensions and complex geometries is straightforward
Post-processing tools are essential for analyzing and interpreting simulation results
Next Steps
Explore other physical phenomena like advection-diffusion
