Simulating Fluids, Fireplace, and Smoke in ActualTime
Notes on the maths, algorithms, and strategies concerned in simulating fluids like hearth and smoke in realtime.
Supply code for this text may be discovered on my GitHub.
Fireplace is an attentiongrabbing graphics drawback. Previous approaches typically faked it. For instance,
Lord of the Rings
(the fluid sim was too costly on the time, even for motion pictures).
Actualtime purposes like video video games have just about completely
used nonphysical approaches.
However within the final 10 years GPUs have made quick fluid simulation straightforward.
Primary fluid dynamics algorithms are simple to implement on the GPU
Harry Potter
And in 2014, NVIDIA launched FlameWorks, a complete system for producing hearth and smoke results for video games.
This submit takes notes on how hearth may be simulated on the GPU. It walks by means of the maths behind fluid dynamics,
parallel algorithms for modeling fluids, and the additional combustion bits that make hearth particular. Readers ought to have an affordable
background in vector calculus and differential equations (know how one can take the gradient of a vector). Demos are carried out with WebGL.
1. Fluid Simulation
Earlier than we simulate hearth, we have to simulate fluid. We assume
our fluid is incompressible and
inviscid, which can vastly simplify our drawback.
1.1 Primary Fluid Dynamics
Suppose (D) is a area in area full of a fluid. At any level ( mathbf{x} in D ) and time (t), the fluid has velocity (mathbf{u}(mathbf{x}, t)).
Computationally, we are able to characterize the 2D velocity area ( mathbf{u} ) with an ( N instances N ) grid, the place the equally spaced grid factors give the worth of
the rate area at that time in area.
Ex: A 16(instances)16 grid representing ( mathbf{u} = (y, x) )
What is going to occur if we put a drop of dye within the fluid?
Let’s outline a scalar area ( psi (mathbf{x}, t) ) because the density of the dye at any level in area and time.
The transport of portions like ( psi ) inside a fluid by the fluid’s velocity is named advection.
Given some fluid’s velocity area and an preliminary density area of our dye, we would prefer to see how the dye’s density in all places
evolves over time by simulating its advection by means of the fluid.
A Naive Technique for Advection
One thought
the advection is to take every grid level, transfer ahead the course and distance that will be traveled
by a particle on the grid level’s velocity and the simulation timestep ( Delta t ), and replace the grid level nearest to the place the particle lands:
$$
psi (mathbf{x} + mathbf{u} (mathbf{x}, t), t + Delta t) = psi (mathbf{x}, t)
$$
However that is tough to parallelize, since 2 grid factors can find yourself in the identical goal level after ahead analysis.
And in observe, the goal level will fall between grid factors, that means it must be interpolated into the encompassing grid factors. Lastly, this course of is
unstable for time steps above some quantity, inflicting ( psi ) to explode.
The Advection Partial Differential Equations
This complete time we have been fixing a partial differential equation! If we will derive a steady methodology for advection, we’ll must
first get an express expression for this PDE. Let’s begin from first rules.
Think about a hard and fast area of area (W) (that’s, (W) doesn’t fluctuate with time). The overall mass of dye inside (W)
is ( int_{W} psi dV). Over time, the change in mass is:
$$
frac{d}{dt} int_{W} psi (mathbf{x}, t) dV = int_{W} frac{partial}{partial t} psi (mathbf{x}, t) dV
$$
Now, letting (S) denote the floor of (W) and ( mathbf{n} ) the outward regular vector outlined alongside the floor, we are able to study
the mass circulate price by means of the floor of (W). Particularly, observe that the amount circulate price – the quantity of fluid that
flows by means of per second – throughout (S) per unit space is ( mathbf{u} cdot mathbf{n} ) and the mass circulate price per unit space
is ( psi mathbf{u} cdot mathbf{n} ).
This provides us the regulation of conservation of mass in integral type:
$$
frac{d}{dt} int_{W} psi dV = – int_{S} psi mathbf{u} cdot mathbf{n} dA
$$
Can we do away with the integrals and say one thing comparable for factors? By divergence theorem, the above is equal to
$$
int_{W} [frac{partial psi}{partial t} + nabla cdot (psi mathbf{u})] dV = 0
$$
Then for a unit subregion ( W = dV ), we are able to say that
$$
frac{partial psi}{partial t} + nabla cdot (psi mathbf{u}) = 0
$$
This provides us an express PDE that we have to resolve for ( psi )!
Hmm… we might cease right here, however we would be capable to simplify this extra. Particularly, it appears to be like like we might isolate out
a time period ( nabla cdot mathbf{u} ) that goes to zero due to incompressibility.
$$
start{aligned}
&frac{partial psi}{partial t} = – nabla cdot (psi mathbf{u})
&= – (frac{partial}{partial x} psi u + frac{partial}{partial y} psi v)
&= – (frac{partial psi}{partial x} u + frac{partial u}{partial x} psi + frac{partial psi}{partial y} v + frac{partial v}{partial y} psi)
&= – (psi nabla cdot mathbf{u} + mathbf{u} cdot nabla psi)
&= – mathbf{u} cdot nabla psi
finish{aligned}
$$
Making use of our incompressibility constraint ( nabla cdot mathbf{u} = 0 ) on the finish yields a scalar PDE, the primary of our incompressible circulate advection equations:
$$
frac{partial psi}{partial t} = textual content{advection} (mathbf{u}, psi) = – mathbf{u} cdot nabla psi tag{1}
$$
$$
frac{partial mathbf{v}}{partial t} = textual content{advection} (mathbf{u}, mathbf{v}) = – mathbf{u} cdot nabla mathbf{v} tag{2}
$$
Eq. 2 for advecting a vector area ( mathbf{v} ) by means of our velocity area may be derived equally to the scalar advection equation.
A Steady Technique for Advection
Let’s look carefully at eqn. (1):
$$
frac{partial psi}{partial t} = – mathbf{u} cdot nabla psi
$$
Discover that the righthand time period is a directional byproduct within the ( mathbf{u} ) course. This provides us a beautiful new methodology for
advecting ( psi ) by an incompressible fluid – beginning at a grid level ( mathbf{x} ), hint the fluid velocity backwards,
changing the worth at our unique level by the worth that we land on (if we land between factors, we are able to interpolate):
$$
psi (mathbf{x}, t + Delta t) = psi (mathbf{x} – mathbf{u} (mathbf{x}, t), t)
$$
In GPU pseudocode:
world Vec2Field u;
world FloatField density;
world float dt;
// Run for every level in our scalar grid that we wish to replace
float advectPoint(vec2 x) {
vec2 coord = x  dt * getVec2At(u, x);
return getFloatAt(density, coord);
}
This methodology is named
SemiLagrangian advection and was invented in 1999 by Jos Stam
Like Euler, it is firstorder correct, however has precisely the extra properties we’d like:
 It is extraordinarily straightforward to parallelize as a result of every grid level solely will get up to date as soon as per iteration.

