//----------------------------------------------------------------------- // // Copyright (c) 2010 Deltares. All rights reserved. // // B.S.T.I.M. The // tom.the@deltares.nl // 17-11-2010 // Perform DAM calculations for a time serie of waterlevels //----------------------------------------------------------------------- using Deltares.Geotechnics; using Deltares.Mathematics; using Deltares.Standard; using Deltares.Standard.Application; using Deltares.Standard.Extensions; using Deltares.Standard.Logging; namespace Deltares.Dam.Data { using System; using System.Collections.Generic; using System.Linq; using System.IO; using Assemblers; using System.Xml.Linq; using System.Text.RegularExpressions; public class TimeSerieStabilityCalculatorException : ApplicationException { public TimeSerieStabilityCalculatorException(string message) : base(message) { } } public class TimeSerieStabilityCalculator { private readonly LogHelper logger = LogHelper.Create(); private readonly Dictionary> stabTimeSerieEntry = new Dictionary>(); private SendMessageDelegate sendMessageDelegate = null; public const string LocationMappingExtension = "_locationmapping.txt"; public string StabilityWorkingPath { get; set; } public string MStabExePath { get; set; } public bool IsMStabCalculationOff { get; set; } public string StatesInputPath { get; set; } public string StatesOutputPath { get; set; } /// /// This boolean is used when the batch is too large and too memory leaks occur, causing the calculation to fail /// Only set this to false when above is a problem, because the calculation time will increase very much /// public bool IsCalculateAllStabilityProjectsAtOnce { get; set; } public SendMessageDelegate SendMessageDelegate { get { return sendMessageDelegate; } set { sendMessageDelegate = value; } } /// /// Constructor /// public TimeSerieStabilityCalculator() { StabilityWorkingPath = ""; IsMStabCalculationOff = false; IsCalculateAllStabilityProjectsAtOnce = true; } /// /// /// private void Validate(TimeSerie timeSerieIn, Location location, MStabParameters mstabParameters) { ThrowHelper.ThrowWhenConditionIsTrue("DikeFlow static calculation not supported", () => location.PLLineCreationMethod == PLLineCreationMethod.DupuitStatic); ThrowHelper.ThrowWhenConditionIsTrue("DikeFlow dynamic restricted to 999 time serie entries", () => (location.PLLineCreationMethod == PLLineCreationMethod.DupuitDynamic) && timeSerieIn.Entries.Count > 999); ThrowHelper.ThrowWhenConditionIsTrue("DikeFlow dynamic combined Bishop Upliftvan calculation not supported", () => (location.PLLineCreationMethod == PLLineCreationMethod.DupuitDynamic) && mstabParameters.IsCombinedBishopUpliftVanCalculation); } /// /// Create SafetyFactor TimeSerie for StabilityInside /// /// /// /// /// /// /// /// /// public TimeSerie CreateStabilityInsideSafetyFactorTimeSerie(TimeSerie timeSerieIn, Dike dike, Location location, int locationCounter, string dataDirectory, MStabParameters mstabParameters, IEnumerable gaugeTimeSerieList) { return CreateStabilitySafetyFactorTimeSerie(FailureMechanismSystemType.StabilityInside, timeSerieIn, dike, location, locationCounter, dataDirectory, mstabParameters, gaugeTimeSerieList); } public TimeSerie CreateStabilityInsideSafetyFactorTimeSerie(TimeSerie timeSerieIn, Dike dike, Location location, int locationCounter, MStabParameters mstabParameters, IEnumerable gaugeTimeSerieList) { return CreateStabilitySafetyFactorTimeSerie(FailureMechanismSystemType.StabilityInside, timeSerieIn, dike, location, locationCounter, null, mstabParameters, gaugeTimeSerieList); } /// /// Create SafetyFactor TimeSerie for StabilityOutside /// public TimeSerie CreateStabilityOutsideSafetyFactorTimeSerie(TimeSerie timeSerieIn, Dike dike, Location location, int locationCounter, string dataDirectory, MStabParameters mstabParameters, IEnumerable gaugeTimeSerieList) { return CreateStabilitySafetyFactorTimeSerie(FailureMechanismSystemType.StabilityOutside, timeSerieIn, dike, location, locationCounter, dataDirectory, mstabParameters, gaugeTimeSerieList); } public TimeSerie CreateStabilityOutsideSafetyFactorTimeSerie(TimeSerie timeSerieIn, Dike dike, Location location, int locationCounter, MStabParameters mstabParameters, IEnumerable gaugeTimeSerieList) { return CreateStabilitySafetyFactorTimeSerie(FailureMechanismSystemType.StabilityOutside, timeSerieIn, dike, location, locationCounter, mstabParameters, gaugeTimeSerieList); } private TimeSerie CreateStabilitySafetyFactorTimeSerie(FailureMechanismSystemType failureMechanismType, TimeSerie timeSerieIn, Dike dike, Location location, int locationCounter, string dataDirectory, MStabParameters mstabParameters, IEnumerable gaugeTimeSerieList) { return CreateStabilitySafetyFactorTimeSerie(failureMechanismType, timeSerieIn, dike, location, locationCounter, mstabParameters, gaugeTimeSerieList); } /// /// Creates the stability safety factor time serie. /// /// Type of the failure mechanism. /// The time serie in. /// The dike. /// The location. /// The location counter. /// The mstab parameters. /// The gauge time serie list. /// private TimeSerie CreateStabilitySafetyFactorTimeSerie(FailureMechanismSystemType failureMechanismType, TimeSerie timeSerieIn, Dike dike, Location location, int locationCounter, MStabParameters mstabParameters, IEnumerable gaugeTimeSerieList) { if (failureMechanismType != FailureMechanismSystemType.StabilityInside && failureMechanismType != FailureMechanismSystemType.StabilityOutside) { throw new TimeSerieStabilityCalculatorException(String.Format("Failurmechanism '{0}' not supported in CreateStabilitySafetyFactorTimeSerie()", failureMechanismType.ToString())); } Validate(timeSerieIn, location, mstabParameters); string soilDatabasePath = location.SoildatabaseName; string serieName = TimeSerieParameters.StabilityInsideFactor.ToString(); TimeSerie serie = TimeSerie.CreateTimeSerie(timeSerieIn, serieName); int timeSerieEntryCounter = 0; SoilGeometryType soilGeometryType; string soilGeometry2DName; CalculationHelper.DetermineSoilGeometryType(location, out soilGeometryType, out soilGeometry2DName); string projectFileNameDikeFlow = CalculationHelper.GetProjectFileName(dike.Name, location, null, null, StabilityWorkingPath, null); var dupuitPLLinesTimeSerie = new DupuitPLLinesTimeSerie(); TimeSerie safetyFactorTimeSerie = null; if (location.PLLineCreationMethod == PLLineCreationMethod.DupuitDynamic) { // Create base geometry for location to calculate DikeFlow pl lines var damCalculationDikeFlow = new DamFailureMechanismeCalculationSpecification(); CalculationHelper.BuildDamCalculation(failureMechanismType, location, soilGeometry2DName, soilGeometryType, null, null, mstabParameters, MStabModelType.Bishop, damCalculationDikeFlow, soilDatabasePath, projectFileNameDikeFlow); damCalculationDikeFlow.FailureMechanismeParamatersMStab.MStabParameters.GeometryCreationOptions.PLLineAssignment = PLLineAssignment.NoPLLines; CalculationHelper.CreateMStabProjectFile(damCalculationDikeFlow.FailureMechanismeParamatersMStab, MStabExePath); // Calculate Dupuit PLLines timeserie var plLinesCreatorDupuit = new PLLinesCreatorDupuit { Geometry2DData = Geometry2DDataCreator.CreateGeometry2DDataFromGeometry2D( projectFileNameDikeFlow), SurfaceLine = location.LocalXZSurfaceLine2, WaterLevelTimeserie = timeSerieIn, SoilList = dike.SoilList, PolderLevel = location.PolderLevel, IsReverseLayerOrder = true }; dupuitPLLinesTimeSerie = plLinesCreatorDupuit.CreateAllPlLines(projectFileNameDikeFlow + ".dupuit"); safetyFactorTimeSerie = plLinesCreatorDupuit.SafetyFactorTimeserie; // Clean up if (File.Exists(projectFileNameDikeFlow)) { File.Delete(projectFileNameDikeFlow); } } //foreach (TimeSerieEntry sourceEntry in timeSerieIn.Entries) for (int timeSerieEntryIndex = 0; timeSerieEntryIndex < timeSerieIn.Entries.Count; timeSerieEntryIndex++) { TimeSerieEntry sourceEntry = timeSerieIn.Entries[timeSerieEntryIndex]; var entry = new TimeSerieEntry(); stabTimeSerieEntry.Add(entry, new List()); entry.Assign(sourceEntry); entry.Value = timeSerieIn.MissVal; // With DupuitDynamic the waterlevel will be interpolated between the available timestep waterlevel points // With DGeoStability calculation a waterlevel entry is mandatory for each timestep bool isWaterLevelAvailable = sourceEntry.Value.IsNearEqual(timeSerieIn.MissVal) && (location.PLLineCreationMethod != PLLineCreationMethod.DupuitDynamic); if (isWaterLevelAvailable) { var message = "No phreatic water level (Phreatic level = missVal): " + location.Name + " ID: " + locationCounter + " for time serie entry: " + timeSerieEntryCounter; logger.LogError(message); if (SendMessageDelegate != null) SendMessageDelegate(new LogMessage(LogMessageType.Error, dike, message)); } else { try { if (location.UsesGauges) { ProcessGauges(location, sourceEntry, gaugeTimeSerieList); } location.GaugeMissVal = timeSerieIn.MissVal; PLLines plLines = null; DupuitPLLines dupuitPLLines = null; if (location.PLLineCreationMethod != PLLineCreationMethod.DupuitDynamic) { plLines = CalculationHelper.CreateAllPLLines(sourceEntry.Value, location, soilGeometry2DName, soilGeometryType); } else { // The PL-Lines were created above with plLinesCreatorDupuit if (timeSerieEntryIndex < dupuitPLLinesTimeSerie.Entries.Count) { dupuitPLLines = dupuitPLLinesTimeSerie.Entries[timeSerieEntryIndex].DupuitPlLines; } else { dupuitPLLines = null; } } if (this.IsMStabCalculationOff) { if (location.PLLineCreationMethod == PLLineCreationMethod.DupuitDynamic && safetyFactorTimeSerie != null) { // With Dupuit (DikeFlow) calculation the stablity factors are also calculated // If no D-GeoStability calculation is made, then we use the values from DikeFlow entry.Value = safetyFactorTimeSerie.Entries[timeSerieEntryIndex].Value; } else { Random random = new Random(); entry.Value = random.NextDouble(); } } else { IList models; if (location.PLLineCreationMethod != PLLineCreationMethod.DupuitDynamic) { double? upliftFactor = CalculationHelper.GetLowestUpliftFactor(location.LocalXZSurfaceLine2, location.GetMostProbableProfile( FailureMechanismSystemType.StabilityInside), soilGeometry2DName, plLines, location); if (mstabParameters != null && !mstabParameters.IsCombinedBishopUpliftVanCalculation) models = new List { mstabParameters.Model }; else models = CalculationHelper.GetMStabModelsToCalculate(upliftFactor); } else { if (mstabParameters != null && !mstabParameters.IsCombinedBishopUpliftVanCalculation) models = new List { mstabParameters.Model }; else { // TODO: do something with ulpiftfator Dupuit models = new List { MStabModelType.Bishop }; } } foreach (MStabModelType model in models) { string projectFileName = CalculationHelper.GetProjectFileName(dike.Name, location, null, model, StabilityWorkingPath, entry.DateTime); // get all the parameters needed for the calculation var damCalculation = new DamFailureMechanismeCalculationSpecification(); CalculationHelper.BuildDamCalculation(failureMechanismType, location, soilGeometry2DName, soilGeometryType, plLines, dupuitPLLines, mstabParameters, model, damCalculation, soilDatabasePath, projectFileName); // Create the project file CalculationHelper.CreateMStabProjectFile(damCalculation.FailureMechanismeParamatersMStab, MStabExePath); stabTimeSerieEntry[entry].Add(projectFileName); //((TimeSerieEntry)entry).MStabProjectPaths.Add(projectFileName); // write line to locationmapping file var dbgMessage = "Location: " + location.Name + " ID: " + locationCounter + " for timeserie entry: " + timeSerieEntryCounter; logger.LogDebug(dbgMessage); } entry.Value = -1 * timeSerieIn.MissVal; } } catch (Exception e) { var message = "Could not create MStab project file for location: " + location.Name + " ID: " + locationCounter + " for time serie entry: " + timeSerieEntryCounter + e.Message; logger.LogError(message, e); if (SendMessageDelegate != null) SendMessageDelegate(new LogMessage(LogMessageType.Error, dike, message)); } } serie.Entries.Add(entry); timeSerieEntryCounter++; } return serie; } /// /// Process the gauge timeserie list and store in dike.gauge /// /// /// /// /// private void ProcessGauges(Location location, TimeSerieEntry sourceEntry, IEnumerable gaugeTimeSerieList) { foreach (TimeSerie gaugeTimeSerie in gaugeTimeSerieList) { Gauge gauge = null; string[] parts = gaugeTimeSerie.LocationId.Split('/'); if (parts.Count() >= 2) { string gaugeName = parts[1]; gauge = location.Gauges.FirstOrDefault(x => x.Name.Equals(gaugeName) && x.Location.Name.Equals(location.Name)); if (gauge == null) throw new TimeSerieStabilityCalculatorException(String.Format("Gauge {0}, as specified in time series, not found in dike at location {1}.", gaugeName, location.Name)); } else throw new TimeSerieStabilityCalculatorException(String.Format("One of the time series at location {0} lacks a gauge ID.", location.Name)); TimeSerieEntry gaugeEntryForTime = gaugeTimeSerie.Entries.FirstOrDefault(x => x.DateTime == sourceEntry.DateTime); if (gaugeEntryForTime == null) throw new TimeSerieStabilityCalculatorException(String.Format("Gauge water pressure time series is out of sync with water level time series at location/gauge {0}: Could not find value for this gauge at time {1}.", gaugeTimeSerie.LocationId, sourceEntry.DateTime.ToString())); gauge.Value = gaugeEntryForTime.Value; } } /// /// Create time serie with results for safety factors MStab /// /// /// /// public void CalculateSafetyFactorFromTimeSeries(string safetyFactor, Dike dike, TimeSerieCollection fewsOutputTimeSerieCollection) { if (string.IsNullOrWhiteSpace(this.StabilityWorkingPath)) throw new InvalidOperationException("Invalid path. The given stability working path is empty. Please specify the correct path"); try { if (IsCalculateAllStabilityProjectsAtOnce) { // string stabilityWorkingPath = string.IsNullOrWhiteSpace(this.StabilityWorkingPath) ? // Path.GetFullPath(".") : Path.GetFullPath(this.StabilityWorkingPath); string stabilityWorkingPath = Path.GetFullPath(this.StabilityWorkingPath); CalculationHelper.CalculateMStabProjects(stabilityWorkingPath, this.MStabExePath); } } catch (Exception e) { var message = "Error occured in MStab batch calculation; " + e.ToString(); logger.LogError(message, e); if (SendMessageDelegate != null) SendMessageDelegate(new LogMessage(LogMessageType.Error, dike, message)); } int locationCounter = 0; int timeSerieEntryCounter = 0; foreach (Location location in dike.Locations) { foreach (TimeSerie timeSerieOut in fewsOutputTimeSerieCollection.GetSeriesByLocationId(location.Name).Where(x => x.ParameterId.Equals(safetyFactor))) { foreach (TimeSerieEntry entry in timeSerieOut.Entries) { if (entry.Value != timeSerieOut.MissVal) { try { if (!IsCalculateAllStabilityProjectsAtOnce) { CalculationHelper.CalculateMStabProjects(stabTimeSerieEntry[entry], MStabExePath); } string basisFileName = ""; entry.Value = CalculationHelper.DetermineSafetyFactor(stabTimeSerieEntry[entry], ref basisFileName, MStabExePath); entry.BasisFileName = Path.GetFileNameWithoutExtension(basisFileName); if (StabilityWorkingPath.Contains(DamProject.ProjectWorkingPath)) { // This code only works when the ProjectWorkingPath is part of the StabilityWorkingPath // This is not the case in FewsDam.exe var index = DamProject.ProjectWorkingPath.Length; entry.RelativeCalculationPathName = StabilityWorkingPath.Substring(index); } } catch (Exception e) { var message = "Could not determine safety factor for location: " + location.Name + " ID: " + locationCounter + " for timeserie entry: " + timeSerieEntryCounter + e.ToString(); logger.LogError(message, e); entry.Value = 0; if (SendMessageDelegate != null) SendMessageDelegate(new LogMessage(LogMessageType.Error, dike, message)); } } } } locationCounter++; } } } }