macrorugo.ts

import { CalculatorType } from "../compute-node";
import { Dichotomie } from "../dichotomie";
import { MacrorugoCompound } from "../index";
import { ParamCalculability, ParamFamily } from "../param/param-definition";
import { ParamValueMode } from "../param/param-value-mode";
import { SessionSettings } from "../session_settings";
import { Message, MessageCode } from "../util/message";
import { Result } from "../util/result";
import { MacrorugoParams } from "./macrorugo_params";
import { MRCInclination } from "./mrc-inclination";
import { FishPass } from "../fish_pass";

export enum MacroRugoFlowType {
    EMERGENT,
    QUASI_EMERGENT,
    SUBMERGED
}

export class MacroRugo extends FishPass {

    private static readonly g = 9.81;

    /** nu: water kinematic viscosity  */
    private static readonly nu = 1E-6;
    // Water at 20 °C has a kinematic viscosity of about 10−6 m2·s−1
    // (https://en.wikipedia.org/wiki/Viscosity#Kinematic_viscosity,_%CE%BD)

    /** Ratio between the width (perpendicular to flow) and the length (parallel to flow) of a cell (-) */
    private static readonly fracAxAy = 1;

    /** Limit between emergent and submerged flow */
    private static readonly limitSubmerg = 1.1;

    /** Flag for submerged Flow */
    private flowType: MacroRugoFlowType;

    /** Velocity at the bed (m.s-1) */
    private u0: number;

    private _cache: { [key: string]: number };

    constructor(prms: MacrorugoParams, dbg: boolean = false) {
        super(prms, dbg);
        this._cache = {};
        this._calcType = CalculatorType.MacroRugo;
        this._defaultCalculatedParam = this.prms.Q;
        this.resetDefaultCalculatedParam();
    }

    /**
     * paramètres castés au bon type
     */
    get prms(): MacrorugoParams {
        return this._prms as MacrorugoParams;
    }

