MacCormack scheme for riverbed variation model

code.js

let chart, worker;
let Z0, Q, B, n, s0, porosity, a, b, reachLength, simulationTime, delX, delT, delayMs, qsInjected;
let calculatedProfiles;

function updateGraph() {
    if (worker != undefined) {
        worker.terminate();
    }
    worker = new Worker('worker.js');


    Z0 = document.getElementById("z0").valueAsNumber;
    Q = document.getElementById("discharge").valueAsNumber;
    B = document.getElementById("width").valueAsNumber;
    n = document.getElementById("manningsn").valueAsNumber;
    s0 = document.getElementById("bedslope").valueAsNumber;
    porosity = document.getElementById("porosity").valueAsNumber;
    a = document.getElementById("a").valueAsNumber;
    b = document.getElementById("b").valueAsNumber;
    reachLength = document.getElementById("reach").valueAsNumber;
    simulationTime = document.getElementById("maxT").valueAsNumber;
    delX = document.getElementById("delX").valueAsNumber;
    delT = document.getElementById("delT").valueAsNumber;
    delayMs = document.getElementById("delayMs").valueAsNumber;
    let q0 = Q / B; let h0 = Math.pow(q0 * q0 * n * n / s0, 0.3);
    qsInjected = document.getElementById("qsInjected").valueAsNumber * a * Math.pow(q0 / h0, b);


    worker.onmessage = (event) => {
        let data = event.data;
        if (data.msgType == 'progress') {
            document.title = data.msgData;
        } else {
            profileCounter = 0;
            calculatedProfiles = data.msgData;
            let maxZ = findMaxZ(calculatedProfiles);
            setTimeout(function repeat() {
                profile = calculatedProfiles[profileCounter++];
                if (profile != null) {
                    drawChart(profile, calculatedProfiles[0], maxZ);
                    setTimeout(repeat, delayMs);
                }
            }, delayMs);
        }
    };

    worker.postMessage({
        Z0, Q, B, n, s0, porosity, a, b, reachLength, simulationTime, delX, delT, delayMs, qsInjected
    });

}

function findMaxZ(calculatedProfiles) {
    let maxZ = 0;
    for (let i = 0; i < calculatedProfiles.length; i++) {
        let profile = calculatedProfiles[i];
        for (let j = 0; j < profile.x.length; j++) {
            if (profile.z[j] > maxZ) {
                maxZ = profile.z[j];
            }
        }
    }
    return maxZ;
}

function drawChart(profile, initialProfile,maxZ) {
    let ctx = document.getElementById('chart').getContext('2d');
    if (chart != undefined) {
        chart.destroy();
    }
    chart = new Chart(ctx, {
        type: 'line',
        data: {
            labels: profile.x,
            datasets: [{
                label: 'Bed level z (m)',
                data: profile.z,
                borderColor: ['rgba(25, 99, 132, 1)'],
                borderWidth: 1,
                fill: false,
                pointRadius: 1
            },
            {
                label: 'Initial bed level z (m)',
                data: initialProfile.z,
                borderColor: ['rgba(255, 99, 132, 1)'],
                borderWidth: 1,
                fill: false,
                pointRadius: 1
            }/* ,
            {
                label: 'Flow depth h (m)',
                data: profile.h,
                borderColor: ['rgba(255, 99, 132, 1)'],
                borderWidth: 1,
                fill: false,
                pointRadius: 1
            },
            {
                label: 'Specific discharge q',
                data: profile.q,
                borderColor: ['rgba(255, 99, 132, 1)'],
                borderWidth: 1,
                fill: false,
                pointRadius: 1
            } */]
        },
        options: {
            animation: {
                duration: 1
            },
            scales: {
                yAxes: [{
                    ticks: {
                        beginAtZero: false,
                        autoSkip: true,
                        min: Z0 - reachLength * s0,
                        max: maxZ
                    },
                    scaleLabel: {
                        display: true,
                        labelString: 'Meters (m)'
                    }
                }],
                xAxes: [{
                    ticks: {
                        autoSkip: true,
                        min: 0,
                        max: reachLength
                    },
                    scaleLabel: {
                        display: true,
                        labelString: 'x (m)'
                    }
                }]
            },
            title: {
                display: true,
                text: profile.t + " s"
            },
            maintainAspectRatio: false,
            responsive: false
        }
    });
}

