Verification

Verification, the process of ensuring that your project achieves its goals, is incredibly important for any software project, especially for one that presents itself as a quality learning example.

We now have two approaches to computing ψ(x,t). We have a closed form exact expression for a free particle wave function, and a finite difference method to compute its time evolution. Ideally, these produce identical results for the case where V=0. It will be insightful to actually make the comparison.

The FDTD is running while we collect error measurements below.
The exact solution provides a comparison with the FDTD solution.

Visually, these are almost, if not completely, identical. A numerical comparison will be much more precise and show some small differences. We will compare three quantities: ψr,ψi and ψ*ψ. Specifically, we compute the root-mean-square error for each of these.

E = 1N i = 0 N - 1 ( f i FDTD - f i exact ) 2

The table shows error data for every 200th frame from the animation. We use 501 FDTD steps between each frame, so this accounts for over 600,000 steps in the simulation. Over the course of this run, the ψ*ψ error quickly plateaus. However, the ψr,ψi errors both grow. This is actually expected and is an example of numerical dispersion.

The RMS error for ψr, ψi and ψ*ψ with 1600 grid points.
Frame Steps Eψr Eψi Eψ*ψ
0 0 0.0000 0.0000 3.1811e-18
200 100200 0.0021347 0.0021347 0.0021383
400 200400 0.0042693 0.0042693 0.0027217
600 300600 0.0064031 0.0064031 0.0026787
800 400800 0.0085367 0.0085367 0.0025088
1000 501000 0.010669 0.010669 0.0023335
1200 601200 0.012800 0.012800 0.0021778

As we would expect for such a numerical artifact, it decreases rapidly with a finer mesh. Doubling the number of mesh points decreases the error by an order or magnitude.

The RMS error for ψr, ψi and ψ*ψ with 3200 grid points.
Frame Steps Eψr Eψi Eψ*ψ
0 0 0.0000 0.0000 3.3140e-18
200 100200 0.00053378 0.00053367 0.00053275
400 200400 0.0010675 0.0010674 0.00067682
600 300600 0.0016010 0.0016010 0.00066536
800 400800 0.0021350 0.0021349 0.00062275
1000 501000 0.0026688 0.0026687 0.00057909
1200 601200 0.0032021 0.0032019 0.00054041

The Code

An exact ψ(x,t)

Our first step is to extend the initialization from the previous example into a more general class capable of generating the wave function at an arbitrary time.

This reworked shader takes time as an input, and provides the wave function as its output.


    // For wave function evaluation at arbitrary times
    @group(3) @binding(0) var<storage, read> time: f32;
    @group(3) @binding(1) var<storage, read_write> timeWaveFunction : array<vec2f>;
            

The new @compute annotation marks a second entry point to our shader, this time computing a single value for ψ(x,t), and loading it into the timeWaveFunctionBuffer.


    /*
     * Populate an array with the wave function for the given time.
     */
    @compute @workgroup_size(WORKGROUP_SIZE)
    fn computeTimeValues(@builtin(global_invocation_id) global_id : vec3u)
    {
      let index = global_id.x;
      // Skip invocations when work groups exceed the actual problem size
      if (index >= parameters.xResolution) {
        return;
      }
      timeWaveFunction[index] = computePsi(index, time);
    }
            

To execute the shader with this version of the code, create a second pipeline specifying the entry point as "computeTimeValues".


    timeEvolvedPipeline = device.createComputePipeline({
      layout: device.createPipelineLayout({
        bindGroupLayouts: [
          parametersBindGroupLayout, wavePacketParametersBindGroupLayout,
          waveFunctionBindGroupLayout, timeBindGroupLayout
        ]
      }),
      compute: {
        module: wavePacketShaderModule,
        entryPoint: "computeTimeValues"
      }
    });
            

