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@9522: import java.io.IOException;
gernotbelger@9522: 
gernotbelger@8915: import org.apache.commons.lang.math.DoubleRange;
gernotbelger@8997: import org.dive4elements.river.artifacts.common.GeneralResultType;
gernotbelger@8997: import org.dive4elements.river.artifacts.common.ResultRow;
gernotbelger@8915: import org.dive4elements.river.artifacts.model.Calculation;
gernotbelger@8948: import org.dive4elements.river.artifacts.sinfo.common.SInfoResultType;
gernotbelger@8915: import org.dive4elements.river.artifacts.sinfo.tkhstate.BedHeightsFinder;
gernotbelger@9522: import org.dive4elements.river.artifacts.sinfo.tkhstate.BedQualityD50TimeRangeConfig.BedQualityParseException;
gernotbelger@9522: import org.dive4elements.river.artifacts.sinfo.tkhstate.TsvHelper.TsvReaderException;
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: 
mschaefer@9335:     /**
mschaefer@9335:      * Return type of the tkh calculation
mschaefer@9335:      *
mschaefer@9335:      */
mschaefer@9335:     public enum TkhCalculateState {
mschaefer@9335:         SUCCESS, NO_W, NO_BED_HEIGHT, NO_DISCHARGE, NO_SOILKIND, NO_D50, NO_VELOCITY, NO_TKH
mschaefer@9335:     }
mschaefer@9335: 
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@9481:     public static TkhCalculator buildTkhCalculator(final boolean useTkh, final Calculation problems, final String label, final River river,
gernotbelger@9481:             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@9522:         final BedQualityD50KmValueFinder bedMeasurementsFinder = loadBedMeasurementsFinder(problems, river, calcRange, soundingYear);
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@9481:         return new TkhCalculator(bedMeasurementsFinder, waterlevelProvider, dischargeProvider, bedHeightsProvider, soilKindFinder, flowVelocitiesFinder);
gernotbelger@8915:     }
gernotbelger@8915: 
gernotbelger@9522:     private static BedQualityD50KmValueFinder loadBedMeasurementsFinder(final Calculation problems, final River river, final DoubleRange calcRange,
gernotbelger@9522:             final int soundingYear) {
gernotbelger@9522: 
gernotbelger@9522:         try {
gernotbelger@9522:             return BedQualityD50KmValueFinder.loadBedMeasurements(problems, river, calcRange, soundingYear);
gernotbelger@9522:         }
gernotbelger@9522:         catch (final BedQualityParseException | IOException | TsvReaderException e) {
gernotbelger@9522:             e.printStackTrace();
gernotbelger@9522:             return null;
gernotbelger@9522:         }
gernotbelger@9522:     }
gernotbelger@9522: 
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: 
mschaefer@9335:     public TkhCalculateState calculateTkh(final double km, final ResultRow row) {
gernotbelger@8948: 
gernotbelger@8997:         row.putValue(GeneralResultType.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@8978:         if (Double.isNaN(wst))
mschaefer@9335:             return TkhCalculateState.NO_W;
gernotbelger@8946: 
gernotbelger@8915:         final double meanBedHeight = this.bedHeightsProvider.getMeanBedHeight(km);
gernotbelger@8948:         row.putValue(SInfoResultType.meanBedHeight, meanBedHeight);
gernotbelger@8978:         if (Double.isNaN(meanBedHeight))
mschaefer@9335:             return TkhCalculateState.NO_BED_HEIGHT;
gernotbelger@8915: 
mschaefer@9503:         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@8915: 
gernotbelger@8946:         if (!this.hasTkh())
mschaefer@9335:             return TkhCalculateState.SUCCESS;
gernotbelger@8978: 
gernotbelger@8978:         // Missing discharge or kind is only a problem if we want to calculate tkh
gernotbelger@8978:         if (Double.isNaN(discharge))
mschaefer@9335:             return TkhCalculateState.NO_DISCHARGE;
gernotbelger@8978:         if (kind == null)
mschaefer@9335:             return TkhCalculateState.NO_SOILKIND;
gernotbelger@8915: 
gernotbelger@8915:         final double d50 = getBedMeasurement(km);
mschaefer@9530:         row.putValue(SInfoResultType.d50, d50 * 1000);
gernotbelger@8915:         if (Double.isNaN(d50))
mschaefer@9335:             return TkhCalculateState.NO_D50;
gernotbelger@8915: 
gernotbelger@8964:         if (!this.flowVelocitiesFinder.findKmQValues(km, discharge))
mschaefer@9335:             return TkhCalculateState.NO_VELOCITY;
gernotbelger@8915: 
gernotbelger@8940:         final double velocity = this.flowVelocitiesFinder.getFindVmainFound();
gernotbelger@8948:         row.putValue(SInfoResultType.velocity, velocity);
gernotbelger@8978:         if (Double.isNaN(velocity))
mschaefer@9335:             return TkhCalculateState.NO_VELOCITY;
gernotbelger@8948: 
gernotbelger@8940:         final double tau = this.flowVelocitiesFinder.getFindTauFound();
gernotbelger@8948:         row.putValue(SInfoResultType.tau, tau);
gernotbelger@8978:         if (Double.isNaN(tau))
mschaefer@9335:             return TkhCalculateState.NO_VELOCITY;
gernotbelger@8940: 
mschaefer@9335:         final double tkh = calculateTkh(flowDepth, velocity, d50, tau);
gernotbelger@8948:         row.putValue(SInfoResultType.tkh, tkh);
gernotbelger@8978:         if (Double.isNaN(tkh))
mschaefer@9335:             return TkhCalculateState.NO_TKH;
gernotbelger@8948: 
gernotbelger@8915:         switch (kind) {
gernotbelger@8915:         case starr:
mschaefer@9375:             row.putValue(SInfoResultType.tkhup, Math.rint(tkh));
gernotbelger@8948:             row.putValue(SInfoResultType.tkhdown, 0.0);
gernotbelger@8915:             break;
gernotbelger@8915: 
gernotbelger@8915:         case mobil:
gernotbelger@8915:         default:
mschaefer@9375:             row.putValue(SInfoResultType.tkhup, Math.rint(tkh / 2));
mschaefer@9375:             row.putValue(SInfoResultType.tkhdown, -(Math.rint(tkh / 2)));
gernotbelger@8915:             break;
gernotbelger@8915:         }
gernotbelger@8915: 
mschaefer@9335:         final double flowDepthTkh = calculateFlowDepthTkh(tkh, kind, flowDepth);
gernotbelger@8948:         row.putValue(SInfoResultType.flowdepthtkh, flowDepthTkh);
gernotbelger@8978: 
mschaefer@9335:         return TkhCalculateState.SUCCESS;
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;
mschaefer@9121:         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: 
mschaefer@9335:     private double calculateFlowDepthTkh(final double tkhValue, final SoilKind tkhKind, final double flowDepth) {
gernotbelger@8940: 
gernotbelger@8940:         switch (tkhKind) {
gernotbelger@8940:         case starr:
mschaefer@9375:             return flowDepth - Math.rint(tkhValue) / 100.0;
gernotbelger@8940: 
gernotbelger@8940:         case mobil:
gernotbelger@8940:         default:
mschaefer@9375:             return flowDepth - Math.rint(tkhValue / 2) / 100.0;
gernotbelger@8940:         }
gernotbelger@8940:     }
gernotbelger@8915: }