index.html

<html>

<head>
    <script src="https://cdn.jsdelivr.net/npm/chart.js@2.8.0"></script>
    <title>
        Riverbed Variation - Unsteady
    </title>
</head>

<body>

    <div id="chartDiv">
        <canvas id="chart" width="1400" height="600" style="border:1px solid">
        </canvas>
    </div>
    <div>
        <div>
            <label for="z0">Z0 (m)</label>
            <input type="range" min=1 max=10000 value=50 class="slider" id="z0"
                oninput='document.getElementById("z0Val").innerHTML=this.value'>
            <label id="z0Val">50</label>
        </div>
        <div>
            <label for="discharge">Discharge (cumecs)</label>
            <input type="range" min=0 max=100 value=14  class="slider" id="discharge"
                oninput='document.getElementById("dischargeVal").innerHTML=this.value'>
            <label id="dischargeVal">14</label>
        </div>
        <div>
            <label for="width">Width (m)</label>
            <input type="range" min=1 max=100 value=23 class="slider" id="width"
                oninput='document.getElementById("widthVal").innerHTML=this.value'>
            <label id="widthVal">23</label>
        </div>
        <div>
            <label for="manningsn">Manning's n</label>
            <input type="range" min=0.0001 max=0.1 value=0.03 step=0.001 class="slider" id="manningsn"
                oninput='document.getElementById("manningsnVal").innerHTML=this.value'>
            <label id="manningsnVal">0.03</label>
        </div>
        <div>
            <label for="porosity">Porosity</label>
            <input type="range" min=0 max=1 value=0.4 step=0.1 class="slider" id="porosity"
                oninput='document.getElementById("porosityVal").innerHTML=this.value'>
            <label id="porosityVal">0.4</label>
        </div>
        <div>
            <label for="a">a</label>
            <input type="range" min=0.0001 max=0.1 value=0.0015 step=0.0001 class="slider" id="a"
                oninput='document.getElementById("aVal").innerHTML=this.value'>
            <label id="aVal">0.0015</label>
        </div>
        <div>
            <label for="b">b</label>
            <input type="range" min=1 max=10 value=5 class="slider" id="b"
                oninput='document.getElementById("bVal").innerHTML=this.value'>
            <label id="bVal">5</label>
        </div>

        <div>
            <label for="reach">Reachlength (m)</label>
            <input type="range" min=0 max=10000 value=300 class="slider" id="reach" step=10
                oninput='document.getElementById("reachVal").innerHTML=this.value'>
            <label id="reachVal">300</label>
        </div>
        <div>
            <label for="delX">delX (m)</label>
            <input type="range" min=1 max=100 value=30 class="slider" id="delX"
                oninput='document.getElementById("delXVal").innerHTML=this.value'>
            <label id="delXVal">30</label>
        </div>
        <div>
            <label for="maxT">Simulation time (s)</label>
            <input type="range" min=1 max=10000 value=4000 class="slider" id="maxT"
                oninput='document.getElementById("maxTVal").innerHTML=this.value'>
            <label id="maxTVal">4000</label>
        </div>
        <div>
            <label for="delT">delT (s)</label>
            <input type="range" min=1 max=1000 value=3 class="slider" id="delT"
                oninput='document.getElementById("delTVal").innerHTML=this.value'>
            <label id="delTVal">3</label>
        </div>
        <div>
            <label for="bedslope">Bed slope</label>
            <input type="range" min=0 max=0.5 value=0.00225 step=0.00001 class="slider" id="bedslope"
                oninput='document.getElementById("bedslopeVal").innerHTML=this.value'>
            <label id="bedslopeVal">0.00225</label>
        </div>
        <div>
            <label for="delayMs">Delay ms</label>
            <input type="range" min=0 max=5000 value=0 class="slider" id="delayMs"
                oninput='document.getElementById("delayMsVal").innerHTML=this.value'>
            <label id="delayMsVal">0</label>
        </div>
        <div>
            <label for="qsInjected">Injected multiple of normal sediment discharge</label>
            <input type="range" min=0 max=5 value=2 step=0.1 class="slider" id="qsInjected"
                oninput='document.getElementById("qsInjectedVal").innerHTML=this.value'>
            <label id="qsInjectedVal">2</label>
        </div>
        <div>
            <input type="button" value="Calculate" id="calculate" onClick="updateGraph()">
        </div>
    </div>

    <script src="code.js"></script>
