Mercurial > dive4elements > river
view artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/tkhcalculation/TkhCalculator.java @ 8916:5d5d0051723f
Working on outputmodes of tkh calculation
author | gernotbelger |
---|---|
date | Wed, 28 Feb 2018 18:55:39 +0100 |
parents | d9dbf0b74bc2 |
children | 29442c03c6e3 |
line wrap: on
line source
/** 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 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); } }