Load control analyzer with trial steps for residual evaluation with minor changes

src/MGroup.NumericalAnalyzers.Discretization/NonLinear/ArcLengthAnalyzer.cs

@@ -55,7 +55,7 @@ public class ArcLengthAnalyzer : NonLinearAnalyzerBase
  /// <param name="shape">Option for the shape of the constraint - 0 for cylindrical, 1 for spherical, intermediate values for elliptic (default : shape = 0)</param>
  /// <param name="constConstraint">Option for constant radius of the constraint (default : constConstraint = 'true')</param>
  /// <param name="numOfIterations">(only usefull for constConstraint = false) Number of expected iterations within a load increment (default : numOfIterations = 4)</param>
- private ArcLengthAnalyzer(IAlgebraicModel algebraicModel, ISolver solver, INonLinearProvider provider, int numIncrements, int maxIterationsPerIncrement, 
+ private ArcLengthAnalyzer(IAlgebraicModel algebraicModel, ISolver solver, INonLinearProvider provider, int numIncrements, int maxIterationsPerIncrement,
  int numIterationsForMatrixRebuild, double residualTolerance, double shape, bool constConstraint, int numOfIterations, bool stopIfNotConverged)
  : base(algebraicModel, solver, provider, numIncrements, maxIterationsPerIncrement,
  numIterationsForMatrixRebuild, residualTolerance)
@@ -323,7 +323,7 @@ public class Builder : NonLinearAnalyzerBuilderBase
  private int numOfIterations;
  private bool stopIfNotConverged = true;
 
- public Builder(IAlgebraicModel algebraicModel, ISolver solver, INonLinearProvider provider, int numIncrements, double shape = 0, int numOfIterations = 4, 
+ public Builder(IAlgebraicModel algebraicModel, ISolver solver, INonLinearProvider provider, int numIncrements, double shape = 0, int numOfIterations = 4,
  bool constConstraint = true, bool stopIfNotConverged = true)
  : base(algebraicModel, solver, provider, numIncrements)
  {
@@ -336,7 +336,7 @@ public Builder(IAlgebraicModel algebraicModel, ISolver solver, INonLinearProvide
  this.stopIfNotConverged = stopIfNotConverged;
  }
 
- public ArcLengthAnalyzer Build() => new ArcLengthAnalyzer(algebraicModel, solver, provider, 
+ public ArcLengthAnalyzer Build() => new ArcLengthAnalyzer(algebraicModel, solver, provider,
  numIncrements, maxIterationsPerIncrement, numIterationsForMatrixRebuild, residualTolerance, shape, constConstraint, numOfIterations, stopIfNotConverged);
  }
  }

src/MGroup.NumericalAnalyzers.Discretization/NonLinear/LoadControlAnalyzer.cs

@@ -7,7 +7,7 @@
 using MGroup.MSolve.Solution.AlgebraicModel;
 using MGroup.NumericalAnalyzers.Logging;
 using MGroup.NumericalAnalyzers.NonLinear;
-
+using MGroup.LinearAlgebra;
 
 namespace MGroup.NumericalAnalyzers.Discretization.NonLinear
 {

src/MGroup.NumericalAnalyzers.Discretization/NonLinear/LoadControlAnalyzerWithTrialsForResidualEvaluation.cs

@@ -0,0 +1,140 @@
+using System;
+using System.Diagnostics;
+using MGroup.MSolve.AnalysisWorkflow;
+using MGroup.MSolve.AnalysisWorkflow.Providers;
+using MGroup.MSolve.Solution;
+using MGroup.MSolve.Solution.LinearSystem;
+using MGroup.MSolve.Solution.AlgebraicModel;
+using MGroup.NumericalAnalyzers.Logging;
+using MGroup.NumericalAnalyzers.NonLinear;
+using MGroup.LinearAlgebra;
+
+namespace MGroup.NumericalAnalyzers.Discretization.NonLinear
+{
+ public class LoadControlAnalyzerWithTrialsForResidualEvaluation : NonLinearAnalyzerWithTrialsBase
+ {
+ /// <summary>
+ /// This class solves the linearized geoemtrically nonlinear system of equations according to Newton-Raphson's load control incremental-iterative method.
+ /// </summary>
+ /// <param name="model">Instance of the model that will be solved</param>
+ /// <param name="solver">Instance of the solver that will solve the linear system of equations</param>
+ /// <param name="provider">Instance of the problem type to be solved</param>
+ /// <param name="subdomainUpdaters">Instance that updates constraints, right-hand-side vector, updates and resets state</param>
+ /// <param name="numIncrements">Number of total load increments</param>
+ /// <param name="maxIterationsPerIncrement">Number of maximum iterations within a load increment</param>
+ /// <param name="numIterationsForMatrixRebuild">Number of iterations for the rebuild of the siffness matrix within a load increment</param>
+ /// <param name="residualTolerance">Tolerance for the convergence criterion of the residual forces</param>
+ private LoadControlAnalyzerWithTrialsForResidualEvaluation(IAlgebraicModel algebraicModel, ISolver solver, INonLinearProvider provider,
+ int numIncrements, int maxIterationsPerIncrement, int numIterationsForMatrixRebuild, double residualTolerance)
+ : base(algebraicModel, solver, provider, numIncrements, maxIterationsPerIncrement,
+ numIterationsForMatrixRebuild, residualTolerance)
+ { }
+
+
+ /// <summary>
+ /// Solves the nonlinear equations and calculates the displacements vector.
+ /// </summary>
+ public override void Solve()
+ {
+ InitializeLogs();
+
+ DateTime start = DateTime.Now;
+ UpdateInternalVectors();
+ for (int increment = 0; increment < numIncrements; increment++)
+ {
+ double errorNorm = 0;
+ ClearIncrementalSolutionVector();
+ UpdateRhs(increment);
+
+ double firstError = 0;
+ int iteration = 0;
+ for (iteration = 0; iteration < maxIterationsPerIncrement; iteration++)
+ {
+ if (iteration == maxIterationsPerIncrement - 1)
+ {
+ return;
+ }
+
+ if (double.IsNaN(errorNorm))
+ {
+ return;
+ }
+
+ solver.Solve();
+ StoreSystemSolutionState();
+ int nTrials = 4; int ratio = 4; double[] trialErrors = new double[nTrials]; double trialSize = (double)1 / ratio;
+ for (int i1 = 0; i1 < nTrials; i1++)
+ {
+ var trialInternalRhsVector = CalculateInternalTrialRHS(increment, iteration, (i1 + 1) * trialSize);
+ trialErrors[i1] = globalRhsNormInitial != 0 ? UpdateResidualForcesAndNormForTrials(increment, trialInternalRhsVector) / globalRhsNormInitial : 0;
+ RestoreSystemSolutionState();
+ }
+ int minErrorTrialNo = 0;
+ double minErrorTrialValue = trialErrors[0];
+ for (int i1 = 1; i1 < nTrials; i1++)
+ {
+ if (trialErrors[i1] < minErrorTrialValue)
+ {
+ minErrorTrialNo = i1;
+ minErrorTrialValue = trialErrors[i1];
+ }
+ }
+
+ Console.WriteLine($"chosen trial no is {minErrorTrialNo + 1}"); Debug.WriteLine($"chosen trial no is {minErrorTrialNo + 1}");
+
+ var internalRhsVector = CalculateInternalTrialRHS(increment, iteration, (minErrorTrialNo + 1) * trialSize);
+ double residualNormCurrent = UpdateResidualForcesAndNormForTrials(increment, internalRhsVector);
+ errorNorm = globalRhsNormInitial != 0 ? residualNormCurrent / globalRhsNormInitial : 0;
+
+ if (iteration == 0)
+ {
+ firstError = errorNorm;
+ }
+
+ if (IncrementalDisplacementsLog != null)
+ {
+ IncrementalDisplacementsLog.StoreDisplacements(uPlusdu);
+ }
+
+ if (TotalDisplacementsPerIterationLog != null)
+ {
+ TotalDisplacementsPerIterationLog.StoreDisplacements(uPlusdu);
+ }
+
+ if (errorNorm < residualTolerance)
+ {
+ if (IncrementalLog != null)
+ {
+ IncrementalLog.LogTotalDataForIncrement(increment, iteration, errorNorm, uPlusdu, internalRhsVector);
+ }
+ break;
+ }
+
+ if ((iteration + 1) % numIterationsForMatrixRebuild == 0)
+ {
+ provider.Reset();
+ parentAnalyzer.BuildMatrices();
+ }
+ }
+ Debug.WriteLine("NR {0}, first error: {1}, exit error: {2}", iteration, firstError, errorNorm);
+ SaveMaterialStateAndUpdateSolution();
+ }
+ DateTime end = DateTime.Now;
+ StoreLogResults(start, end);
+ }
+
+ public class Builder : NonLinearAnalyzerBuilderBase
+ {
+ public Builder(IAlgebraicModel algebraicModel, ISolver solver, INonLinearProvider provider, int numIncrements)
+ : base(algebraicModel, solver, provider, numIncrements)
+ {
+ MaxIterationsPerIncrement = 1000;
+ NumIterationsForMatrixRebuild = 1;
+ ResidualTolerance = 1E-3;
+ }
+
+ public LoadControlAnalyzerWithTrialsForResidualEvaluation Build() => new LoadControlAnalyzerWithTrialsForResidualEvaluation(algebraicModel, solver, provider,
+ numIncrements, maxIterationsPerIncrement, numIterationsForMatrixRebuild, residualTolerance);
+ }
+ }
+}

src/MGroup.NumericalAnalyzers/Dynamic/BDFDynamicAnalyzer.cs

@@ -253,7 +253,7 @@ private void SolveCurrentTimestep()
  ChildAnalyzer.Solve();
  analysisStatistics[currentStep] = ChildAnalyzer.AnalysisStatistics;
  end = DateTime.Now;
- Debug.WriteLine("BDF elapsed time: {0}", end-start);
+ Debug.WriteLine("BDF elapsed time: {0}", end - start);
  }
 
  /// <summary>

src/MGroup.NumericalAnalyzers/LinearAnalyzer.cs

@@ -1,5 +1,6 @@
 using System;
 using System.Collections.Generic;
+
 using DotNumerics.ODE.Radau5;
 
 using MGroup.LinearAlgebra.Iterative;
@@ -53,8 +54,8 @@ public LinearAnalyzer(IAlgebraicModel algebraicModel, ISolver solver, IAnalyzerP
 
  public IGlobalVector CurrentAnalysisLinearSystemRhs { get => solver.LinearSystem.RhsVector; }
 
- public IterativeStatistics AnalysisStatistics => analysisStatistics; 
- 
+ public IterativeStatistics AnalysisStatistics => analysisStatistics;
+
  public void Initialize(bool isFirstAnalysis)
  {
  //if (isFirstAnalysis)
@@ -67,7 +68,7 @@ public void Solve()
  {
  var start = DateTime.Now;
  solver.Solve();
- 
+
  analysisStatistics.HasConverged = true;
  analysisStatistics.NumIterationsRequired = 1;
  Responses = solver.LinearSystem.Solution.Copy();
@@ -81,7 +82,7 @@ private void AddEquivalentNodalLoadsToRHS()
  //TODO: equivalentNodalLoads will often be 0. Perhaps instead of AddToGlobalVector, we should have AxpyToGlobalVector
  IGlobalVector equivalentNodalLoads = algebraicModel.CreateZeroVector();
  algebraicModel.AddToGlobalVector(provider.EnumerateEquivalentNeumannBoundaryConditions, equivalentNodalLoads);
- solver.LinearSystem.RhsVector.SubtractIntoThis(equivalentNodalLoads); 
+ solver.LinearSystem.RhsVector.SubtractIntoThis(equivalentNodalLoads);
  }
 
  private void InitializeLogs()

src/MGroup.NumericalAnalyzers/NonLinear/NonLinearAnalyzerBase.cs

@@ -75,7 +75,7 @@ public IParentAnalyzer ParentAnalyzer
  public IGlobalVector CurrentAnalysisLinearSystemRhs { get => solver.LinearSystem.RhsVector; }
 
  public IterativeStatistics AnalysisStatistics => analysisStatistics;
- 
+
  GenericAnalyzerState IAnalyzer.CurrentState
  {
  get => currentState;