scirs2-integrate 0.4.2

Numerical integration module for SciRS2 (scirs2-integrate)
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

#version 430

// Runge-Kutta 4th order - Stage 3 compute shader
// Computes k3 = h * f(t + h/2, y + k2/2)

layout(local_size_x = 256, local_size_y = 1, local_size_z = 1) in;

layout(std430, binding = 0) buffer YBuffer {
    float y[];
};

layout(std430, binding = 1) buffer K2Buffer {
    float k2[];
};

layout(std430, binding = 2) buffer K3Buffer {
    float k3[];
};

uniform float t;
uniform float h;
uniform int n;

// Example ODE function: dy/dt = f(t, y)
vec4 ode_function(float time, vec4 state) {
    return vec4(state.y, -state.x, 0.0, 0.0);
}

void main() {
    uint index = gl_GlobalInvocationID.x;
    
    if (index >= n) return;
    
    // Load current state and k2
    vec4 current_state = vec4(
        y[index * 4 + 0],
        y[index * 4 + 1], 
        y[index * 4 + 2],
        y[index * 4 + 3]
    );
    
    vec4 k2_val = vec4(
        k2[index * 4 + 0],
        k2[index * 4 + 1],
        k2[index * 4 + 2],
        k2[index * 4 + 3]
    );
    
    // Compute intermediate state: y + k2/2
    vec4 intermediate_state = current_state + 0.5 * k2_val;
    
    // Compute derivative at intermediate state
    vec4 derivative = ode_function(t + 0.5 * h, intermediate_state);
    
    // Store k3 = h * f(t + h/2, y + k2/2)
    k3[index * 4 + 0] = h * derivative.x;
    k3[index * 4 + 1] = h * derivative.y;
    k3[index * 4 + 2] = h * derivative.z;
    k3[index * 4 + 3] = h * derivative.w;
}