It is unconditionally steady. Why? Observe that for any grid level, the utmost worth it could get up to date to is the utmost worth
of all of the grid factors.
For a hard and fast velocity area fulfilling the incompressibility constraint, it really works nice.
Click on anyplace above to drop some dye within the circulate
1.2 The NavierStokes Equations
Up to now we have discovered a mannequin that describes how scalar properties of a fluid evolve over time, assuming the circulate is mounted. What in regards to the fluid circulate itself –
how does the rate area ( mathbf{u} ) have an effect on itself over time?
The NavierStokes equations
and smoothness drawback, one of many Clay Institute’s seven Millenium Prize issues in math.
$$
frac{partial mathbf{u}}{partial t} =
{underbrace{ – mathbf{u} cdot nabla mathbf{u} }_text{selfadvection}} +
{underbrace{ mu^2 nabla mathbf{u} }_text{diffusion}} –
{underbrace{ nabla p }_text{stress}} +
{underbrace{ textbf{F} }_text{ext. forces}}
$$
$$
textual content{the place~}forall t{~,~} nabla cdot mathbf{u} = 0
$$
Right here, the fixed ( mu ) is the fluid’s viscosity and ( mathbf{F} ) are
exterior forces. However we assumed earlier that our fluid was inviscid,
so ( mu = 0 ), and we are able to simply ignore exterior forces for now. So we’re left with two phrases – selfadvection and stress.
$$
frac{partial mathbf{u}}{partial t} =
{underbrace{ – mathbf{u} cdot nabla mathbf{u} }_text{selfadvection}} –
{underbrace{ nabla p }_text{stress}}
tag{3}
$$
$$
textual content{the place~}forall t{~,~} nabla cdot mathbf{u} = 0
$$
If at each timestep we numerically compute these phrases and add them, we are able to simulate our fluid! In pseudocode:
let u = createVectorGrid();
let density = createScalarGrid();
let p = createScalarGrid();
whereas (true) {
// Clear up for the following velocity area.
u = advect(u, u);
p = computePressure(...);
u = u  gradient(p);
// Advect dye by means of the brand new velocity area.
density = advect(u, density);
}
Let’s take a more indepth have a look at every of those.
SelfAdvection
From our incompressible advection equations, we are able to see that the selfadvection time period is the advection of the fluid’s velocity
area ( mathbf{u} ) by means of itself:
$$
textual content{advection} (mathbf{u}, mathbf{u}) = – mathbf{u} cdot nabla mathbf{u} tag{4}
$$
The place do the opposite phrases come from? Properly, advecting ( mathbf{u} ) by means of itself yields a brand new velocity area
( mathbf{u}^prime ) (computable by way of the SemiLagrangian backtracing algorithm from above):
$$
mathbf{u}^prime = mathbf{u} – mathbf{u} cdot nabla mathbf{u}
$$
Strain
We do not know if this new velocity area follows the incompressibility constraint (e.g. has zero divergence).
So the stress time period ( p ) must appropriate this someway:
$$
nabla cdot (mathbf{u}^prime – nabla p) = 0
$$
We rearrange this to get
$$
nabla^2 p = nabla cdot mathbf{u}^prime tag{5}
$$
which is a kind of equation referred to as a Poisson equation, the place the lefthand facet is
the Laplacian of an unknown scalar area and the
righthand facet is a recognized scalar. Fixing this Poisson equation is basically the slowest computational step
in fluid simulation, for causes we’ll see shortly.
Fixing for Strain
So how will we resolve this specific PDE for ( p )? Properly, we all know the worth of our candidate velocity area ( mathbf{u}^prime )
in any respect of our grid factors, so we are able to roughly compute the righthand facet of the Poisson equation by making use of a discrete
model
$$
nabla cdot mathbf{u}^prime approx
frac{ u_{i+1, j} – u_{i1, j} }{ 2 } +
frac{ v_{i, j+1} – v_{i, j1} }{ 2 }
$$
the place ( mathbf{u}^prime = (u, v) ) in 2 dimensions.
Then we are able to use a discrete model of the Laplacian
$$
nabla^2 p approx
p_{i+1, j} + p_{i1, j} + p_{i, j+1} + p_{i, j1} – 4p_{i, j}
$$
to rework the entire equation right into a linear equation with 5 unknowns.
However actually, we’re fixing the Poisson equation (5) over all of area without delay, so for an ( N instances N ) grid,
we find yourself with a system of ( N^2 ) linear equations with precisely ( N^2 ) unknowns! So we find yourself with
the acquainted previous equation
$$
mathbf{Ax} = mathbf{b}
$$
the place ( mathbf{A} ) is a matrix making use of the Laplacian operator to the entire grid and ( mathbf{b} ) is a vector containing
the rate area’s divergence in any respect grid factors.
There are numerous offtheshelf algorithms for fixing linear techniques precisely. Sadly for us, even the quickest algorithms
scale superlinearly with our grid dimension.
Fixing for Strain… Effectively
If we will make a realtime simulation, we have to go quick. Can we leverage the GPU?
Properly, it is not likely doable to attain an actual resolution to the linear system effectively, however we must always
notice that the linear system is already an approximation to the Poisson equation. And it’s doable to attain
arbitrarily correct approximations with iterative strategies – which start with an estimate and enhance resolution
accuracy each iteration – so we are able to simply choose an iterative algorithm and run it
till we now have one thing that is “ok”. One significantly easy and easytoimplement iterative algorithm for fixing
linear equations is the Jacobi method.
We begin with the very first equation within the system:
$$
A_{11}x_1 + A_{12}x_2 + … + A_{1n}x_n = b_1
$$
On the (okay)th iteration, given some guess ( mathbf{x}^okay ) for the answer ( mathbf{x} ), we now have some error.
We will use this error to replace our guess for ( x_1 ) as follows:
$$
x_1^{okay+1} = frac{ b_1 – A_{12}x_2^okay – … – A_{1n}x_n^okay }{ A_{11} }
$$
In Jacobi, our guesses for all parts of ( mathbf{x} ) are executed in parallel, giving an ideal
match for implementation on the GPU. In pseudocode:
world FloatField divergence;
world FloatField stress;
world float texelSize;
// Run for every level in our stress grid that we wish to replace
float iterateJacobi(vec2 x) {
float div = getFloatAt(divergence, x);
float L = getFloatAt(stress, x + vec2(texelSize, 0.));
float R = getFloatAt(stress, x + vec2(texelSize, 0.));
float T = getFloatAt(stress, x + vec2(0., texelSize));
float B = getFloatAt(stress, x + vec2(0., texelSize));
return (div  L  R  T  B) / 4.;
}
It is price noting that different, fasterconverging solvers can be
carried out on the GPU, just like the Conjugate Gradient methodology and the Multigrid methodology. However relying on
the fluid and utility, stress accuracy is probably not as necessary as advection accuracy or ease of implementation.
For smoke and hearth, modifications in fluid quantity aren’t as obvious as they’re for fluids like water
Abstract: Simulating Navier Stokes
The maths behind NavierStokes generally is a little bit dense, however at a highlevel, simulating a fluid by fixing the equations
comes down to some key replace procedures on a grid per timestep. For our dye drawback, this is our simulation
would possibly look:
let u = createVectorGrid();
let density = createScalarGrid();
let div = createScalarGrid();
let p = createScalarGrid();
whereas (true) {
// Clear up for the following velocity area.
u = advect(u, u);
// Implement incompressibility with stress projection.
div = divergence(u);
for (let i = 0; i < JACOBI_ITERATIONS; i++) {
p = updatePressure(p, div);
}
u = u  gradient(p);
// Advect dye by means of the brand new velocity area.
density = advect(u, density);
}
1.3 Vorticity Confinement
Utilizing a grid to retailer our velocity area is extraordinarily handy, but it surely ends in undesirable
numerical smoothing each time we now have to interpolate values between grid factors.
This mixed with the comparatively coarse approximation of a firstorder SemiLagrangian advection
scheme has the impact of dissipating out turbulent vortices in our circulate. Bodily, the rate area
loses vitality, and the top result’s typically overly clean, “boring” fluid circulate.
One technique to fight misplaced vorticity is to extend the decision of our grid, however this is not actually possible
for realtime simulations which have restricted computational assets. What we’d ideally
love to do is use all of the small particulars that get smoothed over every step of the simulation, and amplify
them. This course of is named vorticity confinement – admittedly, it isn’t completely sensible, however succeeds in
preserving small scale particulars in kind of bodily appropriate places
Certainly, it was initially invented to resolve very advanced circulate fields in engineering simulations of helicopter blades,
the place it simply wasn’t doable so as to add the variety of vital grid factors
The smallest turbulent options we are able to discover are the vortices centered at every grid level in our simulation.
We will measure the depth of those vortices (the vorticity of them) by taking the curl of ( mathbf{u} )
at every level, and amplify them by basically including a round circulate scaled by vorticity about every level.
Mathematically, the vorticity is outlined by
$$
bm{omega} = nabla instances mathbf{u}
$$
For every grid level, we compute a normalized location vector that factors to the best close by vorticity focus:
$$
mathbf{N} = frac{ nabla  bm{omega}  }{  nabla  bm{omega}   }
$$
And eventually, we compute the confined vorticity vector area and add it to our circulate:
$$
mathbf{f_{conf}} = epsilon (mathbf{N} instances bm{omega})
$$
$$
mathbf{u_{conf}} = mathbf{u} + mathbf{f_{conf}}
$$
Right here, the confinement fixed ( epsilon > 0 ) is a parameter controlling the quantity of small scale element added
again to the circulate. Even low confinement ranges (round 015) could make an enormous distinction, particularly for simulations
utilizing SemiLagrangian advection schemes, and better confinement ranges can create extremely stylized, billowing flows.
Click on and drag to drop some dye within the turbulent simulation above
On the GPU, we are able to compute curl and confinement like so:
world Vec2Field u;
world float texelSize;
// Run to get curl for every level in grid
float computeCurl(vec2 x) {
float L = getVec2At(u, x + vec2(texelSize, 0.)).y;
float R = getVec2At(u, x + vec2(texelSize, 0.)).y;
float T = getVec2At(u, x + vec2(0., texelSize)).x;
float B = getVec2At(u, x + vec2(0., texelSize)).x;
return (R  L)  (T  B);
}
world Vec2Field curl;
world float confinement;
// Run to get confinement power for every level in grid
vec2 confinementForce(vec2 x) {
float L = getFloatAt(curl, x + vec2(texelSize, 0.));
float R = getFloatAt(curl, x + vec2(texelSize, 0.));
float T = getFloatAt(curl, x + vec2(0., texelSize));
float B = getFloatAt(curl, x + vec2(0., texelSize));
float C = getFloatAt(curl, x);
vec2 N = vec2(abs(T)  abs(B), abs(R)  abs(L));
N = N / size(N);
return confinement * C;
}
The complete simulation with turbulence:
let u = createVectorGrid();
let density = createScalarGrid();
let div = createScalarGrid();
let p = createScalarGrid();
let curl = createVectorGrid();
whereas (true) {
// Clear up for the following velocity area.
u = advect(u, u);
// Use vorticity confinement to amplify turbulence of velocity area.
curl = computeCurl(u);
u = u + confinementForce(curl, CONFINEMENT);
// Implement incompressibility with stress projection.
div = divergence(u);
for (let i = 0; i < JACOBI_ITERATIONS; i++) {
p = updatePressure(p, div);
}
u = u  gradient(p);
// Advect dye by means of the brand new velocity area.
density = advect(u, density);
}
CurlNoise Turbulence
Curl noise is a technique that basically does the identical factor as vorticity confinement, however as a substitute of measuring
and amplifying the vorticity of the rate area, a scalar vorticity area is produced from scratch utilizing noise capabilities.
Mathematically, we are able to mix vorticity confinement and curlnoise turbulence by synthesizing a random
vorticity area
$$
bm{phi} = textual content{rand} * mathbf{z}
$$
Then computing our remaining vorticity area by
$$
bm{omega}^* = bm{omega} + bm{phi}
$$
Quickmoving, highlyturbulent fluids like smoke and hearth profit essentially the most from vorticity confinement and curl noise,
and in observe the curl noise area ( bm{phi} ) each evolves with time and can also be advected by the fluid circulate.
2. Fireplace Simulation
Should you’ve gotten this far, pat your self on the again! The strategies within the earlier part allow us to effectively
and precisely simulate fluids with various bodily parameters (oil, water, honey)
assuming the fluid area is a hard and fast area. These inquisitive about dealing with various domains (that’s, fluids that occupy totally different areas
throughout the grid, like a halffull cup of water that sloshes round) will wish to discover accounting for various boundary situations
throughout the grid simulation
GPU Gems 3 Chapter 30.2
There are additionally non gridbased strategies, however these are exterior scope right here.
Simulating hearth and smoke requires a pair additions. First, we’ll want so as to add channels representing
gasoline and temperature to our simulation, and mannequin the combustion of gasoline to create warmth. Subsequent we’ll deal with how sizzling pockets
of our fluid rise with a thermal buoyancy mannequin, and eventually, we’ll must render our flames appropriately, bearing in mind
blackbody radiation of the flames, human notion of sunshine, and hearth motion.
2.1 A Primary Combustion Mannequin
Chemically, hearth is attributable to the oxidation of a gasoline materials in a response that releases each warmth and lightweight.
In our case, we are able to assume that any gasoline in our system has already ignited and is actively including warmth; we cannot fear
about the issue of unignited gasoline.
To be extra particular, let’s outline a scalar area ( rho ) the place ( 0 leq rho leq 1 ) represents the density of gasoline and
one other scalar area ( T > 0 ) representing the temperature of the fluid in all places. At each timestep, temperature is
added to the system by the gasoline, which burns at a given burn temperature:
$$
T^prime = textual content{max} ( T, rho * T_{textual content{burn}} )
$$
In fact, temperature is not static – warmth diffuses from sizzling to chilly areas, and with fluids particularly, largescale actions
of molecules transport warmth. The mixture of those 2 processes defines heat convection,
and conveniently, we have already got a mathematical mannequin for the way it works – advection! Simulationwise, we advect our temperature
area alongside our velocity area. Since any reacting molecules are additionally moved by the fluid, we must always advect gasoline as effectively.
The warmth itself additionally impacts the motion of the fluid – we’ll see how one can deal with this shortly.
Moreover, sizzling molecules radiate off temperature as mild
to it when rendering the fireplace coloration. The soot particles current in most hearth radiate like supreme blackbodies.
StefanBoltzmann Law,
in a quintic equation
$$
T^prime = T – sigma_{textual content{cool}} ( frac{T}{T_{textual content{max}}} )^4 * Delta t
$$
Right here, ( sigma_{textual content{cool}} ) is the cooling price parameter. For a bodily appropriate simulation, we’d set it to the StefanBoltzmann fixed,
however for a graphical simulation, it is good for the artist to have the ability to management the speed of cooling.
To finish our combustion mannequin, notice that our gasoline is at all times burning (we are able to think about it because the density of ionized fuel particles
that give off thermal vitality and return to a decrease vitality state), so each timestep we dissipate it by some given burn
price ( gamma_{gasoline} ):
$$
rho^prime = rho (1 – gamma_{gasoline})^{Delta t}
$$
2.2 Thermal Buoyancy
Up to now, our temperature area would not do something to our fluid circulate. Nevertheless it ought to – sizzling pockets of air ought to increase and rise, and cooler pockets ought to fall.
We will mannequin this with a thermal buoyancy power. Since we’re assuming incompressibility, we cannot truly deal with air enlargement, however the fluid circulate
ought to expertise an upward power relying on temperature:
$$
mathbf{u}^prime = mathbf{u} + (beta T Delta t) mathbf{j}
$$
Right here, ( beta ) is a given optimistic buoyancy fixed, and ( mathbf{j} ) is the upward unit vector.
Including a combustion mannequin and thermal buoyancy power offers us a implausible simulator for a decidedly “firelike” fluid –
with the precise values of buoyancy and cooling, we are able to get cumbersome, billowing plumes of fabric.
Not precisely flames, however similar to smoke.
Faucet and drag within the simulation beneath to inject some combusting gasoline.
The displayed pixels characterize density of smoke particles, which dissipate at a relentless price as a substitute of getting used
up throughout combustion, however are nonetheless advected by the fluid simulation.
Click on and drag so as to add smoke above
The simulation code builds off the fundamental fluid routines:
let u = createVectorGrid();
let density = createScalarGrid();
let div = createScalarGrid();
let p = createScalarGrid();
let curl = createVectorGrid();
let gasoline = createScalarGrid();
let temp = createScalarGrid();
whereas (true) {
// Clear up for the following velocity area.
u = advect(u, u);
// Combustion step.
temp = combust(temp, gasoline);
// Use vorticity confinement to amplify turbulence of velocity area.
curl = computeCurl(u);
u = u + confineVorticity(curl, CONFINEMENT);
// Add thermal buoyancy.
u = u + buoyancy(temp);
// Implement incompressibility with stress projection.
div = divergence(u);
for (let i = 0; i < JACOBI_ITERATIONS; i++) {
p = updatePressure(p, div);
}
u = u  gradient(p);
// Advect dye by means of the brand new velocity area.
density = advect(u, density);
}
2.3 Fireplace Rendering
Fireplace is a participating medium,
that means it emits mild by means of blackbody radiation
5.1 in Nguyen et al. 2002
rendering our combusting gasoline simulation utilizing the proper formulation is then all we have to go from smoke to fireplace!
The temperaturetocolor spectrum as described by Planck’s regulation
Planck’s Law describes the spectral density of sunshine radiated by
a black physique at a given temperature ( T ):
$$
M(lambda, T) = frac{ c_1 }{ lambda^5 } frac{ 1 }{ exp{ frac{c_2}{lambda T} } – 1 }
$$
the place
$$
c_1 = 2 pi h c^2
c_2 = frac{hc}{okay}
$$
and (h), (c), and (okay) are Planck’s fixed, the pace of sunshine, and Boltzmann’s fixed, respectively.
After implementing blackbody rendering utilizing fragment shaders, we now have an entire hearth simulation!
Click on and drag so as to add hearth above
That is it for these notes! There’s in fact rather more to fluid and hearth simulation not lined right here, like
totally different (e.g. nongridbased) methods for fixing the identical drawback of simulation inside a hard and fast quantity,
totally different issues to unravel involving various domains or dynamic obstacles, enhancements to rendering like
extra correct blackbody radiation, mild scattering, or postprocessing results. Useful introductions
to those subjects may be present in references beneath.