Boundary Conditions Implementation

Absorbing boundary conditions dramatically improve the fidelity of our simulation.

A realistic simulation requires absorbing boundary conditions.

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 x = L .

Ψ L t + Δ t = Ψ L - Δ x t + v p Δ t - Δ x v p Δ t + Δ x ( Ψ L - Δ x t + Δ t - Ψ L t )

And another equation for the left edge at x = 0 .

Ψ 0 t + Δ t = Ψ Δ x t + v p Δ t - Δ x v p Δ t + Δ x ( Ψ Δ x t + Δ t - Ψ 0 t )

Another challange is that t + Δ t terms appear on both sides of these equations.

We impliment 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 the geometry. We use the geometry directly without any projections or other transformations, so we know that the points will be (-1, 0, 0) for the left edge, and (1, 0, 0) for the right edge. Shift the left edge in a bit to keep it from being trimmed by OpenGL.


  /**
   * Return 1D boundary with texture coordinates for Mur boundary conditions.
   * One point at each end of the 1xn strip we are using in this simulation.
   * (-1, 0) as the left boundary, (1,0) as the right edge.
   *
   * @returns {Float32Array} A set of points and textures containing one point at
   *                         each end of the simulation.
   */
  this.getGeometry = function ()
  {
    // Sets of x,y,z(=0),s,t coordinates.
    return new Float32Array([-1.0+(0.5/xResolution), 0.0, 0.0, 0.0, 0.0,  // left edge
                              1.0,                   0.0, 0.0, 1.0, 0.0]);// right edge
  };
      

The question now is how to write a shader that will evaluate these update equations at these points. The update equations look like slightly different equations, but they are really special cases of a more general formula.

Ψ X t + Δ t = Ψ X + ϵ Δ x t + v p Δ t - Δ x v p Δ t + Δ x ( Ψ X + ϵ Δ x t + Δ t - Ψ X t )

where

ϵ = { + 1 if   X = 0 ; - 1 if   X = L ;

The values for Ψ at t are supplied by texture reads. But what about the value at t + Δ t ? We will use the standard update equations to compute these values.

We just computed the t + Δ t , why not just read them back? These values are stored in the texture that we are writing, and we can not read and write the same texture.


  /**
   * Return 1D boundary with texture coordinates for Mur boundary conditions.
   * One point at each end of the 1xn strip we are using in this simulation.
   * (-1, 0) as the left boundary, (1,0) as the right edge.
   *
   * @returns {Float32Array} A set of points and textures containing one point at
   *                         each end of the simulation.
   */
  this.getGeometry = function ()
  {
    // Sets of x,y,z(=0),s,t coordinates.
    return new Float32Array([-1.0+(0.5/xResolution), 0.0, 0.0, 0.0, 0.0,  // left edge
                              1.0,                   0.0, 0.0, 1.0, 0.0]);// right edge
  };
      

The question now is how to write a shader that will evaluate these update equations at these points. The update equations look like slightly different equations, but they are really special cases of a more general formula.

Ψ X t + Δ t = Ψ X + ϵ Δ x t + v p Δ t - Δ x v p Δ t + Δ x ( Ψ X + ϵ Δ x t + Δ t - Ψ X t )

where

ϵ = { + 1 if   X = 0 ; - 1 if   X = L ;

The values for Ψ at t are supplied by texture reads. But what about the value at t + Δ t ? We will use the standard update equations to compute these values.

We just computed the t + Δ t , why not just read them back? These values are stored in the texture that we are writing, and we can not read and write the same texture.

The boundary conditions use the same Ψ on a grid at t and t - Δ t . It writes results to the same texture for Ψ at t + Δ t and the size and spacing on the grid are provided in the same way. The boundary condition has the same inputs, uniforms, as the Schrödinger solver with the addition of the phase velocity, vp.


  // The delta-t for each timestep.
  uniform float dt;
  // The physical length of the grid in nm.
  uniform float length;
  // An estimate of the phase velocity at the boundary
  uniform float vp;
  // At time t - delta t waveFunction.r is the real part waveFunction.g is the imaginary part.
  uniform sampler2D oldWaveFunction;
  // The number of points along the x axis.
  uniform int xResolution;
  
  // At time t waveFunction.r is the real part waveFunction.g is the imaginary part.
  uniform sampler2D waveFunction;
  // Discrete representation of the potential function.
  uniform sampler2D potential;
  
  // Vector to mix the real and imaginary parts in the wave function update.
  const vec2 mixing = vec2(-1.0, +1.0);
  
  varying vec2 vTextureCoord;
  
  void main()
  {
    float dx;
    vec2  ds;
    vec2  innerTextureCoord;
    vec2  innerValue;
    float offset;
    vec4  value;
  
    dx                = length/float(xResolution);
    ds                = vec2(1.2/float(xResolution), 0.0);
    // On the left edge, compute psi(x+dx), on the right edge, compute psi(x-dx)
    offset            = 1.0 - 2.0*step(0.5, vTextureCoord.s);
    innerTextureCoord = vTextureCoord + offset*ds;
    value             = texture2D(waveFunction, innerTextureCoord);
    // One step in from the edge, at t+Δt
    innerValue        =  texture2D(oldWaveFunction, innerTextureCoord).rg
                         + ((texture2D(waveFunction, innerTextureCoord+ds).gr
                               - 2.0*value.gr
                               + texture2D(waveFunction, innerTextureCoord-ds).gr)/(dx*dx)
                            - 2.0*texture2D(potential, innerTextureCoord).r*value.gr)*mixing*dt;
  
    gl_FragColor.rg = value.rg
                      + ((vp*dt-dx)/(vp*dt+dx))*(innerValue
                                                 - texture2D(waveFunction, vTextureCoord).rg);
  }
      

The innerValue is Ψ X + ϵ Δ x t + Δ t , which must be recomputed.

These boundary conditions allow us to present a wider variety of simulations than was available previously.

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License.