Mercurial > dive4elements > river

comparison artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/tkhcalculation/TkhCalculator.java @ 8915:d9dbf0b74bc2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
8914:e3519c3e7a0a 8915:d9dbf0b74bc2
1 /** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
2 * Software engineering by
3 * Björnsen Beratende Ingenieure GmbH
4 * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
5 *
6 * This file is Free Software under the GNU AGPL (>=v3)
7 * and comes with ABSOLUTELY NO WARRANTY! Check out the
8 * documentation coming with Dive4Elements River for details.
9 */
10 package org.dive4elements.river.artifacts.sinfo.tkhcalculation;
11
12 import org.apache.commons.lang.math.DoubleRange;
13 import org.apache.commons.math.ArgumentOutsideDomainException;
14 import org.apache.commons.math.FunctionEvaluationException;
15 import org.dive4elements.artifacts.CallContext;
16 import org.dive4elements.river.artifacts.model.Calculation;
17 import org.dive4elements.river.artifacts.resources.Resources;
18 import org.dive4elements.river.artifacts.sinfo.tkhstate.BedHeightsFinder;
19 import org.dive4elements.river.model.River;
20
21 /**
22 * @author Gernot Belger
23 */
24 public final class TkhCalculator {
25
26 private static final int VALID_BED_MEASUREMENT_YEARS = 20;
27
28 private final Calculation problems;
29
30 private final String problemLabel;
31
32 private final CallContext context;
33
34 private final BedQualityD50KmValueFinder bedMeasurementsFinder;
35
36 private final SoilKindKmValueFinder soilKindFinder;
37
38 private final BedHeightsFinder bedHeightsProvider;
39
40 private final DischargeValuesFinder dischargeProvider;
41
42 private final FlowVelocityModelKmValueFinder flowVelocitiesFinder;
43
44 public static TkhCalculator buildTkhCalculator(final boolean useTkh, final CallContext context, final Calculation problems, final String label,
45 final River river, final DoubleRange calcRange, final DischargeValuesFinder dischargeProvider, final BedHeightsFinder bedHeightsProvider) {
46
47 if (!useTkh)
48 return null;
49
50 if (!dischargeProvider.isValid()) {
51 final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label);
52 problems.addProblem(message);
53 return null;
54 }
55
56 final Integer soundingYear = bedHeightsProvider.getInfo().getYear();
57 final BedQualityD50KmValueFinder bedMeasurementsFinder = BedQualityD50KmValueFinder.loadBedMeasurements(river, calcRange, soundingYear,
58 VALID_BED_MEASUREMENT_YEARS);
59
60 if (bedMeasurementsFinder == null) {
61 final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label);
62 problems.addProblem(message);
63 return null;
64 }
65
66 // FIXME: wie wird ggf. interpoliert? prüfung ob werte vorhanden?
67 final SoilKindKmValueFinder soilKindFinder = SoilKindKmValueFinder.loadValues(river, calcRange);
68 if (soilKindFinder == null) {
69 final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label);
70 problems.addProblem(message);
71 return null;
72 }
73
74 final DoubleRange qRange = dischargeProvider.getRange();
75 final FlowVelocityModelKmValueFinder flowVelocitiesFinder = FlowVelocityModelKmValueFinder.loadValues(river, calcRange, qRange);
76 if (flowVelocitiesFinder == null) {
77 final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingVelocity", null, label);
78 problems.addProblem(message);
79 return null;
80 }
81
82 return new TkhCalculator(problems, label, context, bedMeasurementsFinder, dischargeProvider, bedHeightsProvider, soilKindFinder, flowVelocitiesFinder);
83 }
84
85 private TkhCalculator(final Calculation problems, final String problemLabel, final CallContext context,
86 final BedQualityD50KmValueFinder bedMeasurementsFinder, final DischargeValuesFinder dischargeProvider, final BedHeightsFinder bedHeightsProvider,
87 final SoilKindKmValueFinder soilKindFinder,
88 final FlowVelocityModelKmValueFinder flowVelocitiesFinder) {
89 this.problems = problems;
90 this.problemLabel = problemLabel;
91 this.context = context;
92 this.bedMeasurementsFinder = bedMeasurementsFinder;
93 this.dischargeProvider = dischargeProvider;
94 this.bedHeightsProvider = bedHeightsProvider;
95 this.soilKindFinder = soilKindFinder;
96 this.flowVelocitiesFinder = flowVelocitiesFinder;
97 }
98
99 private double getDischarge(final double km) {
100
101 try {
102 return this.dischargeProvider.getDischarge(km);
103 }
104 catch (final FunctionEvaluationException e) {
105 // TODO: exceptions nicht komplett schlucken? evtl. mit log.debug(e) ausgeben
106 return Double.NaN;
107 }
108 }
109
110 private SoilKind getSoilKind(final double km) {
111
112 try {
113 return this.soilKindFinder.findSoilKind(km);
114 }
115 catch (final ArgumentOutsideDomainException e) {
116 // FIXME: cumulate problems to one message?
117 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, this.problemLabel);
118 this.problems.addProblem(km, message);
119 return null;
120 }
121 }
122
123 private double getBedMeasurement(final double km) {
124
125 try {
126 return this.bedMeasurementsFinder.findD50(km);
127 }
128 catch (final Exception e) {
129 // FIXME: cumulate problems to one message?
130 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingD50", null, this.problemLabel);
131 this.problems.addProblem(km, message);
132
133 return Double.NaN;
134 }
135 }
136
137 public Tkh getTkh(final double km, final double wst) {
138
139 final SoilKind kind = getSoilKind(km);
140
141 final double meanBedHeight = this.bedHeightsProvider.getMeanBedHeight(km);
142
143 final double discharge = getDischarge(km);
144 if (Double.isNaN(discharge)) {
145
146 // final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null,
147 // this.problemLabel);
148 // this.problems.addProblem(km, message);
149
150 // TODO: nochmal gemeinsam überlegen welche probleme wir loggen, an dieser stelle müsste man ggf. die station
151 // mitausgeben
152
153 return new Tkh(km, wst, meanBedHeight, Double.NaN, kind, Double.NaN, Double.NaN, Double.NaN);
154 }
155
156 final double d50 = getBedMeasurement(km);
157 if (Double.isNaN(d50))
158 return new Tkh(km, wst, meanBedHeight, discharge, kind, Double.NaN, Double.NaN, Double.NaN);
159
160 if (!this.flowVelocitiesFinder.findKmQValues(km, discharge)) {
161 // TODO: ggf. station in Fehlermeldung?
162 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, this.problemLabel);
163 this.problems.addProblem(km, message);
164 // FIXME: cumulate problems to one message?
165 }
166
167 final double tkh = calculateTkh(wst - meanBedHeight, this.flowVelocitiesFinder.getFindVmainFound(), d50, this.flowVelocitiesFinder.getFindTauFound());
168 // FIXME: noch mal prüfen, im alten code wurde hier immer auf 0 gesetzt
169 if (Double.isNaN(tkh) || (tkh < 0)) {
170 // TODO: ggf. station in Fehlermeldung?
171
172 // FIXME: Fehlermeldung nicht korrekt, passiert mit Wasserspiegel 'MHQ' und 'QP-1993': alle Daten (auch Abfluss)
173 // vorhanden, aber tkh negativ...
174 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, this.problemLabel);
175 this.problems.addProblem(km, message);
176
177 return new Tkh(km, wst, meanBedHeight, discharge, kind, Double.NaN, Double.NaN, Double.NaN);
178 }
179
180 /*
181 * log.debug(String.format("calculateTkh km %.3f q %.0f w %.2f mbh %.2f vm %.1f tau %.1f d50(mm) %.1f tkh(cm) %.1f",
182 * km, discharge, wst, meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), flowVelocitiesFinder.getFindTauFound(),
183 * d50*1000, tkh));
184 */
185
186 double tkhUp;
187 double tkhDown;
188 switch (kind) {
189 case starr:
190 tkhUp = tkh;
191 tkhDown = 0;
192 break;
193
194 case mobil:
195 default:
196 tkhUp = tkh / 2;
197 tkhDown = -tkh / 2;
198 break;
199 }
200
201 return new Tkh(km, wst, meanBedHeight, discharge, kind, tkh, tkhUp, tkhDown);
202 }
203
204 /**
205 * Calculates a transport body height
206 *
207 * @param h
208 * flow depth in m
209 * @param vm
210 * flow velocity in m
211 * @param d50
212 * grain diameter D50 in m (!)
213 * @param tau
214 * shear stress in N/m^2
215 * @return transport body height in cm (!)
216 */
217 private double calculateTkh(final double h, final double vm, final double d50, final double tau) {
218 final double PHYS_G = 9.81;
219 final double PHYS_SPECGRAV_S = 2.6;
220 final double PHYS_VELOCCOEFF_N = 6;
221 final double PHYS_FORMCOEFF_ALPHA = 0.7;
222 final double PHYS_VISCOSITY_NUE = 1.3e-6;
223 final double PHYS_GRAIN_DENSITY_RHOS = 2603;
224 final double PHYS_WATER_DENSITY_RHO = 999.97;
225
226 final double froude = vm / Math.sqrt(PHYS_G * h);
227 final double partReynolds = Math.sqrt((PHYS_SPECGRAV_S - 1) * PHYS_G * d50) / PHYS_VISCOSITY_NUE * d50;
228 final double critShields = 0.22 * Math.pow(partReynolds, -0.6) + 0.06 * Math.pow(10, 7.7 * Math.pow(partReynolds, -0.6));
229 final double critTau = critShields * (PHYS_GRAIN_DENSITY_RHOS - PHYS_WATER_DENSITY_RHO) * PHYS_G * d50;
230 return 100 * h * (1 - Math.pow(froude, 2)) / (2 * PHYS_VELOCCOEFF_N * PHYS_FORMCOEFF_ALPHA) * (1 - critTau / tau);
231 }
232 }
http://dive4elements.wald.intevation.org