    /**
     * Calcul du débit total, de la cote amont ou aval ou d'un paramètre d'une structure
     * @param sVarCalc Nom du paramètre à calculer
     * @param rInit Valeur initiale
     */
    public Calc(sVarCalc: string, rInit?: number): Result {
        // Si on ne force pas à CALCUL, les tests à la 5e décimale ne passent plus (macrorugo.spec.ts)
        const originalValueMode = this.getParameter(sVarCalc).valueMode;
        this.getParameter(sVarCalc).valueMode = ParamValueMode.CALCUL;

        const r: Result = super.Calc(sVarCalc, rInit);
        this.getParameter(sVarCalc).valueMode = originalValueMode;

        // La largeur de la rampe est-elle adéquate par rapport à la largeur de motif ax ?
        // (si le parent est une MacroRugoCompound avec radier incliné, on ignore ces problèmes)
        if (
            this.parent === undefined
            || (
                this.parent instanceof MacrorugoCompound
                && this.parent.properties.getPropValue("inclinedApron") === MRCInclination.NOT_INCLINED
            )
        ) {
            const ax: number = this.prms.PBD.v / Math.sqrt(this.prms.C.v);
            if (this.prms.B.v < ax - 0.001) { // B < ax, with a little tolerance
                const m = new Message(MessageCode.WARNING_RAMP_WIDTH_LOWER_THAN_PATTERN_WIDTH);
                m.extraVar.pattern = ax;
                r.resultElement.log.add(m);
            } else if (this.prms.B.v % (ax / 2) > 0.001 && this.prms.B.v % (ax / 2) < ax / 2 - 0.001) {
                const m = new Message(MessageCode.WARNING_RAMP_WIDTH_NOT_MULTIPLE_OF_HALF_PATTERN_WIDTH);
                m.extraVar.halfPattern = ax / 2;
                m.extraVar.lower = Math.floor(this.prms.B.v / ax * 2) * (ax / 2);
                m.extraVar.higher = Math.ceil(this.prms.B.v / ax * 2) * (ax / 2);
                r.resultElement.log.add(m);
            }
        }

        // La concentration est-elle dans les valeurs admissibles 8-20% (#284)
        if (this.parent === undefined) {
            if(this.prms.C.V < 0.08 || this.prms.C.V > 0.2) {
                r.resultElement.log.add(
                    new Message(MessageCode.WARNING_MACRORUGO_CONCENTRATION_OUT_OF_BOUNDS)
                );
            }
        }

        // Ajout des résultats complémentaires
        // Cote de fond aval
        r.resultElement.values.ZF2 = this.prms.ZF1.v - this.prms.If.v * this.prms.L.v;
        // Vitesse débitante
        let resVdeb = this.prms.Q.V / this.prms.B.v / this.prms.Y.v;
        if (isNaN(resVdeb)) {
            resVdeb = 0;
        }
        r.resultElement.values.Vdeb = resVdeb;
        if (this.flowType !== MacroRugoFlowType.SUBMERGED) {
            // Froude
            r.resultElement.values.Vg =  r.resultElement.values.Vdeb / (1 - Math.sqrt(MacroRugo.fracAxAy * this.prms.C.v));
            let resFr = r.resultElement.values.Vg / Math.sqrt(MacroRugo.g * this.prms.Y.v);
            if (isNaN(resFr)) { // if Y == 0
                resFr = 0;
            }
            r.resultElement.values.Fr = resFr;
            // Vitesse maximale
            r.resultElement.values.Vmax = this.r * r.resultElement.values.Vg * this.CalcfFr(resVdeb);
        }
        // Puissance dissipée
        r.resultElement.values.PV = 1000 * MacroRugo.g * r.resultElement.values.Vdeb * this.prms.If.v;
        // Type d'écoulement
        if (this.prms.Y.v / this.prms.PBH.v < 1) {
            r.resultElement.values.ENUM_MacroRugoFlowType = MacroRugoFlowType.EMERGENT;
        } else if (this.prms.Y.v / this.prms.PBH.v < MacroRugo.limitSubmerg) {
            r.resultElement.values.ENUM_MacroRugoFlowType = MacroRugoFlowType.QUASI_EMERGENT;
        } else {
            r.resultElement.values.ENUM_MacroRugoFlowType = MacroRugoFlowType.SUBMERGED;
        }
        // Vitesse et débit du guide technique
        /* let cQ: [number, number, number, number];
        let cV: [number, number, number];
        let hdk: number;
        if (r.resultElement.values.ENUM_MacroRugoFlowType === MacroRugoFlowType.SUBMERGED) {
            cQ = [0.955, 2.282, 0.466, -0.23];
            hdk = this.prms.PBH.v;
        } else {
            hdk = this.prms.PBD.v;
            if (Math.abs(this.prms.Cd0.v - 2) < 0.05) {
                cQ = [0.648, 1.084, 0.56, -0.456];
                cV = [3.35, 0.27, 0.53];
            } else {
                cQ = [0.815, 1.45, 0.557, -0.456];
                cV = [4.54, 0.32, 0.56];
            }
        } */
        if (this.prms.Y.v > 0 && this.prms.If.v > 0) {
            r.resultElement.values.Strickler =
                this.prms.Q.V / (Math.pow(this.prms.Y.v, 5 / 3) * this.prms.B.v * Math.pow(this.prms.If.v, 0.5));
        } else {
            r.resultElement.values.Strickler = 0;
        }
        return r;
    }

    public Equation(sVarCalc: string): Result {
        const Q = this.prms.Q.v;
        const q0 = Math.sqrt(2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v * (1 - (this.sigma * this.prms.C.v)) /
            (this.Cx * this.prms.C.v)) * this.prms.Y.v * this.prms.B.v;
        let r: Result;
        if (q0 > 0) {
            this.setFlowType();
            const dicho = new Dichotomie(this, "Q", false, this.resolveQ);
            r = dicho.Dichotomie(0, SessionSettings.precision, q0);
        } else {
            r = new Result(0);
        }
        this.prms.Q.v = Q;
        return r;
    }

