gernotbelger@8915: /** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
gernotbelger@8915:  * Software engineering by
gernotbelger@8915:  *  Björnsen Beratende Ingenieure GmbH
gernotbelger@8915:  *  Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
gernotbelger@8915:  *
gernotbelger@8915:  * This file is Free Software under the GNU AGPL (>=v3)
gernotbelger@8915:  * and comes with ABSOLUTELY NO WARRANTY! Check out the
gernotbelger@8915:  * documentation coming with Dive4Elements River for details.
gernotbelger@8915:  */
gernotbelger@8915: package org.dive4elements.river.artifacts.sinfo.tkhcalculation;
gernotbelger@8915: 
gernotbelger@8915: import org.apache.commons.lang.math.DoubleRange;
gernotbelger@8915: import org.dive4elements.river.artifacts.model.Calculation;
gernotbelger@8948: import org.dive4elements.river.artifacts.sinfo.common.SInfoResultRow;
gernotbelger@8948: import org.dive4elements.river.artifacts.sinfo.common.SInfoResultType;
gernotbelger@8915: import org.dive4elements.river.artifacts.sinfo.tkhstate.BedHeightsFinder;
gernotbelger@8915: import org.dive4elements.river.model.River;
gernotbelger@8915: 
gernotbelger@8915: /**
gernotbelger@8915:  * @author Gernot Belger
gernotbelger@8915:  */
gernotbelger@8915: public final class TkhCalculator {
gernotbelger@8915: 
gernotbelger@8915:     private static final int VALID_BED_MEASUREMENT_YEARS = 20;
gernotbelger@8915: 
gernotbelger@8915:     private final BedQualityD50KmValueFinder bedMeasurementsFinder;
gernotbelger@8915: 
gernotbelger@8915:     private final SoilKindKmValueFinder soilKindFinder;
gernotbelger@8915: 
gernotbelger@8915:     private final BedHeightsFinder bedHeightsProvider;
gernotbelger@8915: 
gernotbelger@8946:     private final WaterlevelValuesFinder waterlevelProvider;
gernotbelger@8946: 
gernotbelger@8915:     private final DischargeValuesFinder dischargeProvider;
gernotbelger@8915: 
gernotbelger@8915:     private final FlowVelocityModelKmValueFinder flowVelocitiesFinder;
gernotbelger@8915: 
gernotbelger@8964:     public static TkhCalculator buildTkhCalculator(final boolean useTkh, final Calculation problems, final String label,
gernotbelger@8946:             final River river, final DoubleRange calcRange, final WaterlevelValuesFinder waterlevelProvider, final DischargeValuesFinder dischargeProvider,
gernotbelger@8946:             final BedHeightsFinder bedHeightsProvider) {
gernotbelger@8915: 
gernotbelger@8915:         if (!useTkh)
gernotbelger@8964:             return new TkhCalculator(null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);
gernotbelger@8915: 
gernotbelger@8915:         if (!dischargeProvider.isValid()) {
gernotbelger@8964:             problems.addProblem("sinfo_calc_flow_depth.warning.missingQ", label);
gernotbelger@8964:             return new TkhCalculator(null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);
gernotbelger@8915:         }
gernotbelger@8915: 
gernotbelger@8964:         /* access bed quality data */
gernotbelger@8946:         final int soundingYear = bedHeightsProvider.getInfo().getYear();
gernotbelger@8964:         final BedQualityD50KmValueFinder bedMeasurementsFinder = BedQualityD50KmValueFinder.loadBedMeasurements(problems, river, calcRange, soundingYear,
gernotbelger@8915:                 VALID_BED_MEASUREMENT_YEARS);
gernotbelger@8964:         if (bedMeasurementsFinder == null)
gernotbelger@8964:             return new TkhCalculator(null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);
gernotbelger@8915: 
gernotbelger@8964:         /* access bed soil kind data */
gernotbelger@8964:         final SoilKindKmValueFinder soilKindFinder = SoilKindKmValueFinder.loadValues(problems, river, calcRange);
gernotbelger@8964:         if (soilKindFinder == null)
gernotbelger@8964:             return new TkhCalculator(null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);
gernotbelger@8915: 
gernotbelger@8915:         final DoubleRange qRange = dischargeProvider.getRange();
gernotbelger@8964:         final FlowVelocityModelKmValueFinder flowVelocitiesFinder = FlowVelocityModelKmValueFinder.loadValues(problems, river, calcRange, qRange);
gernotbelger@8964:         if (flowVelocitiesFinder == null)
gernotbelger@8964:             return new TkhCalculator(null, waterlevelProvider, dischargeProvider, bedHeightsProvider, null, null);
gernotbelger@8915: 
gernotbelger@8964:         return new TkhCalculator(bedMeasurementsFinder, waterlevelProvider, dischargeProvider, bedHeightsProvider, soilKindFinder,
gernotbelger@8946:                 flowVelocitiesFinder);
gernotbelger@8915:     }
gernotbelger@8915: 
gernotbelger@8964:     private TkhCalculator(final BedQualityD50KmValueFinder bedMeasurementsFinder, final WaterlevelValuesFinder waterlevelProvider,
gernotbelger@8946:             final DischargeValuesFinder dischargeProvider, final BedHeightsFinder bedHeightsProvider, final SoilKindKmValueFinder soilKindFinder,
gernotbelger@8915:             final FlowVelocityModelKmValueFinder flowVelocitiesFinder) {
gernotbelger@8915:         this.bedMeasurementsFinder = bedMeasurementsFinder;
gernotbelger@8946:         this.waterlevelProvider = waterlevelProvider;
gernotbelger@8915:         this.dischargeProvider = dischargeProvider;
gernotbelger@8915:         this.bedHeightsProvider = bedHeightsProvider;
gernotbelger@8915:         this.soilKindFinder = soilKindFinder;
gernotbelger@8915:         this.flowVelocitiesFinder = flowVelocitiesFinder;
gernotbelger@8915:     }
gernotbelger@8915: 
gernotbelger@8946:     public boolean hasTkh() {
gernotbelger@8915: 
gernotbelger@8946:         if (this.dischargeProvider == null || !this.dischargeProvider.isValid())
gernotbelger@8946:             return false;
gernotbelger@8946: 
gernotbelger@8946:         if (this.bedMeasurementsFinder == null)
gernotbelger@8946:             return false;
gernotbelger@8946: 
gernotbelger@8946:         if (this.soilKindFinder == null)
gernotbelger@8946:             return false;
gernotbelger@8946: 
gernotbelger@8946:         if (this.flowVelocitiesFinder == null)
gernotbelger@8946:             return false;
gernotbelger@8946: 
gernotbelger@8946:         return true;
gernotbelger@8915:     }
gernotbelger@8915: 
gernotbelger@8915:     private SoilKind getSoilKind(final double km) {
gernotbelger@8915: 
gernotbelger@8964:         if (this.soilKindFinder == null)
gernotbelger@8964:             return null;
gernotbelger@8954: 
gernotbelger@8964:         return this.soilKindFinder.findSoilKind(km);
gernotbelger@8915:     }
gernotbelger@8915: 
gernotbelger@8915:     private double getBedMeasurement(final double km) {
gernotbelger@8964:         return this.bedMeasurementsFinder.findD50(km);
gernotbelger@8915:     }
gernotbelger@8915: 
gernotbelger@8948:     public void calculateTkh(final double km, final SInfoResultRow row) {
gernotbelger@8948: 
gernotbelger@8948:         row.putValue(SInfoResultType.station, km);
gernotbelger@8915: 
gernotbelger@8915:         final SoilKind kind = getSoilKind(km);
gernotbelger@8948:         row.putValue(SInfoResultType.soilkind, kind);
gernotbelger@8915: 
gernotbelger@8946:         final double wst = this.waterlevelProvider.getWaterlevel(km);
gernotbelger@8948:         row.putValue(SInfoResultType.waterlevel, wst);
gernotbelger@8946: 
gernotbelger@8915:         final double meanBedHeight = this.bedHeightsProvider.getMeanBedHeight(km);
gernotbelger@8948:         row.putValue(SInfoResultType.meanBedHeight, meanBedHeight);
gernotbelger@8915: 
gernotbelger@8940:         final double flowDepth = wst - meanBedHeight;
gernotbelger@8948:         row.putValue(SInfoResultType.flowdepth, flowDepth);
gernotbelger@8940: 
gernotbelger@8946:         final double discharge = this.dischargeProvider.getDischarge(km);
gernotbelger@8948:         row.putValue(SInfoResultType.discharge, discharge);
gernotbelger@8946:         if (Double.isNaN(discharge))
gernotbelger@8948:             return;
gernotbelger@8915: 
gernotbelger@8946:         if (!this.hasTkh())
gernotbelger@8948:             return;
gernotbelger@8915: 
gernotbelger@8915:         final double d50 = getBedMeasurement(km);
gernotbelger@8915:         if (Double.isNaN(d50))
gernotbelger@8948:             return;
gernotbelger@8948:         row.putValue(SInfoResultType.d50, d50);
gernotbelger@8915: 
gernotbelger@8964:         if (!this.flowVelocitiesFinder.findKmQValues(km, discharge))
gernotbelger@8948:             return;
gernotbelger@8915: 
gernotbelger@8940:         final double velocity = this.flowVelocitiesFinder.getFindVmainFound();
gernotbelger@8948:         row.putValue(SInfoResultType.velocity, velocity);
gernotbelger@8948: 
gernotbelger@8940:         final double tau = this.flowVelocitiesFinder.getFindTauFound();
gernotbelger@8948:         row.putValue(SInfoResultType.tau, tau);
gernotbelger@8940: 
gernotbelger@8940:         final double tkh = calculateTkh(wst - meanBedHeight, velocity, d50, tau);
gernotbelger@8948:         row.putValue(SInfoResultType.tkh, tkh);
gernotbelger@8948: 
gernotbelger@8915:         switch (kind) {
gernotbelger@8915:         case starr:
gernotbelger@8948:             row.putValue(SInfoResultType.tkhup, tkh);
gernotbelger@8948:             row.putValue(SInfoResultType.tkhdown, 0.0);
gernotbelger@8915:             break;
gernotbelger@8915: 
gernotbelger@8915:         case mobil:
gernotbelger@8915:         default:
gernotbelger@8948:             row.putValue(SInfoResultType.tkhup, tkh / 2);
gernotbelger@8948:             row.putValue(SInfoResultType.tkhdown, -tkh / 2);
gernotbelger@8915:             break;
gernotbelger@8915:         }
gernotbelger@8915: 
gernotbelger@8948:         final double flowDepthTkh = calculateFlowDepthTkh(tkh, kind, wst, meanBedHeight);
gernotbelger@8948:         row.putValue(SInfoResultType.flowdepthtkh, flowDepthTkh);
gernotbelger@8915:     }
gernotbelger@8915: 
gernotbelger@8915:     /**
gernotbelger@8915:      * Calculates a transport body height
gernotbelger@8915:      *
gernotbelger@8915:      * @param h
gernotbelger@8915:      *            flow depth in m
gernotbelger@8915:      * @param vm
gernotbelger@8915:      *            flow velocity in m
gernotbelger@8915:      * @param d50
gernotbelger@8915:      *            grain diameter D50 in m (!)
gernotbelger@8915:      * @param tau
gernotbelger@8915:      *            shear stress in N/m^2
mschaefer@8920:      * @return transport body height in cm (!), never negative
gernotbelger@8915:      */
gernotbelger@8915:     private double calculateTkh(final double h, final double vm, final double d50, final double tau) {
gernotbelger@8915:         final double PHYS_G = 9.81;
gernotbelger@8915:         final double PHYS_SPECGRAV_S = 2.6;
gernotbelger@8915:         final double PHYS_VELOCCOEFF_N = 6;
gernotbelger@8915:         final double PHYS_FORMCOEFF_ALPHA = 0.7;
gernotbelger@8915:         final double PHYS_VISCOSITY_NUE = 1.3e-6;
gernotbelger@8915:         final double PHYS_GRAIN_DENSITY_RHOS = 2603;
gernotbelger@8915:         final double PHYS_WATER_DENSITY_RHO = 999.97;
gernotbelger@8915: 
gernotbelger@8915:         final double froude = vm / Math.sqrt(PHYS_G * h);
gernotbelger@8915:         final double partReynolds = Math.sqrt((PHYS_SPECGRAV_S - 1) * PHYS_G * d50) / PHYS_VISCOSITY_NUE * d50;
gernotbelger@8915:         final double critShields = 0.22 * Math.pow(partReynolds, -0.6) + 0.06 * Math.pow(10, 7.7 * Math.pow(partReynolds, -0.6));
gernotbelger@8915:         final double critTau = critShields * (PHYS_GRAIN_DENSITY_RHOS - PHYS_WATER_DENSITY_RHO) * PHYS_G * d50;
mschaefer@8920:         final double tkh = 100 * h * (1 - Math.pow(froude, 2)) / (2 * PHYS_VELOCCOEFF_N * PHYS_FORMCOEFF_ALPHA) * (1 - critTau / tau);
mschaefer@8920:         // Some regular input values may give a negative calculation result; that is unwanted
mschaefer@8920:         if (tkh < 0.0)
mschaefer@8920:             return 0.0;
gernotbelger@8938: 
gernotbelger@8938:         return tkh;
gernotbelger@8915:     }
gernotbelger@8940: 
gernotbelger@8940:     private double calculateFlowDepthTkh(final double tkhValue, final SoilKind tkhKind, final double wst, final double meanBedHeight) {
gernotbelger@8940: 
gernotbelger@8940:         if (Double.isNaN(tkhValue) || tkhKind == null)
gernotbelger@8940:             return Double.NaN;
gernotbelger@8940: 
gernotbelger@8940:         switch (tkhKind) {
gernotbelger@8940:         case starr:
gernotbelger@8940:             return wst - (meanBedHeight + tkhValue / 100);
gernotbelger@8940: 
gernotbelger@8940:         case mobil:
gernotbelger@8940:         default:
gernotbelger@8940:             return wst - (meanBedHeight + tkhValue / 200);
gernotbelger@8940:         }
gernotbelger@8940:     }
gernotbelger@8915: }