</body>

</html>

worker.js

let reportCounter, timerId;
let Z0, Q, B, n, s0, porosity, a, b, reachLength, simulationTime, delX, delT, delayMs, qsInjected;

onmessage = function (ev) {
    compute(ev.data);
}


function compute(data) {

    ({ Z0, Q, B, n, s0, porosity, a, b, reachLength, simulationTime, delX, delT, delayMs, qsInjected } = data);//javascript destructuring FTW!

    let profile = {
        t: 0,
        x: [],
        z: [],
        h: [],
        q: []
    };
    let profiles = [];

    let q0 = Q / B;
    let h0 = Math.pow(q0 * q0 * n * n / s0, 0.3);

    /*initialize boundary profile*/
    profile.t = 0;
    for (let x = 0; x <= reachLength + (simulationTime / delT + 1) * delX; x += delX) {
        profile.x.push(x);
        profile.z.push(Z0 - x * s0);
        profile.h.push(h0);
        profile.q.push(q0);
    }
    profiles.push(copyObject(profile));    //push finished profile to the profiles array
    /*end initialization*/





    for (let cnt = 0; cnt <= simulationTime / delT; cnt++) {
        let currentProfile = copyObject(profiles[cnt]); //current profile is used to calculate future profiles
        clearProfile(profile);
        profile.t = (cnt + 1) * delT;

        let suggested_delT = delX / (Math.sqrt(9.81 * currentProfile.h[0]) + currentProfile.q[0] / currentProfile.h[0]);
        console.log("Suggested delT = " + suggested_delT + "\nUsing delT = " + delT);
        if (suggested_delT < delT) {
            console.log("Unstability warning @ !" + cnt * delT);
        }

        let hi_star, qi_star, zi_star, himinus1_star, qiminus1_star, ziminus1_star;
        for (let i = 0; i <= Math.trunc(reachLength / delX + simulationTime / delT + 1); i++) {
            //information simulationTime / delT * delX downstream of the downstream boundary is required for this scheme

            //predictor
            hi_star = find_hi_star(currentProfile, i);
            if (cnt > 0) {
                qi_star = find_qi_star(currentProfile, i);
            } else {
                qi_star = q0;
            }
            if (cnt > 0) {
                zi_star = find_zi_star(currentProfile, i, qi_star, hi_star);
            } else {
                zi_star = currentProfile.z[i] - 1 / (1 - porosity) * (qs(qi_star, hi_star, a, b) * hi_star / qi_star - 0) - delT / delX * 0
            }


            himinus1_star = find_himinus1_star(currentProfile, i);
            if (cnt > 0) {
                qiminus1_star = find_qiminus1_star(currentProfile, i);
            } else {
                qiminus1_star = q0;
            }
            if (cnt > 0) {
                ziminus1_star = find_ziminus1_star(currentProfile, i, qiminus1_star, himinus1_star);
            } else {
                if (i == 0) {
                    ziminus1_star = currentProfile.z[0];
                } else {
                    ziminus1_star = currentProfile.z[i - 1] - 1 / (1 - porosity) * (qs(qiminus1_star, himinus1_star, a, b) * himinus1_star / qiminus1_star - 0) - delT / delX * 0
                }
            }


            //corrector
            let hi_star2 = hi_star - delT / delX * (qi_star - qiminus1_star);
            let qi_star2;
            if (cnt > 0) {
                qi_star2 = qi_star - delT * 9.81 * i_e(qi_star, hi_star) * hi_star - delT / delX * (p(qi_star, hi_star) - p(qiminus1_star, himinus1_star) + (zi_star - ziminus1_star) * 9.81 * hi_star);
            } else {
                qi_star2 = q0;
            }
            let zi_star2 = zi_star - 1 / (1 - porosity) * (qs(qi_star2, hi_star2, a, b) * hi_star2 / qi_star2 - qs(qi_star, hi_star, a, b) * hi_star / qi_star) - delT / delX * (1 / (1 - porosity)) * (qs(qi_star, hi_star, a, b) - qs(qiminus1_star, himinus1_star, a, b));

            //averaged result
            let hi_kplus1 = (currentProfile.h[i] + hi_star2) / 2;
            let qi_kplus1 = (currentProfile.q[i] + qi_star2) / 2;
            let zi_kplus1 = (currentProfile.z[i] + zi_star2) / 2;

            //console.log(hi_kplus1, qi_kplus1, zi_kplus1);

            profile.x.push(i * delX);
            profile.q.push(qi_kplus1);
            profile.h.push(hi_kplus1);
            profile.z.push(zi_kplus1);
        }

        profiles.push(copyObject(profile));

        reportBack({ msgType: 'progress', msgData: Math.round(profiles[cnt].t / simulationTime * 100) });


    }

    console.log('Calculation finished!');

    //prune the profiles array
    for (let i = 0; i < profiles.length; i++) {
        let tmp = profiles[i];
        let usefulLength = Math.trunc(reachLength / delX + 1);
        tmp.x.length = usefulLength;
        tmp.z.length = usefulLength;
        tmp.h.length = usefulLength;
        tmp.q.length = usefulLength;
    }

    reportBack({ msgType: 'data', msgData: profiles });
}