    /**
     * paramétrage de la calculabilité des paramètres
     */
    protected setParametersCalculability() {
        this.prms.ZF1.calculability = ParamCalculability.FREE;
        this.prms.L.calculability = ParamCalculability.FREE;
        this.prms.Ks.calculability = ParamCalculability.FREE;
        this.prms.B.calculability = ParamCalculability.DICHO;
        this.prms.If.calculability = ParamCalculability.DICHO;
        this.prms.Q.calculability = ParamCalculability.EQUATION;
        this.prms.Y.calculability = ParamCalculability.DICHO;
        this.prms.C.calculability = ParamCalculability.DICHO;
        this.prms.PBD.calculability = ParamCalculability.FREE;
        this.prms.PBH.calculability = ParamCalculability.FREE;
        this.prms.Cd0.calculability = ParamCalculability.FREE;
    }

    protected setResultsUnits() {
        this._resultsUnits = {
            PV: "W/m³",
            Vdeb: "m/s",
            Vmax: "m/s",
            Vg: "m/s",
            ZF2: "m",
            Strickler: "SI"
        }
    }

    protected exposeResults() {
        this._resultsFamilies = {
            ZF2: ParamFamily.ELEVATIONS,
            Vdeb: ParamFamily.SPEEDS,
            Vmax: ParamFamily.SPEEDS,
            Vg: ParamFamily.SPEEDS,
            Fr: undefined,
            PV: undefined,
            Strickler: ParamFamily.STRICKLERS
        };
    }

    private setFlowType() {
        const hstar: number = this.prms.Y.v / this.prms.PBH.v;
        if (hstar > MacroRugo.limitSubmerg) {
            this.flowType = MacroRugoFlowType.SUBMERGED;
        } else if (hstar < 1) {
            this.flowType = MacroRugoFlowType.EMERGENT;
        } else {
            this.flowType = MacroRugoFlowType.QUASI_EMERGENT;
        }
    }

    /**
     * Equation from Cassan, L., Laurens, P., 2016. Design of emergent and submerged rock-ramp fish passes.
     * Knowledge & Management of Aquatic Ecosystems 45.
     * @param sVarCalc Variable à calculer
     */
    private resolveQ(): number {
        // Reset cached variables depending on Q (or not...)
        this._cache = {};

        /** adimensional water depth */
        const hstar: number = this.prms.Y.v / this.prms.PBH.v;

        switch (this.flowType) {
            case MacroRugoFlowType.SUBMERGED:
                return this.resolveQSubmerged();
            case MacroRugoFlowType.EMERGENT:
                return this.resolveQEmergent();
            case MacroRugoFlowType.QUASI_EMERGENT:
                const a = (hstar - 1) / (MacroRugo.limitSubmerg - 1);
                return (1 - a) * this.resolveQEmergent() + a * this.resolveQSubmerged();
        }
    }

    /**
     * Averaged velocity (m.s-1)
     */
    private get U0(): number {
        if (this._cache.U0 === undefined) {
            this._cache.U0 = this.prms.Q.v / this.prms.B.v / this.prms.Y.v;
        }
        return this._cache.U0;
    }

    private get CdChD(): number {
        if (this._cache.CdChD === undefined) {
            this._cache.CdChD = this.Cd * this.prms.C.v * this.prms.PBH.v / this.prms.PBD.v;
        }
        return this._cache.CdChD;
    }

    /**
     * sigma ratio between the block area in the x, y plane and D2
     */
    private get sigma(): number {
        if (this._cache.sigma === undefined) {
            if (this.prms.Cd0.v >= 2) {
                this._cache.sigma = 1;
            } else {
                this._cache.sigma = Math.PI / 4;
            }
        }
        return this._cache.sigma;
    }