When we execute the shader we provide layouts and bindings for all of the input variables from both entry points. While it might be possible to skip some of them for one entry point or another, on some platforms, the spec implies that not providing bindings and layouts for all inputs should fail.


    /**
     * Evaluate the free particle wave function at a given time. At the conclusion, the timeWaveFunctionBuffer
     * is populated with the wave function values.
     *
     * @param time The time at which to evaluate the wave function.
     * @see #getTimeWaveFunctionBuffer
     */
    getTimeWaveFunction(time)
    {
      device.queue.writeBuffer(timeBuffer, 0, new Float32Array([time]), 0, 1);
      const commandEncoder = device.createCommandEncoder();

      const passEncoder = commandEncoder.beginComputePass();
      passEncoder.setPipeline(timeEvolvedPipeline);
      passEncoder.setBindGroup(0, parametersBindGroup);
      passEncoder.setBindGroup(1, wavePacketParametersBindGroup);
      passEncoder.setBindGroup(2, waveFunctionBindGroup);
      passEncoder.setBindGroup(3, timeBindGroup);
      const workgroupCountX = Math.ceil(xResolution / 64);
      passEncoder.dispatchWorkgroups(workgroupCountX);
      passEncoder.end();

      // Submit GPU commands.
      const gpuCommands = commandEncoder.finish();
      device.queue.submit([gpuCommands]);
    }
            

A visual comparison

The next order of business is to present the exact wave function for visual comparison with the finite difference solution. For this we reuse the SchrodingerRenderer class. Yes, object reuse is a real thing. This second instance renders to a different canvas, but uses the same coloring and scaling to represent the same data.


    const exactRenderer = await SchrodingerRenderer.getInstance(schrodinger, exactCanvasID,
      reColor, imColor, psiColor, psiMax, vColor, vMax, E, yResolution,
      lineWidth);
          

At every frame we get and render the exact wave function value for the same time that we render the FDTD solution.


    const time = Math.trunc(((nframes+1)*nsteps)/3.0) * 3.0 * dt;
    exactWaveFunction.getTimeWaveFunction(time);
    exactRenderer.render(exactWaveFunction.getTimeWaveFunctionBuffer());
          

This visual comparison is a great assurance to the reader, not to mention the author, of the fidelity of the finite difference code.

A numerical comparison

We continue on to make a direct numerical comparison between our two versions of the wave function. This is a much more sensitive comparison and will show differences that are far smaller than a single pixel, but that are still important.

To do this numerical comparison, we need to retrieve the wave function values from the GPU. One of our earliest discussions of WebGPU showed how to copy data from the GPU. We put that code into a JavaScript class for ease of use.


    const computed = await dumpSchrodinger.dumpWavefunction(schrodinger.getWaveFunctionBuffer0());
    const exact = await dumpSchrodinger.dumpWavefunction(exactWaveFunction.getTimeWaveFunctionBuffer());
          

Now that we have the wave function data, we loop over it and accumulate the error measurements.


    let Er = 0;
    let Ei = 0;
    let E = 0;
    let tmp;
    for (let i = 0; i < xResolution; i++) {
      tmp = computed[2 * i] - exact[2 * i];
      Er += tmp*tmp;
      tmp = (computed[2 * i + 1] - exact[2 * i + 1]);
      Ei += tmp*tmp;
      tmp = computed[2 * i] * computed[2 * i] + computed[2 * i + 1] * computed[2 * i + 1]
            - exact[2 * i] * exact[2 * i] - exact[2 * i + 1] * exact[2 * i + 1];
      E += tmp*tmp;
    }
    Er = Math.sqrt(Er / xResolution);
    Ei = Math.sqrt(Ei / xResolution);
    E = Math.sqrt(E / xResolution);
          

Finally, we log it to the console where you will find the results for your system.


    console.log("<td>" + Er.toPrecision(5) + "</td>");
    console.log("<td>" + Ei.toPrecision(5) + "</td>");
    console.log("<td>" + E.toPrecision(5) + "</td>");
          
Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License.