Mercurial > dive4elements > river

comparison artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flowdepth/FlowDepthCalculation.java @ 8914:e3519c3e7a0a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Workflow for SINFO-Transport bodies heights inclduing winfo calculation
author gernotbelger
date Tue, 27 Feb 2018 18:06:52 +0100
parents 37ff7f435912
children d9dbf0b74bc2
comparison
equal deleted inserted replaced
8913:924cd9943337 8914:e3519c3e7a0a
17 import java.util.List; 17 import java.util.List;
18 18
19 import org.apache.commons.lang.math.DoubleRange; 19 import org.apache.commons.lang.math.DoubleRange;
20 import org.apache.commons.math.FunctionEvaluationException; 20 import org.apache.commons.math.FunctionEvaluationException;
21 import org.apache.commons.math.analysis.UnivariateRealFunction; 21 import org.apache.commons.math.analysis.UnivariateRealFunction;
22 import org.apache.log4j.Logger;
23 import org.dive4elements.artifacts.ArtifactDatabase; 22 import org.dive4elements.artifacts.ArtifactDatabase;
24 import org.dive4elements.artifacts.CallContext; 23 import org.dive4elements.artifacts.CallContext;
25 import org.dive4elements.river.artifacts.BedHeightsArtifact; 24 import org.dive4elements.river.artifacts.BedHeightsArtifact;
26 import org.dive4elements.river.artifacts.model.Calculation; 25 import org.dive4elements.river.artifacts.model.Calculation;
27 import org.dive4elements.river.artifacts.model.CalculationResult; 26 import org.dive4elements.river.artifacts.model.CalculationResult;
45 import org.dive4elements.river.utils.GaugeIndex; 44 import org.dive4elements.river.utils.GaugeIndex;
46 import org.dive4elements.river.utils.RiverUtils; 45 import org.dive4elements.river.utils.RiverUtils;
47 46
48 class FlowDepthCalculation { 47 class FlowDepthCalculation {
49 48
50 private static Logger log = Logger.getLogger(FlowDepthCalculation.class); 49 // private static Logger log = Logger.getLogger(FlowDepthCalculation.class);
51 50
52 private static final int VALID_BED_MEASUREMENT_YEARS = 20; 51 private static final int VALID_BED_MEASUREMENT_YEARS = 20;
53 52
54 private static final String CSV_NOT_IN_GAUGE_RANGE = "export.waterlevel.csv.not.in.gauge.range"; 53 private static final String CSV_NOT_IN_GAUGE_RANGE = "export.waterlevel.csv.not.in.gauge.range";
55 54
170 UnivariateRealFunction qInterpolator = null; 169 UnivariateRealFunction qInterpolator = null;
171 DoubleRange qRange = null; 170 DoubleRange qRange = null;
172 if (doCalcTkh) { 171 if (doCalcTkh) {
173 qInterpolator = DoubleUtil.getLinearInterpolator(((QKms) wstKms).allKms(), ((QKms) wstKms).allQs()); 172 qInterpolator = DoubleUtil.getLinearInterpolator(((QKms) wstKms).allKms(), ((QKms) wstKms).allQs());
174 if (qInterpolator != null) 173 if (qInterpolator != null)
175 qRange = new DoubleRange( ((QKms) wstKms).allQs().min(), ((QKms) wstKms).allQs().max()); 174 qRange = new DoubleRange(((QKms) wstKms).allQs().min(), ((QKms) wstKms).allQs().max());
176 else { 175 else {
177 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label); 176 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label);
178 problems.addProblem(message); 177 problems.addProblem(message);
179 doCalcTkh = false; 178 doCalcTkh = false;
180 } 179 }
181 } 180 }
182 181
183 // FIXME: sort by station first, but in what direction? 182 // FIXME: sort by station first, but in what direction?
184 // FIXME: using river.getKmUp()? 183 // FIXME: using river.getKmUp()?
185 final List<BedHeightValue> values = bedHeight.getValues(); 184 final List<BedHeightValue> values = bedHeight.getValues();
186 185
187 final List<BedHeightValue> sortedValues = new ArrayList<>(values); 186 final List<BedHeightValue> sortedValues = new ArrayList<>(values);
196 doCalcTkh = false; 195 doCalcTkh = false;
197 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label); 196 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label);
198 problems.addProblem(message); 197 problems.addProblem(message);
199 } 198 }
200 } 199 }
201 200
202 FlowVelocityModelKmValueFinder flowVelocitiesFinder = null; 201 FlowVelocityModelKmValueFinder flowVelocitiesFinder = null;
203 if (doCalcTkh) { 202 if (doCalcTkh) {
204 flowVelocitiesFinder = new FlowVelocityModelKmValueFinder(); 203 flowVelocitiesFinder = new FlowVelocityModelKmValueFinder();
205 if (!flowVelocitiesFinder.loadValues(river, calcRange, qRange)) { 204 if (!flowVelocitiesFinder.loadValues(river, calcRange, qRange)) {
206 doCalcTkh = false; 205 doCalcTkh = false;
243 double tkh = 0; 242 double tkh = 0;
244 if (doCalcTkh) { 243 if (doCalcTkh) {
245 double d50 = Double.NaN; 244 double d50 = Double.NaN;
246 try { 245 try {
247 d50 = bedMeasurementsFinder.findD50(km); 246 d50 = bedMeasurementsFinder.findD50(km);
248 } catch (Exception e) { 247 }
248 catch (final Exception e) {
249 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingD50", null, label); 249 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingD50", null, label);
250 problems.addProblem(km, message); 250 problems.addProblem(km, message);
251 //FIXME: cumulate problems to one message? 251 // FIXME: cumulate problems to one message?
252 } 252 }
253 if (!Double.isNaN(d50)) { 253 if (!Double.isNaN(d50)) {
254 if (flowVelocitiesFinder.findKmQValues(km, discharge)) { 254 if (flowVelocitiesFinder.findKmQValues(km, discharge)) {
255 tkh = calculateTkh(wst - meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), d50, flowVelocitiesFinder.getFindTauFound()); 255 tkh = calculateTkh(wst - meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), d50, flowVelocitiesFinder.getFindTauFound());
256 if (!Double.isNaN(tkh) && (tkh < 0)) 256 if (!Double.isNaN(tkh) && (tkh < 0))
257 tkh = 0; 257 tkh = 0;
258 /* log.debug(String.format("calculateTkh km %.3f q %.0f w %.2f mbh %.2f vm %.1f tau %.1f d50(mm) %.1f tkh(cm) %.1f", 258 /*
259 km, discharge, wst, meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), flowVelocitiesFinder.getFindTauFound(), d50*1000, tkh)); */ 259 * log.debug(String.format("calculateTkh km %.3f q %.0f w %.2f mbh %.2f vm %.1f tau %.1f d50(mm) %.1f tkh(cm) %.1f",
260 } 260 * km, discharge, wst, meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), flowVelocitiesFinder.getFindTauFound(),
261 else { 261 * d50*1000, tkh));
262 */
263 } else {
262 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label); 264 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, label);
263 problems.addProblem(km, message); 265 problems.addProblem(km, message);
264 //FIXME: cumulate problems to one message? 266 // FIXME: cumulate problems to one message?
265 } 267 }
266 } 268 } else
267 else
268 tkh = Double.NaN; 269 tkh = Double.NaN;
269 } 270 }
270 271
271 // Soil kind 272 // Soil kind
272 SoilKind kind = SoilKind.mobil; 273 SoilKind kind = SoilKind.mobil;
273 if (doCalcTkh) { 274 if (doCalcTkh) {
274 try { 275 try {
275 kind = soilKindFinder.findSoilKind(km); 276 kind = soilKindFinder.findSoilKind(km);
276 } catch (Exception e) { 277 }
278 catch (final Exception e) {
277 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label); 279 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label);
278 problems.addProblem(km, message); 280 problems.addProblem(km, message);
279 //FIXME: cumulate problems to one message? 281 // FIXME: cumulate problems to one message?
280 } 282 }
281 } 283 }
282
283 /* // REMARK: bissl spielerei zum testen damit die sohlart nicht zu schnell wechselt
284 SoilKind kind;
285 final boolean changeKind = Math.random() > 0.95;
286 if (changeKind) {
287 switch (lastKind) {
288 case starr:
289 kind = SoilKind.mobil;
290 break;
291 case mobil:
292 default:
293 kind = SoilKind.starr;
294 break;
295 }
296 } else
297 kind = lastKind;
298 lastKind = kind;
299 */
300
301 //tkh = 100 + 10 * (Math.random() - 0.5);
302 284
303 final double flowDepthTkh; 285 final double flowDepthTkh;
304 final double tkhUp; 286 final double tkhUp;
305 final double tkhDown; 287 final double tkhDown;
306 switch (kind) { 288 switch (kind) {
324 // REMARK: access the gauge once only during calculation 306 // REMARK: access the gauge once only during calculation
325 final Gauge gauge = findGauge(waterlevel, refGauge, gaugeIndex, km); 307 final Gauge gauge = findGauge(waterlevel, refGauge, gaugeIndex, km);
326 308
327 final String gaugeLabel = gauge == null ? notinrange : gauge.getName(); 309 final String gaugeLabel = gauge == null ? notinrange : gauge.getName();
328 310
329 resultData.addRow(km, flowDepth, flowDepthTkh, kind, tkh, tkhUp, tkhDown, wst, discharge, wstLabel, 311 resultData.addRow(km, flowDepth, flowDepthTkh, kind, tkh, tkhUp, tkhDown, wst, discharge, wstLabel, gaugeLabel, meanBedHeight, bedHeightLabel,
330 gaugeLabel, meanBedHeight, bedHeightLabel, location); 312 location);
331 } 313 }
332 catch (final FunctionEvaluationException e) { 314 catch (final FunctionEvaluationException e) {
333 /* should only happen if out of range */ 315 /* should only happen if out of range */
334 e.printStackTrace(); 316 e.printStackTrace();
335 /* simply ignore */ 317 /* simply ignore */
336 } 318 }
337
338 } 319 }
339 320
340 return resultData; 321 return resultData;
341 } 322 }
342 323
468 // BedHeightFactory uses its own (direct) way of accessing the data, with its own implemented data classes. 449 // BedHeightFactory uses its own (direct) way of accessing the data, with its own implemented data classes.
469 // return BedHeightFactory.getHeight(bedheightType, bedheightId, from, to); 450 // return BedHeightFactory.getHeight(bedheightType, bedheightId, from, to);
470 451
471 return BedHeight.getBedHeightById(bedheightId); 452 return BedHeight.getBedHeightById(bedheightId);
472 } 453 }
473 454
474 /** 455 /**
475 * Calculates a transport body height 456 * Calculates a transport body height
476 * @param h flow depth in m 457 *
477 * @param vm flow velocity in m 458 * @param h
478 * @param d50 grain diameter D50 in m (!) 459 * flow depth in m
479 * @param tau shear stress in N/m^2 460 * @param vm
461 * flow velocity in m
462 * @param d50
463 * grain diameter D50 in m (!)
464 * @param tau
465 * shear stress in N/m^2
480 * @return transport body height in cm (!) 466 * @return transport body height in cm (!)
481 */ 467 */
482 private double calculateTkh(double h, double vm, double d50, double tau) { 468 private double calculateTkh(final double h, final double vm, final double d50, final double tau) {
483 final double PHYS_G = 9.81; 469 final double PHYS_G = 9.81;
484 final double PHYS_SPECGRAV_S = 2.6; 470 final double PHYS_SPECGRAV_S = 2.6;
485 final double PHYS_VELOCCOEFF_N = 6; 471 final double PHYS_VELOCCOEFF_N = 6;
486 final double PHYS_FORMCOEFF_ALPHA = 0.7; 472 final double PHYS_FORMCOEFF_ALPHA = 0.7;
487 final double PHYS_VISCOSITY_NUE = 1.3e-6; 473 final double PHYS_VISCOSITY_NUE = 1.3e-6;
488 final double PHYS_GRAIN_DENSITY_RHOS = 2603; 474 final double PHYS_GRAIN_DENSITY_RHOS = 2603;
489 final double PHYS_WATER_DENSITY_RHO = 999.97; 475 final double PHYS_WATER_DENSITY_RHO = 999.97;
490 476
491 final double froude = vm / Math.sqrt(PHYS_G * h); 477 final double froude = vm / Math.sqrt(PHYS_G * h);
492 final double partReynolds = Math.sqrt((PHYS_SPECGRAV_S - 1) * PHYS_G * d50) / PHYS_VISCOSITY_NUE * d50; 478 final double partReynolds = Math.sqrt((PHYS_SPECGRAV_S - 1) * PHYS_G * d50) / PHYS_VISCOSITY_NUE * d50;
493 final double critShields = 0.22 * Math.pow(partReynolds, -0.6) + 0.06 * Math.pow(10, 7.7 * Math.pow(partReynolds, -0.6)); 479 final double critShields = 0.22 * Math.pow(partReynolds, -0.6) + 0.06 * Math.pow(10, 7.7 * Math.pow(partReynolds, -0.6));
494 final double critTau = critShields * (PHYS_GRAIN_DENSITY_RHOS - PHYS_WATER_DENSITY_RHO) * PHYS_G * d50; 480 final double critTau = critShields * (PHYS_GRAIN_DENSITY_RHOS - PHYS_WATER_DENSITY_RHO) * PHYS_G * d50;
495 return 100 * h * (1 - Math.pow(froude, 2)) / (2 * PHYS_VELOCCOEFF_N * PHYS_FORMCOEFF_ALPHA) * (1 - critTau / tau); 481 return 100 * h * (1 - Math.pow(froude, 2)) / (2 * PHYS_VELOCCOEFF_N * PHYS_FORMCOEFF_ALPHA) * (1 - critTau / tau);
496 } 482 }
497 } 483 }
http://dive4elements.wald.intevation.org