Skip to content
GitLab
Commit d418d0fb authored by Dorchies David's avatar Dorchies David
Browse files

#22 Début des tests et du débuggage

parent d35f198d
devel 296-macrorugo-test-ludovic-cassan-formulation-of-2021-03-24 325-courbe-de-remous-calculer-une-courbe-de-remous-de-passe-a-macro-rugosites 335-utiliser-des-identifiants-symboliques-plutot-que-numeriques-dans-les-fichiers 336-macrorugo-mise-a-jour-de-la-formule-f-h-a-partir-de-cassan-2023 355-cloisons-le-champ-cote-de-l-eau-amont-n-est-pas-vide-par-defaut master v2021-03-15 v4.15.0 stable nghyd_4.17.0 nghyd_4.16.3 nghyd_4.16.2 nghyd_4.16.1 nghyd_4.16.0 nghyd_4.15.1 nghyd_4.15.0 nghyd_4.14.2 nghyd_4.14.1 nghyd_4.14.0 nghyd_4.13.1 nghyd_4.13.0 nghyd_4.12.1 nghyd_4.12.0 nghyd_4.11.1 nghyd_4.11.0 nghyd_4.10.6 nghyd_4.10.5 nghyd_4.10.4 nghyd_4.10.3 nghyd_4.10.2 nghyd_4.10.1 nghyd_4.10.0 nghyd_4.9.0 nghyd_4.8.1 nghyd_4.8.0 nghyd_4.7.0 nghyd_4.6.1 nghyd_4.6.0 nghyd_4.5.0 nghyd_4.4.2 nghyd_4.4.1 nghyd_4.4.0 nghyd_4.4.0a nghyd_4.3.0 nghyd_4.2.0 nghyd_4.1.0
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 141 additions and 77 deletions
+141 -77
spec/macrorugo/macrorugo.spec.ts 0 → 100644
+ 54
− 0
View file @ d418d0fb
/**
* IMPORTANT !
* Décommenter temporairement la ligne suivante (import { } from "./mock_jasmine")
* Pour exécuter ce code dans le débugger.
* Faire de même avec le fichier test_func.ts
*/
// import { describe, expect, it, xdescribe, xit } from "../mock_jasmine";
import { ParamCalculability } from "../../src";
import { MacroRugo, MacrorugoParams } from "../../src/macrorugo/macrorugo";
import { checkResult } from "../test_func";
function macroRugoInstance(): MacroRugo {
return new MacroRugo(
new MacrorugoParams(
12.5, // ZF1
6, // L
1, // B
0.05, // If
0.5, // Q
0.6, // h
0.01, // Ks
0.05, // C
0.5, // D
0.8, // k
1.5 // Cd0
)
);
}
function testMacroRugo(varTest: string) {
describe("Calc(): ", () => {
it("V should be 1", () => {
const nub = macroRugoInstance();
nub.prms.Q.v = nub.Calc("Q").vCalc;
const res: number = nub.prms[varTest].v;
nub.prms[varTest].v = undefined;
checkResult(nub.Calc(varTest, 0), res);
});
});
}
describe("Class MacroRugo: ", () => {
const nub = macroRugoInstance();
const vCalc = nub.Calc("Q").vCalc;
for (const prm of nub.prms) {
if ([ParamCalculability.DICHO, ParamCalculability.EQUATION].includes(prm.calculability)) {
testMacroRugo(prm.symbol);
}
}
});
spec/pab/pab_dimension.spec.ts
+ 0
− 1
View file @ d418d0fb
import { PabDimension, PabDimensionParams } from "../../src/pab/pab_dimension";
import { Result } from "../../src/util/result";
import { checkResult } from "../test_func";
function pabDimensionTest(varTest: string) {
......
spec/test_func.ts
+ 1
− 1
View file @ d418d0fb
......@@ -3,7 +3,7 @@
* décommenter la ligne suivante (import { expect } from "./mock_jasmine") dans le cas d'une fonction
* qui utilise la fonction expect() de Jasmine mais qui doit être exécutée en dehors de tests Jasmine
*/
// import { describe, expect, it, xdescribe } from "./mock_jasmine";
// import { expect } from "./mock_jasmine";
import { cLog } from "../src/util/log";
import { Message, MessageCode } from "../src/util/message";
......
src/macrorugo/macrorugo.ts
+ 85
− 75
View file @ d418d0fb
......@@ -3,23 +3,25 @@ import { ParamCalculability } from "../param/param-definition";
import { Result } from "../util/result";
import { MacrorugoParams } from "./macrorugo_params";
export { MacrorugoParams };
export class MacroRugo extends Nub {
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 lenght (parallel to flow) and the width (perpendicular to flow) of a cell (-) */
private static readonly fracAyAx = 1;
/** Rugosité de fond (m) */
private ks: number;
/** Averaged velocity at the bed (m.s-1) */
private u0: number;
private _cache: { [key: string]: number }
static readonly g = 9.81;
/** nu: water kinematic viscosity */
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 lenght (parallel to flow) and the width (perpendicular to flow) of a cell (-) */
static readonly frac_ay_ax = 1;
private _cache: { [key: string]: number };
constructor(prms: MacrorugoParams, dbg: boolean = false) {
super(prms, dbg);
......@@ -33,27 +35,27 @@ export class MacroRugo extends Nub {
}
public Equation(sVarCalc: string): Result {
const Q = uniroot(this.calcQ, 0, 1E7);
const Q = uniroot(this.calcQ, this, 0, 1E7);
return new Result(Q);
}
/**
* Equation from Cassan, L., Laurens, P., 2016. Design of emergent and submerged rock-ramp fish passes. Knowledge & Management of Aquatic Ecosystems 45.
* 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
*/
public calcQ(Q: number): number {
public calcQ(this: MacroRugo, Q: number): number {
// Reset cached variables depending on Q (or not...)
this._cache = {};
/** Longueur (m) */
const L:number = this.prms.L.v;
const L: number = this.prms.L.v;
/** Tirant d'eau (m) */
const h: number = this.prms.Y.v;
/** Paramètre de bloc : Forme (1 pour rond, 2 pour carré)
* drag coefficient of a block considering a single block
infinitely high with F ≪ 1;
* infinitely high with F ≪ 1;
*/
// tslint:disable-next-line:variable-name
const Cd0: number = this.prms.Cd0.v;
/** Concentration de blocs (-) */
const C: number = this.prms.C.v;
......@@ -62,12 +64,11 @@ infinitely high with F ≪ 1;
/** Paramètre de bloc : Hauteur (m) */
const k: number = this.prms.PBH.v;
/** Pente (m/m) */
const S:number = this.prms.If.v;
const S: number = this.prms.If.v;
const g = MacroRugo.g;
const kappa = 0.41; // von Karman constant
/** Calulated average velocity */
let u: number;
if (h / k > 1.1) {
......@@ -76,7 +77,7 @@ infinitely high with F ≪ 1;
/** Velocity at the bed §2.3.2 Cassan et al., 2016 */
this.u0 = Math.sqrt(2 * g * S * D * this.R / (Cd0 * C));
/** turbulent length scale (m) within the blocks layer (alpha_t) */
const alpha = uniroot(this.calcAlpha_t,0, 100);
const alpha = uniroot(this.calcAlpha_t, this, 0, 100);
/** averaged velocity at the top of blocks (m.s-1) */
const uk = this.calcUz(alpha);
/** Equation (13) Cassan et al., 2016 */
......@@ -84,17 +85,22 @@ infinitely high with F ≪ 1;
/** Equation (14) Cassan et al., 2016 */
const z0 = (k - d) * Math.exp(- kappa * uk / this.ustar);
/** Integral of Equation (12) Cassan et al., 2016 */
const Qsup = this.ustar / kappa * ((h - d) * (Math.log((h - d) / z0) - 1) - ((k - d) * (Math.log((k - d) / z0) - 1)));
// tslint:disable-next-line:variable-name
const Qsup = this.ustar / kappa * (
(h - d) * (Math.log((h - d) / z0) - 1)
- ((k - d) * (Math.log((k - d) / z0) - 1))
);
// calcul intégrale dans la canopée----
// tslint:disable-next-line:variable-name
let Qinf: number = this.u0;
u = this.u0;
let u_old: number;
let uOld: number;
const step = 0.1;
for (let z = step; z <= 1; z += step) {
u_old = u;
uOld = u;
u = this.calcUz(alpha, z);
Qinf += (u_old + u) / 2
Qinf += (uOld + u) / 2;
}
Qinf = Qinf * step * k;
......@@ -107,14 +113,14 @@ infinitely high with F ≪ 1;
// u0 = Averaged velocity at the bed (m.s-1)
this.u0 = Q / this.prms.L.v / this.prms.Y.v;
//
u = uniroot(this.calcU0, 0, 1E7)
// Resolve equation (4) Cassan et al., 2016
u = uniroot(this.calcU0, this, 0, 1E7);
}
return this.u0 - u;
}
private get CdChD(): number {
if( this._cache.CdChD !== undefined) {
if (this._cache.CdChD !== undefined) {
return this._cache.CdChD;
}
return this.calcCd(1) * this.prms.C.v * this.prms.PBH.v / this.prms.PBD.v;
......@@ -128,7 +134,7 @@ infinitely high with F ≪ 1;
return this._cache.sigma;
}
if (this.prms.Cd0.v === 2) {
return 1
return 1;
} else {
return Math.PI / 4;
}
......@@ -146,13 +152,14 @@ infinitely high with F ≪ 1;
*/
private get Cf(): number {
// Between Eq (8) and (9) (Cassan et al., 2016)
const Re = this.u0 * this.prms.Y.v / MacroRugo.nu
// tslint:disable-next-line:variable-name
const Re = this.u0 * this.prms.Y.v / MacroRugo.nu;
if (this.prms.Ks.v < 1E-6) {
return 0.3164/4.* Math.pow(Re, -0.25);
return 0.3164 / 4. * Math.pow(Re, -0.25);
} else {
// Equation (3) (Cassan et al., 2016)
return 2 / Math.pow(5.1*Math.log10(this.prms.Y.v / this.prms.Ks.v - 1) + 6, 2);
return 2 / Math.pow(5.1 * Math.log10(this.prms.Y.v / this.prms.Ks.v - 1) + 6, 2);
}
}
......@@ -189,7 +196,7 @@ infinitely high with F ≪ 1;
* \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 {
private calcBeta(alpha: number): number {
return Math.sqrt(this.prms.PBH.v * this.CdChD / alpha / this.R);
}
......@@ -201,7 +208,9 @@ infinitely high with F ≪ 1;
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)
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 {
......@@ -211,7 +220,7 @@ infinitely high with F ≪ 1;
return Math.sqrt(MacroRugo.g * this.prms.If.v * (this.prms.Y.v - this.prms.PBH.v));
}
private calcAlpha_t(alpha: number):number {
private calcAlpha_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 */
......@@ -221,16 +230,20 @@ infinitely high with F ≪ 1;
return alpha * this.calcUz(alpha) - l0 * this.ustar;
}
private calcU0(U01: number, C: number, h: number, r: number, Cd0: number, D: number, cfmean: number, S: number, sigma: number) {
private calcU0(U01: number) {
const g = MacroRugo.g;
const alpha = 1 - (1 * C);
const alpha = 1 - (1 * this.prms.C.v);
// tslint:disable-next-line:variable-name
const Cd = this.calcCd(this.calc_fFr(U01));
/** N from Cassan 2016 eq(2) et Cassan 2014 eq(12) */
const N = (alpha * this.Cf) / (h / D * Cd * C);
const N = (alpha * this.Cf) / (this.prms.Y.v / this.prms.PBD.v * Cd * this.prms.C.v);
return U01 - Math.sqrt(2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v * (1 - this.sigma * this.prms.C.v) / (Cd * C * (1 + N)));
return U01 - 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))
);
}
/**
......@@ -239,21 +252,21 @@ infinitely high with F ≪ 1;
*/
private calc_fFr(u0: number): number {
const Fr = u0 / (1 - Math.sqrt(MacroRugo.frac_ay_ax * this.prms.C.v)) / Math.sqrt(MacroRugo.g * this.prms.Y.v);
// tslint:disable-next-line:variable-name
const Fr = u0 / (1 - Math.sqrt(MacroRugo.fracAyAx * this.prms.C.v)) / Math.sqrt(MacroRugo.g * this.prms.Y.v);
/** Interpolation linéaire entre le bloc rond (Cd0=1) et le carré (Cd0=2) */
const r = 0.4 * this.prms.Cd0.v + 0.7;
if (Fr < 1.3) {
return Math.pow(Math.min(r / (1 - Math.pow(Fr, 2) / 4), Math.pow(Fr,-2/3)), 2);
return Math.pow(Math.min(r / (1 - Math.pow(Fr, 2) / 4), Math.pow(Fr, -2 / 3)), 2);
} else {
return Math.pow(Fr,-4/3);
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
......@@ -272,27 +285,24 @@ infinitely high with F ≪ 1;
* @returns an estimate for the root within accuracy.
*
*/
function uniroot(func: Function, lowerLimit: number, upperLimit: number,
errorTol: number = 0, maxIter: number = 1000
function uniroot<T>(func: (param: number) => number, thisArg: T, lowerLimit: number, upperLimit: number,
errorTol: number = 0, maxIter: number = 1000
) {
let a = lowerLimit
, b = upperLimit
, c = a
, fa = func(a)
, fb = func(b)
, fc = fa
, s = 0
, fs = 0
, tol_act // Actual tolerance
, new_step // Step at this iteration
, prev_step // Distance from the last but one to the last approximation
, p // Interpolation step is calculated in the form p/q; division is delayed until the last moment
, q
;
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) {
prev_step = b - a;
prevStep = b - a;
if (Math.abs(fc) < Math.abs(fb)) {
// Swap data for b to be the best approximation
......@@ -300,24 +310,25 @@ function uniroot(func: Function, lowerLimit: number, upperLimit: number,
fa = fb, fb = fc, fc = fa;
}
tol_act = 1e-15 * Math.abs(b) + errorTol / 2;
new_step = (c - b) / 2;
tolAct = 1e-15 * Math.abs(b) + errorTol / 2;
newStep = (c - b) / 2;
if (Math.abs(new_step) <= tol_act || fb === 0) {
if (Math.abs(newStep) <= tolAct || fb === 0) {
return b; // Acceptable approx. is found
}
// Decide if the interpolation can be tried
if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) {
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
var t1, cb, t2;
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
} 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);
......@@ -325,26 +336,25 @@ function uniroot(func: Function, lowerLimit: number, upperLimit: number,
if (p > 0) {
q = -q; // p was calculated with the opposite sign; make p positive
}
else {
} else {
p = -p; // and assign possible minus to q
}
if (p < (0.75 * cb * q - Math.abs(tol_act * q) / 2) &&
p < Math.abs(prev_step * q / 2)) {
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
new_step = p / q;
newStep = p / q;
}
// If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
}
if (Math.abs(new_step) < tol_act) { // Adjust the step to be not less than tolerance
new_step = (new_step > 0) ? tol_act : -tol_act;
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 += new_step, fb = func(b); // Do step to a new approxim.
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
......@@ -352,4 +362,4 @@ function uniroot(func: Function, lowerLimit: number, upperLimit: number,
}
return undefined;
}
\ No newline at end of file
}
src/macrorugo/macrorugo_params.ts
+ 1
− 0
View file @ d418d0fb
......@@ -3,6 +3,7 @@ import { ParamDomainValue } from "../param/param-domain";
import { ParamsEquation } from "../param/params-equation";
export class MacrorugoParams extends ParamsEquation {
[key: string]: any; // pour pouvoir faire this['methode']();
/** Cote de fond amont (m) */
private _ZF1: ParamDefinition;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment