scirs2-integrate 0.4.1

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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85

#version 430

// Runge-Kutta 4th order - Final combination compute shader
// Computes y_new = y + (k1 + 2*k2 + 2*k3 + k4) / 6

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 K1Buffer {
    float k1[];
};

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

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

layout(std430, binding = 4) buffer K4Buffer {
    float k4[];
};

layout(std430, binding = 5) buffer ResultBuffer {
    float result[];
};

uniform int n;

void main() {
    uint index = gl_GlobalInvocationID.x;
    
    if (index >= n) return;
    
    // Load current state
    vec4 current_state = vec4(
        y[index * 4 + 0],
        y[index * 4 + 1], 
        y[index * 4 + 2],
        y[index * 4 + 3]
    );
    
    // Load all k values
    vec4 k1_val = vec4(
        k1[index * 4 + 0],
        k1[index * 4 + 1],
        k1[index * 4 + 2],
        k1[index * 4 + 3]
    );
    
    vec4 k2_val = vec4(
        k2[index * 4 + 0],
        k2[index * 4 + 1],
        k2[index * 4 + 2],
        k2[index * 4 + 3]
    );
    
    vec4 k3_val = vec4(
        k3[index * 4 + 0],
        k3[index * 4 + 1],
        k3[index * 4 + 2],
        k3[index * 4 + 3]
    );
    
    vec4 k4_val = vec4(
        k4[index * 4 + 0],
        k4[index * 4 + 1],
        k4[index * 4 + 2],
        k4[index * 4 + 3]
    );
    
    // RK4 combination: y_new = y + (k1 + 2*k2 + 2*k3 + k4) / 6
    vec4 weighted_sum = k1_val + 2.0 * k2_val + 2.0 * k3_val + k4_val;
    vec4 new_state = current_state + weighted_sum / 6.0;
    
    // Store result
    result[index * 4 + 0] = new_state.x;
    result[index * 4 + 1] = new_state.y;
    result[index * 4 + 2] = new_state.z;
    result[index * 4 + 3] = new_state.w;
}