From 90c97994213b4d9897c7adcbe1ed84f7c487af46 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 19 Jan 2015 17:01:37 +0100
Subject: [PATCH] [Noise] Clean up initial conditions
---
src/sand-wedge-data/boundaryconditions.py | 1 +
src/sand-wedge.cc | 80 ++++++++++-------------
2 files changed, 36 insertions(+), 45 deletions(-)
diff --git a/src/sand-wedge-data/boundaryconditions.py b/src/sand-wedge-data/boundaryconditions.py
index 05035c80..8db53fe0 100644
--- a/src/sand-wedge-data/boundaryconditions.py
+++ b/src/sand-wedge-data/boundaryconditions.py
@@ -6,6 +6,7 @@ class neumannCondition:
class velocityDirichletCondition:
def __call__(self, relativeTime):
+ # Assumed to vanish at time zero
finalVelocity = -5e-5
if relativeTime <= 0.1:
return finalVelocity * ( 1.0 - math.cos(relativeTime * math.pi / 0.1) ) / 2.0
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index 6615a7a1..6c75bb17 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -194,7 +194,6 @@ int main(int argc, char *argv[]) {
myAssembler.assembleViscosity(body.getShearViscosityField(),
body.getBulkViscosityField(), C);
myAssembler.assembleMass(body.getDensityField(), M);
- EnergyNorm<Matrix, Vector> const ANorm(A), MNorm(M);
ScalarMatrix frictionalBoundaryMass;
myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
@@ -213,21 +212,18 @@ int main(int argc, char *argv[]) {
_relativeTime);
_ell += gravityFunctional;
};
- Vector ell0(leafVertexCount);
- computeExternalForces(0.0, ell0);
-
- // {{{ Initial conditions
- using LinearFactory =
- SolverFactory<dims,
- BlockNonlinearTNNMGProblem<ConvexProblem<
- ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
- Grid>;
- ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
+ // helper
auto const solveLinearProblem = [&](
Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
- Vector const &_rhs, Vector &_x, EnergyNorm<Matrix, Vector> const &_norm,
+ Vector const &_rhs, Vector &_x,
Dune::ParameterTree const &_localParset) {
+
+ using LinearFactory = SolverFactory<
+ dims, BlockNonlinearTNNMGProblem<ConvexProblem<
+ ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
+ Grid>;
+ ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
refinements, *grid, _dirichletNodes);
@@ -237,9 +233,10 @@ int main(int argc, char *argv[]) {
auto multigridStep = factory.getStep();
multigridStep->setProblem(_x, problem);
+ EnergyNorm<Matrix, Vector> const norm(_matrix);
LoopSolver<Vector> solver(
multigridStep.get(), _localParset.get<size_t>("maximumIterations"),
- _localParset.get<double>("tolerance"), &_norm,
+ _localParset.get<double>("tolerance"), &norm,
_localParset.get<Solver::VerbosityMode>("verbosity"),
false); // absolute error
@@ -247,52 +244,45 @@ int main(int argc, char *argv[]) {
solver.solve();
};
- // Solve the stationary problem
+ // {{{ Initial conditions
+ Vector ell0(leafVertexCount);
+ computeExternalForces(0.0, ell0);
+
+ // Start from a situation of minimal stress, which is
+ // automatically attained in the case [v = 0 = a]. Assuming
+ // dPhi(v = 0) = 0, we thus only have to solve Au = ell0
Vector u_initial(leafVertexCount);
- u_initial = 0.0;
- solveLinearProblem(dirichletNodes, A, ell0, u_initial, ANorm,
+ solveLinearProblem(dirichletNodes, A, ell0, u_initial,
parset.sub("u0.solver"));
- ScalarVector alpha_initial(leafVertexCount);
- alpha_initial = parset.get<double>("boundary.friction.initialAlpha");
ScalarVector normalStress(leafVertexCount);
myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
body.getYoungModulus(),
body.getPoissonRatio(), u_initial);
- MyGlobalFrictionData<LocalVector> frictionInfo(
- parset.sub("boundary.friction"), weakPatch);
- auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
- parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
- frictionalBoundary, frictionInfo, normalStress);
- myGlobalFriction->updateAlpha(alpha_initial);
+ ScalarVector alpha_initial(leafVertexCount);
+ alpha_initial = parset.get<double>("boundary.friction.initialAlpha");
+
+ // Start from zero velocity
Vector v_initial(leafVertexCount);
v_initial = 0.0;
- {
- double v_initial_const;
- velocityDirichletFunction.evaluate(0.0, v_initial_const);
- assert(std::abs(v_initial_const) < 1e-14);
- }
+ // Compute the acceleration from Ma = ell0 - Au (without Dirichlet
+ // constraints), again assuming dPhi(v = 0) = 0
Vector a_initial(leafVertexCount);
- a_initial = 0.0;
- {
- // We solve Ma = ell0 - [Au + Cv + DPsi(v)]
- Vector accelerationRHS(leafVertexCount);
- {
- accelerationRHS = 0.0;
- Arithmetic::addProduct(accelerationRHS, A, u_initial);
- Arithmetic::addProduct(accelerationRHS, C, v_initial);
- // NOTE: We assume differentiability of Psi at 0 here!
- myGlobalFriction->addGradient(v_initial, accelerationRHS);
- accelerationRHS *= -1.0;
- accelerationRHS += ell0;
- }
- solveLinearProblem(noNodes, M, accelerationRHS, a_initial, MNorm,
- parset.sub("a0.solver"));
- }
+ Vector accelerationRHS = ell0;
+ Arithmetic::subtractProduct(accelerationRHS, A, u_initial);
+ solveLinearProblem(noNodes, M, accelerationRHS, a_initial,
+ parset.sub("a0.solver"));
// }}}
+ MyGlobalFrictionData<LocalVector> frictionInfo(
+ parset.sub("boundary.friction"), weakPatch);
+ auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
+ parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
+ frictionalBoundary, frictionInfo, normalStress);
+ myGlobalFriction->updateAlpha(alpha_initial);
+
Vector vertexCoordinates(leafVertexCount);
{
Dune::MultipleCodimMultipleGeomTypeMapper<
--
GitLab