    private get R(): number {
        if (this._cache.R === undefined) {
            this._cache.R = (1 - this.sigma * this.prms.C.v);
        }
        return this._cache.R;
    }

    /**
     * Bed friction coefficient Equation (3) (Cassan et al., 2016)
     * @param Y Water depth (m)
     */
    private calcCf(Y: number): number {

        if (this.prms.Ks.v < 1E-9) {
            // Between Eq (8) and (9) (Cassan et al., 2016)
            const reynolds = this.U0 * this.prms.Y.v / MacroRugo.nu;
            return 0.3164 / 4. * Math.pow(reynolds, -0.25);
        } else {
            // Equation (3) (Cassan et al., 2016)
            return 2 / Math.pow(5.1 * Math.log10(Y / this.prms.Ks.v) + 6, 2);
        }
    }

    /**
     * Interpolation of Cd0 for Cd from calibrated Cd0 (See #291)
     * Cx = 1.1 for Cd0 = 1 and Cx = 2.592 for Cd0 = 2
     */
    private get Cx(): number {
        if (this._cache.Cx === undefined) {
            this._cache.Cx = this.prms.Cd0.v * 1.4917 -0.3914;
        }
        return this._cache.Cx;
    }

    /**
     * Calculation of Cd : drag coefficient of a block under the actual flow conditions
     */
    private get Cd(): number {
        if (this._cache.Cd === undefined) {
            this._cache.Cd = this.Cx * Math.min(3, (1 + 1 / Math.pow(this.prms.Y.v / this.prms.PBD.v, 2)));
        }
        return this._cache.Cd;
    }

    /**
     * Calcul de Beta force ratio between drag and turbulent stress (Cassan et al. 2016 eq(8))
     * \Beta = (k / alpha_t) (C_d C k / D) / (1 - \sigma C)
     * @param alpha \alpha_t turbulent length scale (m) within the blocks layer
     */
    private calcBeta(alpha: number): number {
        return Math.min(100, Math.sqrt(this.prms.PBH.v * this.CdChD / alpha / this.R));
    }

    /**
     * Averaged velocity at a given vertical position (m.s-1)
     * @param alpha turbulent length scale (m) within the blocks layer
     * @param z dimensionless vertical position z / k
     */
    private calcUz(alpha: number, z: number = 1): number {
        const beta = this.calcBeta(alpha);
        // Equation (9) Cassan et al., 2016
        return this.u0 * Math.sqrt(
            beta * (this.prms.Y.v / this.prms.PBH.v - 1) * Math.sinh(beta * z) / Math.cosh(beta) + 1
        );
    }

    private get ustar(): number {
        if (this._cache.ustar === undefined) {
            this._cache.ustar = Math.sqrt(MacroRugo.g * this.prms.If.v * (this.prms.Y.v - this.prms.PBH.v));
        }
        return this._cache.ustar;
    }

    private resolveAlpha_t(alpha: number): number {
        /** s: minimum distance between blocks */
        const s = this.prms.PBD.v * (1 / Math.sqrt(this.prms.C.v) - 1);
        /** Equation(11) Cassan et al., 2016 */
        const l0 = Math.min(s, 0.15 * this.prms.PBH.v);
        // Equation(10) Cassan et al., 2016
        return alpha * this.calcUz(alpha) - l0 * this.ustar;
    }

