// Copyright (C) Stichting Deltares 2017. All rights reserved.
//
// This file is part of Ringtoets.
//
// Ringtoets is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// All names, logos, and references to "Deltares" are registered trademarks of
// Stichting Deltares and remain full property of Stichting Deltares at all times.
// All rights reserved.
using System;
using System.Collections.Generic;
using System.Linq;
using Core.Common.Base.Data;
using Ringtoets.Common.Data.Probabilistics;
using Ringtoets.MacroStabilityInwards.InputParameterCalculation;
using Ringtoets.MacroStabilityInwards.Primitives;
using Ringtoets.MacroStabilityInwards.Primitives.Exceptions;
namespace Ringtoets.MacroStabilityInwards.Data
{
///
/// Class responsible for calculating the derived macro stability inwards input.
///
public class DerivedMacroStabilityInwardsInput
{
private const double seepageLengthCoefficientOfVariation = 0.1;
private readonly MacroStabilityInwardsInput input;
///
/// Creates a new instance of .
///
/// The input to calculate the derived macro stability inwards input.
/// Thrown when is null.
public DerivedMacroStabilityInwardsInput(MacroStabilityInwardsInput input)
{
if (input == null)
{
throw new ArgumentNullException(nameof(input), @"Cannot create DerivedMacroStabilityInwardsInput without MacroStabilityInwardsInput.");
}
this.input = input;
}
///
/// Gets the piezometric head at the exit point.
/// [m]
///
public RoundedDouble PiezometricHeadExit
{
get
{
RoundedDouble dampingFactorExit = MacroStabilityInwardsSemiProbabilisticDesignValueFactory.GetDampingFactorExit(input).GetDesignValue();
RoundedDouble phreaticLevelExit = MacroStabilityInwardsSemiProbabilisticDesignValueFactory.GetPhreaticLevelExit(input).GetDesignValue();
return new RoundedDouble(2, InputParameterCalculationService.CalculatePiezometricHeadAtExit(input.AssessmentLevel,
dampingFactorExit,
phreaticLevelExit));
}
}
///
/// Gets the horizontal distance between entry and exit point.
/// [m]
///
public VariationCoefficientLogNormalDistribution SeepageLength
{
get
{
var seepageLength = new VariationCoefficientLogNormalDistribution(2);
double seepageLengthMean = input.ExitPointL - input.EntryPointL;
seepageLength.Mean = (RoundedDouble) seepageLengthMean;
seepageLength.CoefficientOfVariation = (RoundedDouble) seepageLengthCoefficientOfVariation;
return seepageLength;
}
}
///
/// Gets the total thickness of the coverage layers at the exit point.
/// [m]
///
public LogNormalDistribution ThicknessCoverageLayer
{
get
{
var thicknessCoverageLayer = new LogNormalDistribution(2)
{
Mean = RoundedDouble.NaN,
StandardDeviation = (RoundedDouble) 0.5
};
UpdateThicknessCoverageLayerMean(thicknessCoverageLayer);
return thicknessCoverageLayer;
}
}
///
/// Gets the effective thickness of the coverage layers at the exit point.
/// [m]
///
public LogNormalDistribution EffectiveThicknessCoverageLayer
{
get
{
var thicknessCoverageLayer = new LogNormalDistribution(2)
{
Mean = RoundedDouble.NaN,
StandardDeviation = (RoundedDouble) 0.5
};
UpdateEffectiveThicknessCoverageLayerMean(thicknessCoverageLayer);
return thicknessCoverageLayer;
}
}
///
/// Gets the total thickness of the aquifer layers at the exit point.
/// [m]
///
public LogNormalDistribution ThicknessAquiferLayer
{
get
{
var thicknessAquiferLayer = new LogNormalDistribution(2)
{
Mean = RoundedDouble.NaN,
StandardDeviation = (RoundedDouble) 0.5
};
UpdateThicknessAquiferLayerMean(thicknessAquiferLayer);
return thicknessAquiferLayer;
}
}
///
/// Gets the sieve size through which 70% of the grains of the top part of the aquifer pass.
/// [m]
///
public VariationCoefficientLogNormalDistribution DiameterD70
{
get
{
var distribution = new VariationCoefficientLogNormalDistribution(6)
{
Mean = RoundedDouble.NaN,
CoefficientOfVariation = RoundedDouble.NaN
};
UpdateDiameterD70Parameters(distribution);
return distribution;
}
}
///
/// Gets the Darcy-speed with which water flows through the aquifer layer.
/// [m/s]
///
public VariationCoefficientLogNormalDistribution DarcyPermeability
{
get
{
var distribution = new VariationCoefficientLogNormalDistribution(6)
{
Mean = RoundedDouble.NaN,
CoefficientOfVariation = RoundedDouble.NaN
};
UpdateDarcyPermeabilityParameters(distribution);
return distribution;
}
}
///
/// Gets the volumic weight of the saturated coverage layer.
///
public LogNormalDistribution SaturatedVolumicWeightOfCoverageLayer
{
get
{
var distribution = new LogNormalDistribution(2)
{
Mean = RoundedDouble.NaN,
StandardDeviation = RoundedDouble.NaN,
Shift = RoundedDouble.NaN
};
UpdateSaturatedVolumicWeightOfCoverageLayerParameters(distribution);
return distribution;
}
}
private void UpdateThicknessAquiferLayerMean(LogNormalDistribution thicknessAquiferLayer)
{
StochasticSoilProfile stochasticSoilProfile = input.StochasticSoilProfile;
RingtoetsMacroStabilityInwardsSurfaceLine surfaceLine = input.SurfaceLine;
RoundedDouble exitPointL = input.ExitPointL;
if (stochasticSoilProfile?.SoilProfile != null && surfaceLine != null && !double.IsNaN(exitPointL))
{
var thicknessTopAquiferLayer = new RoundedDouble(GetNumberOfDecimals(thicknessAquiferLayer),
GetThicknessTopAquiferLayer(stochasticSoilProfile.SoilProfile, surfaceLine, exitPointL));
if (thicknessTopAquiferLayer > 0)
{
thicknessAquiferLayer.Mean = thicknessTopAquiferLayer;
}
}
}
private void UpdateThicknessCoverageLayerMean(LogNormalDistribution thicknessCoverageLayerDistribution)
{
StochasticSoilProfile stochasticSoilProfile = input.StochasticSoilProfile;
RingtoetsMacroStabilityInwardsSurfaceLine surfaceLine = input.SurfaceLine;
RoundedDouble exitPointL = input.ExitPointL;
if (stochasticSoilProfile?.SoilProfile != null && surfaceLine != null && !double.IsNaN(exitPointL))
{
var weightedMean = new RoundedDouble(GetNumberOfDecimals(thicknessCoverageLayerDistribution),
GetThicknessCoverageLayers(stochasticSoilProfile.SoilProfile, surfaceLine, exitPointL));
if (weightedMean > 0)
{
thicknessCoverageLayerDistribution.Mean = weightedMean;
}
}
}
private void UpdateEffectiveThicknessCoverageLayerMean(LogNormalDistribution effectiveThicknessCoverageLayerDistribution)
{
if (input.SurfaceLine != null && input.StochasticSoilProfile?.SoilProfile != null && !double.IsNaN(input.ExitPointL))
{
var weightedMean = new RoundedDouble(GetNumberOfDecimals(effectiveThicknessCoverageLayerDistribution),
InputParameterCalculationService.CalculateEffectiveThicknessCoverageLayer(
input.WaterVolumetricWeight,
MacroStabilityInwardsSemiProbabilisticDesignValueFactory.GetPhreaticLevelExit(input).GetDesignValue(),
input.ExitPointL,
input.SurfaceLine,
input.StochasticSoilProfile.SoilProfile));
if (weightedMean > 0)
{
effectiveThicknessCoverageLayerDistribution.Mean = weightedMean;
}
}
}
private void UpdateDiameterD70Parameters(VariationCoefficientLogNormalDistribution diameterD70Distribution)
{
MacroStabilityInwardsSoilLayer topMostAquiferLayer = GetConsecutiveAquiferLayers().FirstOrDefault();
if (topMostAquiferLayer != null)
{
var diameterD70Mean = new RoundedDouble(GetNumberOfDecimals(diameterD70Distribution), topMostAquiferLayer.DiameterD70Mean);
if (diameterD70Mean > 0)
{
diameterD70Distribution.Mean = diameterD70Mean;
}
diameterD70Distribution.CoefficientOfVariation = (RoundedDouble) topMostAquiferLayer.DiameterD70CoefficientOfVariation;
}
}
private void UpdateDarcyPermeabilityParameters(VariationCoefficientLogNormalDistribution darcyPermeabilityDistribution)
{
MacroStabilityInwardsSoilLayer[] aquiferLayers = GetConsecutiveAquiferLayers();
int numberOfDecimals = GetNumberOfDecimals(darcyPermeabilityDistribution);
if (HasCorrectDarcyPermeabilityWeightDistributionParameterDefinition(
aquiferLayers,
numberOfDecimals))
{
MacroStabilityInwardsSoilLayer topMostAquiferLayer = aquiferLayers.First();
var permeabilityCoefficientOfVariation = new RoundedDouble(numberOfDecimals, topMostAquiferLayer.PermeabilityCoefficientOfVariation);
var weightedMean = new RoundedDouble(numberOfDecimals,
GetWeightedMeanForDarcyPermeabilityOfAquiferLayer(aquiferLayers,
input.StochasticSoilProfile.SoilProfile,
input.SurfaceLine.GetZAtL(input.ExitPointL)));
if (weightedMean > 0)
{
darcyPermeabilityDistribution.Mean = weightedMean;
}
darcyPermeabilityDistribution.CoefficientOfVariation = permeabilityCoefficientOfVariation;
}
}
private void UpdateSaturatedVolumicWeightOfCoverageLayerParameters(LogNormalDistribution volumicWeightDistribution)
{
MacroStabilityInwardsSoilLayer[] coverageLayers = GetConsecutiveCoverageLayers();
int numberOfDecimals = GetNumberOfDecimals(volumicWeightDistribution);
if (HasCorrectSaturatedWeightDistributionParameterDefinition(
coverageLayers,
numberOfDecimals))
{
MacroStabilityInwardsSoilLayer topMostAquitardLayer = coverageLayers.First();
volumicWeightDistribution.Shift = (RoundedDouble) topMostAquitardLayer.BelowPhreaticLevelShift;
volumicWeightDistribution.StandardDeviation = (RoundedDouble) topMostAquitardLayer.BelowPhreaticLevelDeviation;
var weightedMean = new RoundedDouble(numberOfDecimals,
GetWeightedMeanForVolumicWeightOfCoverageLayer(
coverageLayers,
input.StochasticSoilProfile.SoilProfile,
input.SurfaceLine.GetZAtL(input.ExitPointL)));
if (weightedMean > 0)
{
volumicWeightDistribution.Mean = weightedMean;
}
}
}
private static int GetNumberOfDecimals(IDistribution distribution)
{
return distribution.Mean.NumberOfDecimalPlaces;
}
private static int GetNumberOfDecimals(IVariationCoefficientDistribution distribution)
{
return distribution.Mean.NumberOfDecimalPlaces;
}
private static bool HasCorrectSaturatedWeightDistributionParameterDefinition(IList consecutiveAquitardLayers, int numberOfDecimals)
{
if (!consecutiveAquitardLayers.Any())
{
return false;
}
LogNormalDistribution[] distributions = GetLayerSaturatedVolumicWeightDistributionDefinitions(consecutiveAquitardLayers, numberOfDecimals);
if (distributions == null)
{
return false;
}
if (distributions.Length == 1)
{
return true;
}
return distributions.All(currentLayerDistribution => AreShiftAndDeviationEqual(
currentLayerDistribution,
distributions[0]));
}
private static bool HasCorrectDarcyPermeabilityWeightDistributionParameterDefinition(IList consecutiveAquitardLayers, int numberOfDecimals)
{
if (!consecutiveAquitardLayers.Any())
{
return false;
}
VariationCoefficientLogNormalDistribution[] distributions = GetLayerPermeabilityDistributionDefinitions(consecutiveAquitardLayers, numberOfDecimals);
if (distributions == null)
{
return false;
}
if (distributions.Length == 1)
{
return true;
}
return distributions.All(currentLayerDistribution => AreCoefficientEqual(
currentLayerDistribution,
distributions[0]));
}
private static LogNormalDistribution[] GetLayerSaturatedVolumicWeightDistributionDefinitions(IList consecutiveAquitardLayers, int numberOfDecimals)
{
try
{
return consecutiveAquitardLayers.Select(layer => new LogNormalDistribution(numberOfDecimals)
{
Mean = (RoundedDouble) layer.BelowPhreaticLevelMean,
StandardDeviation = (RoundedDouble) layer.BelowPhreaticLevelDeviation,
Shift = (RoundedDouble) layer.BelowPhreaticLevelShift
}).ToArray();
}
catch (ArgumentOutOfRangeException)
{
return null;
}
}
private static VariationCoefficientLogNormalDistribution[] GetLayerPermeabilityDistributionDefinitions(IList consecutiveAquitardLayers, int numberOfDecimals)
{
try
{
return consecutiveAquitardLayers.Select(layer => new VariationCoefficientLogNormalDistribution(numberOfDecimals)
{
Mean = (RoundedDouble) layer.PermeabilityMean,
CoefficientOfVariation = (RoundedDouble) layer.PermeabilityCoefficientOfVariation
}).ToArray();
}
catch (ArgumentOutOfRangeException)
{
return null;
}
}
private static bool AreShiftAndDeviationEqual(LogNormalDistribution currentLayerDistribution, LogNormalDistribution baseLayerDistribution)
{
return currentLayerDistribution.StandardDeviation == baseLayerDistribution.StandardDeviation &&
currentLayerDistribution.Shift == baseLayerDistribution.Shift;
}
private static bool AreCoefficientEqual(
VariationCoefficientLogNormalDistribution currentLayerDistribution,
VariationCoefficientLogNormalDistribution baseLayerDistribution)
{
return Math.Abs(baseLayerDistribution.CoefficientOfVariation - currentLayerDistribution.CoefficientOfVariation) < 1e-6;
}
private static double GetWeightedMeanForVolumicWeightOfCoverageLayer(MacroStabilityInwardsSoilLayer[] aquitardLayers, MacroStabilityInwardsSoilProfile profile, double surfaceLevel)
{
var totalThickness = 0.0;
var weighedTotal = 0.0;
foreach (MacroStabilityInwardsSoilLayer layer in aquitardLayers)
{
double layerThickness = profile.GetLayerThickness(layer);
double bottom = layer.Top - layerThickness;
double thicknessUnderSurface = Math.Min(layer.Top, surfaceLevel) - bottom;
totalThickness += thicknessUnderSurface;
weighedTotal += layer.BelowPhreaticLevelMean * thicknessUnderSurface;
}
return weighedTotal / totalThickness;
}
private static double GetWeightedMeanForDarcyPermeabilityOfAquiferLayer(MacroStabilityInwardsSoilLayer[] aquitardLayers, MacroStabilityInwardsSoilProfile profile, double surfaceLevel)
{
var totalThickness = 0.0;
var weighedTotal = 0.0;
foreach (MacroStabilityInwardsSoilLayer layer in aquitardLayers)
{
double layerThickness = profile.GetLayerThickness(layer);
double bottom = layer.Top - layerThickness;
double thicknessUnderSurface = Math.Min(layer.Top, surfaceLevel) - bottom;
totalThickness += thicknessUnderSurface;
weighedTotal += layer.PermeabilityMean * thicknessUnderSurface;
}
return weighedTotal / totalThickness;
}
private MacroStabilityInwardsSoilLayer[] GetConsecutiveAquiferLayers()
{
RingtoetsMacroStabilityInwardsSurfaceLine surfaceLine = input.SurfaceLine;
MacroStabilityInwardsSoilProfile soilProfile = input.StochasticSoilProfile?.SoilProfile;
RoundedDouble exitPointL = input.ExitPointL;
if (surfaceLine != null && soilProfile != null && !double.IsNaN(exitPointL))
{
return soilProfile.GetConsecutiveAquiferLayersBelowLevel(surfaceLine.GetZAtL(exitPointL)).ToArray();
}
return new MacroStabilityInwardsSoilLayer[0];
}
private MacroStabilityInwardsSoilLayer[] GetConsecutiveCoverageLayers()
{
RingtoetsMacroStabilityInwardsSurfaceLine surfaceLine = input.SurfaceLine;
MacroStabilityInwardsSoilProfile soilProfile = input.StochasticSoilProfile?.SoilProfile;
RoundedDouble exitPointL = input.ExitPointL;
if (surfaceLine != null && soilProfile != null && !double.IsNaN(exitPointL))
{
MacroStabilityInwardsSoilLayer[] consecutiveAquitardLayersBelowLevel = soilProfile
.GetConsecutiveCoverageLayersBelowLevel(surfaceLine.GetZAtL(exitPointL))
.ToArray();
return consecutiveAquitardLayersBelowLevel;
}
return new MacroStabilityInwardsSoilLayer[0];
}
private static double GetThicknessTopAquiferLayer(MacroStabilityInwardsSoilProfile soilProfile, RingtoetsMacroStabilityInwardsSurfaceLine surfaceLine, RoundedDouble exitPointL)
{
try
{
double zAtL = surfaceLine.GetZAtL(exitPointL);
return soilProfile.GetTopmostConsecutiveAquiferLayerThicknessBelowLevel(zAtL);
}
catch (Exception e)
{
if (e is RingtoetsMacroStabilityInwardsSurfaceLineException || e is InvalidOperationException || e is ArgumentException)
{
return double.NaN;
}
throw;
}
}
private static double GetThicknessCoverageLayers(MacroStabilityInwardsSoilProfile soilProfile, RingtoetsMacroStabilityInwardsSurfaceLine surfaceLine, RoundedDouble exitPointL)
{
try
{
double zAtL = surfaceLine.GetZAtL(exitPointL);
return soilProfile.GetConsecutiveCoverageLayerThicknessBelowLevel(zAtL);
}
catch (Exception e)
{
if (e is RingtoetsMacroStabilityInwardsSurfaceLineException || e is InvalidOperationException || e is ArgumentException)
{
return double.NaN;
}
throw;
}
}
}
}