Mercurial > dive4elements > river
diff artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/tkhcalculation/TkhCalculator.java @ 8915:d9dbf0b74bc2
Refaktoring of flow depth calculation, extracting tkh part. First implementation of tkh calculation.
author | gernotbelger |
---|---|
date | Wed, 28 Feb 2018 17:27:15 +0100 |
parents | |
children | 5d5d0051723f |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/tkhcalculation/TkhCalculator.java Wed Feb 28 17:27:15 2018 +0100 @@ -0,0 +1,232 @@ +/** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde + * Software engineering by + * Björnsen Beratende Ingenieure GmbH + * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt + * + * This file is Free Software under the GNU AGPL (>=v3) + * and comes with ABSOLUTELY NO WARRANTY! Check out the + * documentation coming with Dive4Elements River for details. + */ +package org.dive4elements.river.artifacts.sinfo.tkhcalculation; + +import org.apache.commons.lang.math.DoubleRange; +import org.apache.commons.math.ArgumentOutsideDomainException; +import org.apache.commons.math.FunctionEvaluationException; +import org.dive4elements.artifacts.CallContext; +import org.dive4elements.river.artifacts.model.Calculation; +import org.dive4elements.river.artifacts.resources.Resources; +import org.dive4elements.river.artifacts.sinfo.tkhstate.BedHeightsFinder; +import org.dive4elements.river.model.River; + +/** + * @author Gernot Belger + */ +public final class TkhCalculator { + + private static final int VALID_BED_MEASUREMENT_YEARS = 20; + + private final Calculation problems; + + private final String problemLabel; + + private final CallContext context; + + private final BedQualityD50KmValueFinder bedMeasurementsFinder; + + private final SoilKindKmValueFinder soilKindFinder; + + private final BedHeightsFinder bedHeightsProvider; + + private final DischargeValuesFinder dischargeProvider; + + private final FlowVelocityModelKmValueFinder flowVelocitiesFinder; + + public static TkhCalculator buildTkhCalculator(final boolean useTkh, final CallContext context, final Calculation problems, final String label, + final River river, final DoubleRange calcRange, final DischargeValuesFinder dischargeProvider, final BedHeightsFinder bedHeightsProvider) { + + if (!useTkh) + return null; + + if (!dischargeProvider.isValid()) { + final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label); + problems.addProblem(message); + return null; + } + + final Integer soundingYear = bedHeightsProvider.getInfo().getYear(); + final BedQualityD50KmValueFinder bedMeasurementsFinder = BedQualityD50KmValueFinder.loadBedMeasurements(river, calcRange, soundingYear, + VALID_BED_MEASUREMENT_YEARS); + + if (bedMeasurementsFinder == null) { + final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label); + problems.addProblem(message); + return null; + } + + // FIXME: wie wird ggf. interpoliert? prüfung ob werte vorhanden? + final SoilKindKmValueFinder soilKindFinder = SoilKindKmValueFinder.loadValues(river, calcRange); + if (soilKindFinder == null) { + final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label); + problems.addProblem(message); + return null; + } + + final DoubleRange qRange = dischargeProvider.getRange(); + final FlowVelocityModelKmValueFinder flowVelocitiesFinder = FlowVelocityModelKmValueFinder.loadValues(river, calcRange, qRange); + if (flowVelocitiesFinder == null) { + final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingVelocity", null, label); + problems.addProblem(message); + return null; + } + + return new TkhCalculator(problems, label, context, bedMeasurementsFinder, dischargeProvider, bedHeightsProvider, soilKindFinder, flowVelocitiesFinder); + } + + private TkhCalculator(final Calculation problems, final String problemLabel, final CallContext context, + final BedQualityD50KmValueFinder bedMeasurementsFinder, final DischargeValuesFinder dischargeProvider, final BedHeightsFinder bedHeightsProvider, + final SoilKindKmValueFinder soilKindFinder, + final FlowVelocityModelKmValueFinder flowVelocitiesFinder) { + this.problems = problems; + this.problemLabel = problemLabel; + this.context = context; + this.bedMeasurementsFinder = bedMeasurementsFinder; + this.dischargeProvider = dischargeProvider; + this.bedHeightsProvider = bedHeightsProvider; + this.soilKindFinder = soilKindFinder; + this.flowVelocitiesFinder = flowVelocitiesFinder; + } + + private double getDischarge(final double km) { + + try { + return this.dischargeProvider.getDischarge(km); + } + catch (final FunctionEvaluationException e) { + // TODO: exceptions nicht komplett schlucken? evtl. mit log.debug(e) ausgeben + return Double.NaN; + } + } + + private SoilKind getSoilKind(final double km) { + + try { + return this.soilKindFinder.findSoilKind(km); + } + catch (final ArgumentOutsideDomainException e) { + // FIXME: cumulate problems to one message? + final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, this.problemLabel); + this.problems.addProblem(km, message); + return null; + } + } + + private double getBedMeasurement(final double km) { + + try { + return this.bedMeasurementsFinder.findD50(km); + } + catch (final Exception e) { + // FIXME: cumulate problems to one message? + final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingD50", null, this.problemLabel); + this.problems.addProblem(km, message); + + return Double.NaN; + } + } + + public Tkh getTkh(final double km, final double wst) { + + final SoilKind kind = getSoilKind(km); + + final double meanBedHeight = this.bedHeightsProvider.getMeanBedHeight(km); + + final double discharge = getDischarge(km); + if (Double.isNaN(discharge)) { + + // final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, + // this.problemLabel); + // this.problems.addProblem(km, message); + + // TODO: nochmal gemeinsam überlegen welche probleme wir loggen, an dieser stelle müsste man ggf. die station + // mitausgeben + + return new Tkh(km, wst, meanBedHeight, Double.NaN, kind, Double.NaN, Double.NaN, Double.NaN); + } + + final double d50 = getBedMeasurement(km); + if (Double.isNaN(d50)) + return new Tkh(km, wst, meanBedHeight, discharge, kind, Double.NaN, Double.NaN, Double.NaN); + + if (!this.flowVelocitiesFinder.findKmQValues(km, discharge)) { + // TODO: ggf. station in Fehlermeldung? + final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, this.problemLabel); + this.problems.addProblem(km, message); + // FIXME: cumulate problems to one message? + } + + final double tkh = calculateTkh(wst - meanBedHeight, this.flowVelocitiesFinder.getFindVmainFound(), d50, this.flowVelocitiesFinder.getFindTauFound()); + // FIXME: noch mal prüfen, im alten code wurde hier immer auf 0 gesetzt + if (Double.isNaN(tkh) || (tkh < 0)) { + // TODO: ggf. station in Fehlermeldung? + + // FIXME: Fehlermeldung nicht korrekt, passiert mit Wasserspiegel 'MHQ' und 'QP-1993': alle Daten (auch Abfluss) + // vorhanden, aber tkh negativ... + final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, this.problemLabel); + this.problems.addProblem(km, message); + + return new Tkh(km, wst, meanBedHeight, discharge, kind, Double.NaN, Double.NaN, Double.NaN); + } + + /* + * log.debug(String.format("calculateTkh km %.3f q %.0f w %.2f mbh %.2f vm %.1f tau %.1f d50(mm) %.1f tkh(cm) %.1f", + * km, discharge, wst, meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), flowVelocitiesFinder.getFindTauFound(), + * d50*1000, tkh)); + */ + + double tkhUp; + double tkhDown; + switch (kind) { + case starr: + tkhUp = tkh; + tkhDown = 0; + break; + + case mobil: + default: + tkhUp = tkh / 2; + tkhDown = -tkh / 2; + break; + } + + return new Tkh(km, wst, meanBedHeight, discharge, kind, tkh, tkhUp, tkhDown); + } + + /** + * Calculates a transport body height + * + * @param h + * flow depth in m + * @param vm + * flow velocity in m + * @param d50 + * grain diameter D50 in m (!) + * @param tau + * shear stress in N/m^2 + * @return transport body height in cm (!) + */ + private double calculateTkh(final double h, final double vm, final double d50, final double tau) { + final double PHYS_G = 9.81; + final double PHYS_SPECGRAV_S = 2.6; + final double PHYS_VELOCCOEFF_N = 6; + final double PHYS_FORMCOEFF_ALPHA = 0.7; + final double PHYS_VISCOSITY_NUE = 1.3e-6; + final double PHYS_GRAIN_DENSITY_RHOS = 2603; + final double PHYS_WATER_DENSITY_RHO = 999.97; + + final double froude = vm / Math.sqrt(PHYS_G * h); + final double partReynolds = Math.sqrt((PHYS_SPECGRAV_S - 1) * PHYS_G * d50) / PHYS_VISCOSITY_NUE * d50; + final double critShields = 0.22 * Math.pow(partReynolds, -0.6) + 0.06 * Math.pow(10, 7.7 * Math.pow(partReynolds, -0.6)); + final double critTau = critShields * (PHYS_GRAIN_DENSITY_RHOS - PHYS_WATER_DENSITY_RHO) * PHYS_G * d50; + return 100 * h * (1 - Math.pow(froude, 2)) / (2 * PHYS_VELOCCOEFF_N * PHYS_FORMCOEFF_ALPHA) * (1 - critTau / tau); + } +} \ No newline at end of file