    private resolveQSubmerged(): number {
        /** Tirant d'eau (m) */
        const h: number = this.prms.Y.v;

        /** Concentration de blocs (-) */
        const C: number = this.prms.C.v;
        /** Paramètre de bloc : Diamètre (m) */
        const D: number = this.prms.PBD.v;
        /** Paramètre de bloc : Hauteur (m) */
        const k: number = this.prms.PBH.v;
        /** Slope (m/m) */
        const S: number = this.prms.If.v;
        /** Accélération gravité (m/s²) */
        const g = MacroRugo.g;
        /** Constante von Karman */
        const kappa = 0.41;
        /** Velocity at the bed §2.3.2 Cassan et al., 2016 */
        this.u0 = Math.sqrt(k * 2 * g * S * this.R
            / (this.Cd * C * k / D + this.calcCf(k) * this.R));
        /** turbulent length scale (m) within the blocks layer (alpha_t) */
        const alpha = uniroot(this.resolveAlpha_t, this, 1E-3, 10);
        /** averaged velocity at the top of blocks (m.s-1) */
        const uk = this.calcUz(alpha);
        /** Equation (13) Cassan et al., 2016 */
        const d = k - 1 / kappa * alpha * uk / this.ustar;
        /** Equation (14) Cassan et al., 2016 */
        const z0 = (k - d) * Math.exp(- kappa * uk / this.ustar);
        /** Integral of Equation (12) Cassan et al., 2016 */
        // tslint:disable-next-line:variable-name
        let Qsup: number;
        if (z0 > 0) {
            Qsup = this.ustar / kappa * (
                (h - d) * (Math.log((h - d) / z0) - 1)
                - ((k - d) * (Math.log((k - d) / z0) - 1))
            );
        } else {
            Qsup = 0;
        }

        // calcul intégrale dans la canopée----
        // tslint:disable-next-line:variable-name
        let Qinf: number = this.u0;
        let u = 0;
        let uOld: number;
        const step = 0.01;
        const zMax = 1 + step / 2;
        for (let z = step; z < zMax; z += step) {
            uOld = u;
            u = this.calcUz(alpha, z);
            Qinf += (uOld + u);
        }
        Qinf = Qinf / 2 * step * k;

        // Calcul de u moyen
        return this.U0 - (Qinf + Qsup) / h;
    }

    private resolveQEmergent(): number {
        // tslint:disable-next-line: variable-name
        const Cd = this.Cd * this.fFr;
        /** N from Cassan 2016 eq(2) et Cassan 2014 eq(12) */
        const N = (1 * this.calcCf(this.prms.Y.v)) / (this.prms.Y.v / this.prms.PBD.v * Cd * this.prms.C.v);

        // const U0i = Math.sqrt(
        //     2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v *
        //     (1 - this.sigma * this.prms.C.v) / (Cd * this.prms.C.v * (1 + N))
        // );

        return this.U0 - uniroot(this.resolveU0Complete, this, 1E-6, 100);
    }

    private resolveU0Complete(U0: number): number {

        const fFr = this.CalcfFr(U0);

        const alpha = 1 - Math.pow(1 * this.prms.C.v, 0.5) - 0.5 * this.sigma * this.prms.C.v;

        const N = (alpha * this.calcCf(this.prms.Y.v)) /
            (this.prms.Y.v / this.prms.PBD.v * this.Cd * fFr * this.prms.C.v);

        return U0 - Math.sqrt(
            2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v *
            (1 - this.sigma * this.prms.C.v) / (this.Cd * fFr * this.prms.C.v * (1 + N))
        );
    }

    /**
     * Calcul du ratio entre la vitesse moyenne à l'aval d'un block et la vitesse maximale
     * r = 1.1 pour un plot circulaire Cd0​=1.1 et r = 1.5 pour un plot à face plane Cd0​=2.6
     * Voir #283
     */
    private get r(): number {
        if (this._cache.r === undefined) {
            this._cache.r = 0.4 * this.prms.Cd0.v + 0.7;
        }
        return this._cache.r;
    }

    /**
     * Froude correction function (Cassan et al. 2014, Eq. 19)
     */
    private get fFr(): number {
        if (this._cache.fFr === undefined) {
            this._cache.fFr = this.CalcfFr(this.U0);
        }
        return this._cache.fFr;
    }

    /**
     * Calculation of Froude correction function (Cassan et al. 2014, Eq. 19)
     */
    private CalcfFr(U0: number): number {
        // tslint:disable-next-line:variable-name
        const Fr = U0 /
            (1 - Math.sqrt(MacroRugo.fracAxAy * this.prms.C.v)) /
            Math.sqrt(MacroRugo.g * this.prms.Y.v);

        if (Fr < 1.3) {
            return Math.pow(Math.min(1 / (1 - Math.pow(Fr, 2) / 4), Math.pow(Fr, -2 / 3)), 2);
        } else {
            return Math.pow(Fr, -4 / 3);
        }
    }
}

/**
 * Searches the interval from <tt>lowerLimit</tt> to <tt>upperLimit</tt>
 * for a root (i.e., zero) of the function <tt>func</tt> with respect to
 * its first argument using Brent's method root-finding algorithm.
 *
 * Translated from zeroin.c in http://www.netlib.org/c/brent.shar.
 *
 * Copyright (c) 2012 Borgar Thorsteinsson <borgar@borgar.net>
 * MIT License, http://www.opensource.org/licenses/mit-license.php
 *
 * @param {function} func function for which the root is sought.
 * @param {number} lowerlimit the lower point of the interval to be searched.
 * @param {number} upperlimit the upper point of the interval to be searched.
 * @param {number} errorTol the desired accuracy (convergence tolerance).
 * @param {number} maxIter the maximum number of iterations.
 * @returns an estimate for the root within accuracy.
 *
 */
function uniroot<T>(func: (param: number) => number, thisArg: T, lowerLimit: number, upperLimit: number,
                    errorTol: number = 0, maxIter: number = 1000
) {
    let a = lowerLimit;
    let b = upperLimit;
    let c = a;
    let fa = func.call(thisArg, a);
    let fb = func.call(thisArg, b);
    let fc = fa;
    let tolAct;   // Actual tolerance
    let newStep;  // Step at this iteration
    let prevStep; // Distance from the last but one to the last approximation
    let p;        // Interpolation step is calculated in the form p/q; division is delayed until the last moment
    let q;

    while (maxIter-- > 0) {

        prevStep = b - a;

        if (Math.abs(fc) < Math.abs(fb)) {
            // Swap data for b to be the best approximation
            a = b, b = c, c = a;
            fa = fb, fb = fc, fc = fa;
        }

        tolAct = 1e-15 * Math.abs(b) + errorTol / 2;
        newStep = (c - b) / 2;

        if (Math.abs(newStep) <= tolAct || fb === 0) {
            return b; // Acceptable approx. is found
        }

        // Decide if the interpolation can be tried
        if (Math.abs(prevStep) >= tolAct && Math.abs(fa) > Math.abs(fb)) {
            // If prev_step was large enough and was in true direction, Interpolatiom may be tried
            let t1;
            let cb;
            let t2;
            cb = c - b;
            if (a === c) { // If we have only two distinct points linear interpolation can only be applied
                t1 = fb / fa;
                p = cb * t1;
                q = 1.0 - t1;
            } else { // Quadric inverse interpolation
                q = fa / fc, t1 = fb / fc, t2 = fb / fa;
                p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1));
                q = (q - 1) * (t1 - 1) * (t2 - 1);
            }

            if (p > 0) {
                q = -q;  // p was calculated with the opposite sign; make p positive
            } else {
                p = -p;  // and assign possible minus to q
            }

            if (p < (0.75 * cb * q - Math.abs(tolAct * q) / 2) &&
                p < Math.abs(prevStep * q / 2)) {
                // If (b + p / q) falls in [b,c] and isn't too large it is accepted
                newStep = p / q;
            }

            // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
        }

        if (Math.abs(newStep) < tolAct) { // Adjust the step to be not less than tolerance
            newStep = (newStep > 0) ? tolAct : -tolAct;
        }

        a = b, fa = fb;     // Save the previous approx.
        b += newStep, fb = func.call(thisArg, b);  // Do step to a new approxim.

        if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) {
            c = a, fc = fa; // Adjust c for it to have a sign opposite to that of b
        }
    }
    return undefined;

}