From 9b1fba6e3a83c3dced8eb6fd72b7af6044e176e4 Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Fri, 31 Mar 2023 13:15:20 +0200
Subject: [PATCH] Add scaling to Riemannian PN-Solver and use it in
cosserat-continuum and film-on-substrate
---
dune/gfe/riemannianpnsolver.cc | 2 +-
dune/gfe/riemannianpnsolver.hh | 12 +++++++++++-
problems/cosserat-continuum-cantilever.parset | 2 +-
problems/cosserat-continuum-twisted-strip.parset | 2 +-
...sserat-continuum-wong-pellegrino-wrinkling.parset | 2 ++
problems/cosserat-continuum-wriggers-l-shape.parset | 4 ++--
problems/film-on-substrate.parset | 2 +-
problems/gradient-flow-curve-shortening.parset | 2 ++
problems/harmonicmaps-skyrmions-hexagon.parset | 2 ++
problems/simofox-cantilever.parset | 2 +-
problems/staticrod.parset | 2 ++
src/cosserat-continuum.cc | 10 ++++++----
src/film-on-substrate.cc | 5 +++--
13 files changed, 35 insertions(+), 14 deletions(-)
diff --git a/dune/gfe/riemannianpnsolver.cc b/dune/gfe/riemannianpnsolver.cc
index a10252d6..480c1ae9 100644
--- a/dune/gfe/riemannianpnsolver.cc
+++ b/dune/gfe/riemannianpnsolver.cc
@@ -253,7 +253,7 @@ void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
// Add the regularization - Identity Matrix for now
for (std::size_t i=0; i<stiffnessMatrix.N(); i++)
for(int j=0; j<blocksize; j++)
- stiffnessMatrix[i][i][j][j] += regularization;
+ stiffnessMatrix[i][i][j][j] += regularization*scaling_[j];
innerSolver_->setProblem(stiffnessMatrix,corr_global,rhs_global);
innerSolver_->preprocess();
diff --git a/dune/gfe/riemannianpnsolver.hh b/dune/gfe/riemannianpnsolver.hh
index b3584fdf..7ea10fe6 100644
--- a/dune/gfe/riemannianpnsolver.hh
+++ b/dune/gfe/riemannianpnsolver.hh
@@ -69,7 +69,9 @@ public:
RiemannianProximalNewtonSolver()
: IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
hessianMatrix_(nullptr), h1SemiNorm_(NULL)
- {}
+ {
+ std::fill(scaling_.begin(), scaling_.end(), 1.0);
+ }
/** \brief Set up the solver using a choldmod or umfpack solver as the inner solver */
void setup(const GridType& grid,
@@ -81,6 +83,11 @@ public:
double initialRegularization,
bool instrumented);
+ void setScaling(const Dune::FieldVector<double,blocksize>& scaling)
+ {
+ scaling_ = scaling;
+ }
+
void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
{
ignoreNodes_ = &ignoreNodes;
@@ -118,6 +125,9 @@ protected:
double initialRegularization_;
double tolerance_;
+ /** \brief Regularization scaling */
+ Dune::FieldVector<double,blocksize> scaling_;
+
/** \brief Maximum number of proximal-newton steps */
std::size_t maxProximalNewtonSteps_;
diff --git a/problems/cosserat-continuum-cantilever.parset b/problems/cosserat-continuum-cantilever.parset
index 9440301a..bfa21a84 100644
--- a/problems/cosserat-continuum-cantilever.parset
+++ b/problems/cosserat-continuum-cantilever.parset
@@ -26,7 +26,7 @@ tolerance = 1e-3
# Max number of steps of the trust region solver
maxSolverSteps = 200
-trustRegionScaling = 1 1 1 0.01 0.01 0.01
+solverScaling = 1 1 1 0.01 0.01 0.01
# Initial trust-region radius
initialTrustRegionRadius = 3.125
diff --git a/problems/cosserat-continuum-twisted-strip.parset b/problems/cosserat-continuum-twisted-strip.parset
index 04c738f4..a10d3197 100644
--- a/problems/cosserat-continuum-twisted-strip.parset
+++ b/problems/cosserat-continuum-twisted-strip.parset
@@ -23,7 +23,7 @@ tolerance = 1e-8
# Max number of steps of the trust region solver
maxTrustRegionSteps = 200
-trustRegionScaling = 1 1 1 1 1 1
+solverScaling = 1 1 1 1 1 1
# Initial trust-region radius
initialTrustRegionRadius = 0.1
diff --git a/problems/cosserat-continuum-wong-pellegrino-wrinkling.parset b/problems/cosserat-continuum-wong-pellegrino-wrinkling.parset
index 790a91cd..4af4a38d 100644
--- a/problems/cosserat-continuum-wong-pellegrino-wrinkling.parset
+++ b/problems/cosserat-continuum-wong-pellegrino-wrinkling.parset
@@ -23,6 +23,8 @@ tolerance = 1e-8
# Max number of steps of the trust region solver
maxTrustRegionSteps = 500
+solverScaling = 1 1 1 0.01 0.01 0.01
+
# Initial trust-region radius
initialTrustRegionRadius = 0.001
diff --git a/problems/cosserat-continuum-wriggers-l-shape.parset b/problems/cosserat-continuum-wriggers-l-shape.parset
index 12c715bd..745f3336 100644
--- a/problems/cosserat-continuum-wriggers-l-shape.parset
+++ b/problems/cosserat-continuum-wriggers-l-shape.parset
@@ -22,7 +22,7 @@ tolerance = 1e-8
# Max number of steps of the trust region solver
maxSolverSteps = 1000
-trustRegionScaling = 1 1 1 0.01 0.01 0.01
+solverScaling = 1 1 1 0.01 0.01 0.01
# Initial trust-region radius
initialTrustRegionRadius = 1
@@ -105,4 +105,4 @@ neumannValues = 0.09 0 0
startFromFile = yes
initialIterateGridFilename = wriggers-L-shape_99_mm.msh
-initialIterateFilename = initial-wriggers-l-shape.vtu
\ No newline at end of file
+initialIterateFilename = initial-wriggers-l-shape.vtu
diff --git a/problems/film-on-substrate.parset b/problems/film-on-substrate.parset
index e6cf6ce9..1ce5db4b 100644
--- a/problems/film-on-substrate.parset
+++ b/problems/film-on-substrate.parset
@@ -87,7 +87,7 @@ initialRegularization = 1000000
# Solver parameters specific for trust-region solver using multigrid solver
#############################################
-trustRegionScaling = 1 1 1 0.01 0.01 0.01
+solverScaling = 1 1 1 0.01 0.01 0.01
# Initial trust-region radius
initialTrustRegionRadius = 3.125
diff --git a/problems/gradient-flow-curve-shortening.parset b/problems/gradient-flow-curve-shortening.parset
index cf5f31f4..857ec53c 100644
--- a/problems/gradient-flow-curve-shortening.parset
+++ b/problems/gradient-flow-curve-shortening.parset
@@ -20,6 +20,8 @@ tolerance = 1e-8
# Max number of steps of the trust region solver
maxTrustRegionSteps = 100
+solverScaling = 1 1 1 0.01 0.01 0.01
+
# Initial trust-region radius
initialTrustRegionRadius = 0.25
diff --git a/problems/harmonicmaps-skyrmions-hexagon.parset b/problems/harmonicmaps-skyrmions-hexagon.parset
index 1394c0ec..4cd7c41b 100644
--- a/problems/harmonicmaps-skyrmions-hexagon.parset
+++ b/problems/harmonicmaps-skyrmions-hexagon.parset
@@ -21,6 +21,8 @@ tolerance = 1e-12
# Max number of steps of the trust region solver
maxTrustRegionSteps = 100
+solverScaling = 1 1 1 0.01 0.01 0.01
+
# Initial trust-region radius
initialTrustRegionRadius = 1
diff --git a/problems/simofox-cantilever.parset b/problems/simofox-cantilever.parset
index 1eab7e5b..4351b50e 100644
--- a/problems/simofox-cantilever.parset
+++ b/problems/simofox-cantilever.parset
@@ -26,7 +26,7 @@ tolerance = 1e-6
# Max number of steps of the trust region solver
maxTrustRegionSteps = 200
-trustRegionScaling = 1 1 1 0.01 0.01
+solverScaling = 1 1 1 0.01 0.01
# Initial trust-region radius
initialTrustRegionRadius = 10
diff --git a/problems/staticrod.parset b/problems/staticrod.parset
index 2d001eca..7f8e5fe5 100644
--- a/problems/staticrod.parset
+++ b/problems/staticrod.parset
@@ -21,6 +21,8 @@ tolerance = 1e-6
# Max number of steps of the trust region solver
maxTrustRegionSteps = 100
+solverScaling = 1 1 1 0.01 0.01 0.01
+
# Initial trust-region radius
initialTrustRegionRadius = 1
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 17333187..23941144 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -477,7 +477,7 @@ int main (int argc, char *argv[]) try
baseTolerance,
instrumented);
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(x);
solver.solve();
@@ -516,7 +516,7 @@ int main (int argc, char *argv[]) try
baseTolerance,
instrumented);
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
@@ -530,6 +530,7 @@ int main (int argc, char *argv[]) try
maxSolverSteps,
initialRegularization,
instrumented);
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
@@ -577,7 +578,7 @@ int main (int argc, char *argv[]) try
baseTolerance,
instrumented);
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(x);
solver.solve();
@@ -616,7 +617,7 @@ int main (int argc, char *argv[]) try
baseTolerance,
instrumented);
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
@@ -630,6 +631,7 @@ int main (int argc, char *argv[]) try
maxSolverSteps,
initialRegularization,
instrumented);
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 8978172a..2dd06845 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -584,7 +584,7 @@ int main (int argc, char *argv[]) try
baseTolerance,
instrumented);
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(x);
solver.solve();
x = solver.getSol();
@@ -604,7 +604,7 @@ int main (int argc, char *argv[]) try
baseTolerance,
instrumented);
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xRBM);
solver.solve();
xRBM = solver.getSol();
@@ -627,6 +627,7 @@ int main (int argc, char *argv[]) try
maxSolverSteps,
initialRegularization,
instrumented);
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xRBM);
solver.solve();
xRBM = solver.getSol();
--
GitLab