function reportBack(dataParcel) {
    postMessage(dataParcel);
}

function copyObject(object) {
    return JSON.parse(JSON.stringify(object));
}

function p(q, h) {
    return q * q / h + 0.5 * 9.81 * h * h;
}

function qs(q, h, a, b) {
    return a * Math.pow(q / h, b);
}

function i_e(q, h) {//calculates energy slope
    return q * q * n * n / Math.pow(h, 10 / 3);
}

function find_hi_star(profile, i) {
    return profile.h[i] - delT / delX * (profile.q[i + 1] - profile.q[i]);
}

function find_qi_star(profile, i) {
    return profile.q[i] - 9.81 * profile.h[i] * i_e(profile.q[i], profile.h[i]) * delT - delT / delX * (p(profile.q[i + 1], profile.h[i + 1]) - p(profile.q[i], profile.h[i]) + 9.81 * profile.h[i] * (profile.z[i + 1] - profile.z[i]));
}

function find_zi_star(profile, i, qi_star, hi_star) {
    if (i > 0) {
        return profile.z[i] - (1 / (1 - porosity)) * (qs(qi_star, hi_star, a, b) * hi_star / qi_star - qs(profile.q[i], profile.h[i], a, b) * profile.h[i] / profile.q[i]) - delT / delX * (1 / (1 - porosity)) * (qs(profile.q[i + 1], profile.h[i + 1], a, b) - qs(profile.q[i], profile.h[i], a, b));
    } else if (i == 0) {
        return profile.z[i] - (1 / (1 - porosity)) * (qsInjected * hi_star / qi_star - qsInjected * profile.h[i] / profile.q[i]) - delT / delX * (1 / (1 - porosity)) * (qs(profile.q[i + 1], profile.h[i + 1], a, b) - qsInjected);
    }

}

function find_qiminus1_star(profile, i) {
    if (i == 0) {
        return profile.q[0];
    } else {
        return find_qi_star(profile, i - 1);
    }
}

function find_himinus1_star(profile, i) {
    if (i == 0) {
        return profile.h[0];
    } else {
        return find_hi_star(profile, i - 1);
    }
}

function find_ziminus1_star(profile, i, qiminus1_star, himinus1_star) {
    if (i == 0) {
        return profile.z[0] + delX * (profile.z[0] - profile.z[1]) / delX;
    } else {
        return find_zi_star(profile, i - 1, qiminus1_star, himinus1_star);
    }
}

function clearProfile(profile) {
    profile.x = [];
    profile.h = [];
    profile.q = [];
    profile.z = [];
}