Boundary Conditions Implementation
Absorbing boundary conditions dramatically improve the fidelity of our simulation.
Absorbing boundary conditions are provided by new update equations applied at the boundary of the grid. We have one equation for the right edge at .
And another equation for the left edge at .
Another challenge is that terms appear on both sides of these equations.
The update equations look like slightly different equations, but they are really special cases of a more general formula.
where
We implement these update equations in a MurBoundary class. This separate class isolates the code for the boundary conditions, and makes it easier to evolve and adapt the code to other simulations. We incorporate the MurBoundary into our process by computing the update for the full domain, then invoke the boundary computations to overwrite the boundary values.
We compute boundary values only on the boundary. For a 1D simulation such as ours, the boundary is a pair of points, one at each end of the domain. To render these points we start with a two element workgroup.
@compute @workgroup_size(2)
fn recomputeBoundary(@builtin(global_invocation_id) global_id : vec3u)
Then the invocation with global_id.x
= 0 computes the left edge boundary, and
the invocation with global_id.x
= 1 computes the right edge boundary.
// Index will be 0 at the left edge, and 1 at the right edge.
let index = global_id.x;
let sotrageBufferIndex = index*(parameters.xResolution-1);
let offset = 1 - 2*index;
The values for at and are readily available from the storage arrays. We also pass in the parameters from the simulation. The only new field is the phase velocity.
struct Parameters {
dt: f32, // The time step, Δt.
xResolution: u32, // The number of points along the x-axis, the number of elements in the array.
length: f32, // The physical length for our simulation.
potential: array<f32> // The potential the particle moves through.
}
// group 0, things that never change within a simulation.
// The parameters for the simulation
@group(0) @binding(0) var<storage, read> parameters: Parameters;
// Group 1, changes on each iteration - the same as in the main solver to keep the same bindings.
// Older wave function at t-Δt.
@group(1) @binding(0) var<storage, read> oldWaveFunction : array<vec2f>;
// Current wave function at t.
@group(1) @binding(1) var<storage, read> waveFunction : array<vec2f>;
// The updated wave function at t+Δt.
@group(1) @binding(2) var<storage, read_write> updatedWaveFunction : array<vec2f>;
// Group 2, boundary value specific data.
@group(2) @binding(0) var<uniform> phaseVelocity : f32;
@compute @workgroup_size(2)
fn recomputeBoundary(@builtin(global_invocation_id) global_id : vec3u)
{
// Index will be 0 at the left edge, and 1 at the right edge.
let index = global_id.x;
let sotrageBufferIndex = index*(parameters.xResolution-1);
let offset = 1 - 2*index;
let dx = parameters.length / f32(parameters.xResolution-1);
updatedWaveFunction[sotrageBufferIndex] = waveFunction[sotrageBufferIndex+offset]
+ ((phaseVelocity*parameters.dt-dx)/(phaseVelocity*parameters.dt+dx))
*(updatedWaveFunction[sotrageBufferIndex+offset]-waveFunction[sotrageBufferIndex]);
}
These boundary conditions allow us to present a wider variety of simulations than was available previously.
