From c8ed9a4b806fa8a1424dff602a9560d27ab7f8b7 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Fri, 16 Apr 2021 18:50:16 +0200
Subject: [PATCH] .
---
data/tools/config.ini | 4 ++--
data/tools/main.py | 4 ++--
data/tools/support/friction_stats.py | 18 +++++++++---------
data/tools/support/iterations.py | 10 +++++-----
.../friction/frictionpotential.hh | 4 ++--
.../factories/stackedblocksfactory.cc | 19 ++++++++++++-------
dune/tectonic/factories/strikeslipfactory.cc | 6 +++---
.../multi-body-problem-2D.cfg | 2 +-
src/multi-body-problem/multi-body-problem.cc | 2 +-
src/multi-body-problem/multi-body-problem.cfg | 9 ++++++---
src/strikeslip/strikeslip.cfg | 6 +++---
11 files changed, 46 insertions(+), 38 deletions(-)
diff --git a/data/tools/config.ini b/data/tools/config.ini
index 79c1ec79..0d5e4320 100644
--- a/data/tools/config.ini
+++ b/data/tools/config.ini
@@ -1,9 +1,9 @@
-#/home/joscha/software/dune/build-release/dune-tectonic/src/foam/output/test/
+# /home/joscha/software/dune/build-release/dune-tectonic/src/foam/output/pipping-2013-euler/
#/home/joscha/software/dune/build-debug/dune-tectonic/src/foam/output/test/
#/home/joscha/Desktop/strikeslip/viscosity-1e3/
#/home/joscha/Downloads/
[directories]
-simulation = /home/joscha/software/dune/build-release/dune-tectonic/src/foam/output/pipping-2013-euler/
+simulation = /home/joscha/software/dune/build-release/dune-tectonic/src/multi-body-problem/output/newmark/
experiment = ~/group/publications/2016-RosenauCorbiDominguezRudolfRitterPipping
output = generated
diff --git a/data/tools/main.py b/data/tools/main.py
index d18a1fa7..31aabba6 100644
--- a/data/tools/main.py
+++ b/data/tools/main.py
@@ -36,7 +36,7 @@ def build_patch(coords, percentage):
return [i for i in range(len(x_coords)) if x_coords[i]>=xmin and x_coords[i]<=xmax]
# read problem parameters
-params = read_params('foam.cfg')
+params = read_params('multi-body-problem.cfg')
TANGENTIAL_COORDS = 1
FINAL_TIME = params['finalTime']
@@ -45,7 +45,7 @@ NBODIES = params['bodyCount']
h5file = read_h5file()
print(list(h5file.keys()))
-interval = [23, 24] #[0.75*FINAL_TIME, FINAL_TIME]
+interval = [0, FINAL_TIME] #[0.75*FINAL_TIME, FINAL_TIME]
iterations(h5file, FINAL_TIME, interval)
diff --git a/data/tools/support/friction_stats.py b/data/tools/support/friction_stats.py
index 0ee1dea6..7f0a13c4 100644
--- a/data/tools/support/friction_stats.py
+++ b/data/tools/support/friction_stats.py
@@ -27,7 +27,7 @@ def friction_stats(h5file, body_ID, FINAL_TIME, patch=[], interval=[], TANGENTIA
min_v = np.min(v_tx, axis=1)
max_v = np.max(v_tx, axis=1)
# plot
- ax_slip = fig.add_subplot(2, 1, 1)
+ ax_slip = fig.add_subplot(1, 1, 1)
#ax_slip.plot(time, min_v, color='gray', linestyle='--')
ax_slip.plot(time, avg_v, color='black', linestyle='-')
#ax_slip.plot(time, max_v, color='gray', linestyle='--')
@@ -41,20 +41,20 @@ def friction_stats(h5file, body_ID, FINAL_TIME, patch=[], interval=[], TANGENTIA
print(np.max(max_v))
# state
- states = np.array(h5file[body + '/state'])
- states_t = states[t,:]
- states_tx = states_t[:,patch]
+ #states = np.array(h5file[body + '/state'])
+ #states_t = states[t,:]
+ #states_tx = states_t[:,patch]
# statistics
- avg_states = np.average(states_tx, axis=1)
+ #avg_states = np.average(states_tx, axis=1)
#min_states = np.min(states_tx, axis=1)
#max_states = np.max(states_tx, axis=1)
# plot
- ax_state = fig.add_subplot(2, 1, 2)
+ #ax_state = fig.add_subplot(2, 1, 2)
#ax_state.plot(time, min_states, color='gray', linestyle='--')
- ax_state.plot(time, avg_states, color='black', linestyle='-')
+ #ax_state.plot(time, avg_states, color='black', linestyle='-')
#ax_state.plot(time, max_states, color='gray', linestyle='--')
- ax_state.set_ylabel('state')
- ax_state.set_xlabel('time [s]')
+ #ax_state.set_ylabel('state')
+ #ax_state.set_xlabel('time [s]')
#-------------------------
# friction coefficient
diff --git a/data/tools/support/iterations.py b/data/tools/support/iterations.py
index 1e231bb9..adf0550a 100644
--- a/data/tools/support/iterations.py
+++ b/data/tools/support/iterations.py
@@ -27,20 +27,20 @@ def iterations(h5file, FINAL_TIME, interval = []):
# plot
fig = plt.figure()
- ax_fpi = fig.add_subplot(2, 1, 1)
+ ax_fpi = fig.add_subplot(3, 1, 1)
ax_fpi.plot(time, fpi_final[t], color='black', linestyle='-')
ax_fpi.set_ylabel('fpi')
#-------------------------
- ax_mg = fig.add_subplot(2, 1, 2)
+ ax_mg = fig.add_subplot(3, 1, 2)
ax_mg.plot(time, multigrid_final[t], color='black', linestyle='-')
ax_mg.set_ylabel('multigrid iter.')
ax_mg.set_xlabel('time [s]')
#-------------------------
- #ax_tau = fig.add_subplot(3, 1, 3)
- #ax_tau.plot(time, tau[t], color='black', linestyle='-')
- #ax_tau.set_ylabel('tau')
+ ax_tau = fig.add_subplot(3, 1, 3)
+ ax_tau.plot(time, tau[t], color='black', linestyle='-')
+ ax_tau.set_ylabel('tau')
#-------------------------
fig.canvas.draw()
\ No newline at end of file
diff --git a/dune/tectonic/data-structures/friction/frictionpotential.hh b/dune/tectonic/data-structures/friction/frictionpotential.hh
index 8e742ea7..0ab5c5eb 100644
--- a/dune/tectonic/data-structures/friction/frictionpotential.hh
+++ b/dune/tectonic/data-structures/friction/frictionpotential.hh
@@ -49,8 +49,8 @@ class TruncatedRateState : public FrictionPotential {
double differential(double V) const override {
//std::cout << "differential: " << weight * fd.C - weightedNormalStress * coefficientOfFriction(V) << std::endl;
- //if (V <= Vmin or regularity(V)>10e8)
- // return 0.0;
+ if (V <= Vmin)
+ return 0.0;
return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
}
diff --git a/dune/tectonic/factories/stackedblocksfactory.cc b/dune/tectonic/factories/stackedblocksfactory.cc
index 877e32cf..bd32036f 100644
--- a/dune/tectonic/factories/stackedblocksfactory.cc
+++ b/dune/tectonic/factories/stackedblocksfactory.cc
@@ -19,9 +19,12 @@ template <class HostGridType, class VectorType>
void StackedBlocksFactory<HostGridType, VectorType>::setBodies() {
// set up cuboid geometries
- double const length = 1.00;
- double const height = 0.27;
- double const weakLength = 0.60;
+ const auto& bodyParset = this->parset_.sub("body");
+
+ double const length = bodyParset.template get<double>("length");
+ double const height = bodyParset.template get<double>("height");
+ double const heightScaling = bodyParset.template get<double>("heightScaling");
+ double const weakLength = 1.00;
std::vector<GlobalCoords> origins(this->bodyCount_);
@@ -60,20 +63,22 @@ void StackedBlocksFactory<HostGridType, VectorType>::setBodies() {
cuboidGeometries_[idx]->addWeakeningPatch(subParset, lowerWeakOrigin, weakBound(lowerWeakOrigin));
#elif MY_DIM == 2
auto weakBound = [&] (const GlobalCoords& origin) {
- GlobalCoords h = {origin[0] + weakLength, origin[1]};
+ GlobalCoords h = {origin[0] + weakLength*length, origin[1]};
return h;
};
origins[0] = {0,0};
- GlobalCoords lowerWeakOrigin = {0.2, 0};
- GlobalCoords upperWeakOrigin = {0.2, height};
+ GlobalCoords lowerWeakOrigin = {0.0, 0};
+ GlobalCoords upperWeakOrigin = {0.0, height};
cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], length, height);
cuboidGeometries_[0]->addWeakeningPatch(subParset, upperWeakOrigin, weakBound(upperWeakOrigin));
for (size_t i=1; i<this->bodyCount_-1; i++) {
origins[i] = cuboidGeometries_[i-1]->upperLeft();
- auto height_i = height/3.0;
+
+ double power = (i<this->bodyCount_/2.0) ? i : (this->bodyCount_-i-1);
+ auto height_i = height*std::pow(heightScaling, power);
GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
GlobalCoords upperWeakOrigin_i = {upperWeakOrigin[0], height_i};
diff --git a/dune/tectonic/factories/strikeslipfactory.cc b/dune/tectonic/factories/strikeslipfactory.cc
index d10478cd..9449ff95 100644
--- a/dune/tectonic/factories/strikeslipfactory.cc
+++ b/dune/tectonic/factories/strikeslipfactory.cc
@@ -78,10 +78,10 @@ void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setBodies() {
bodyData_[1] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body1"), 0.0, zenith_());
this->bodies_[1] = std::make_shared<typename Base::LeafBody>(bodyData_[1], grids[1]);
#else
- bodyData_[0] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body0"), this->parset_.template get<double>("gravity"), zenith_());
+ bodyData_[0] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body0"), this->parset_.template get<double>("general.gravity"), zenith_());
this->bodies_[0] = std::make_shared<typename Base::LeafBody>(bodyData_[0], grids[0]);
- bodyData_[1] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body1"), this->parset_.template get<double>("gravity"), zenith_());
+ bodyData_[1] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body1"), this->parset_.template get<double>("general.gravity"), zenith_());
this->bodies_[1] = std::make_shared<typename Base::LeafBody>(bodyData_[1], grids[1]);
#endif
}
@@ -124,7 +124,7 @@ void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
using Function = Dune::VirtualFunction<double, double>;
std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
- std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>(-1.0*this->parset_.template get<double>("boundary.dirichlet.finalVelocity"));
+ std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>(-1.0*this->parset_.template get<double>("boundary.dirichlet.finalVelocity"), 0.25);
const double lengthScale = TriangleGeometry::lengthScale();
diff --git a/src/multi-body-problem/multi-body-problem-2D.cfg b/src/multi-body-problem/multi-body-problem-2D.cfg
index f6110f74..a47d46dc 100644
--- a/src/multi-body-problem/multi-body-problem-2D.cfg
+++ b/src/multi-body-problem/multi-body-problem-2D.cfg
@@ -1,6 +1,6 @@
# -*- mode:conf -*-
[boundary.friction]
-smallestDiameter = 2e-2 # 0.05 2e-3 [m]
+smallestDiameter = 6.25e-2 # 0.05 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-5 # 1e-5
diff --git a/src/multi-body-problem/multi-body-problem.cc b/src/multi-body-problem/multi-body-problem.cc
index 58da326b..b4e40038 100644
--- a/src/multi-body-problem/multi-body-problem.cc
+++ b/src/multi-body-problem/multi-body-problem.cc
@@ -360,7 +360,7 @@ int main(int argc, char *argv[]) {
/*UniformTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
timeStepper(stepBase, contactNetwork, current,
- programState.relativeTime, programState.relativeTau); */
+ programState.relativeTime, programState.relativeTau);*/
AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
timeStepper(stepBase, contactNetwork, current,
diff --git a/src/multi-body-problem/multi-body-problem.cfg b/src/multi-body-problem/multi-body-problem.cfg
index 384b352c..a8160fb9 100644
--- a/src/multi-body-problem/multi-body-problem.cfg
+++ b/src/multi-body-problem/multi-body-problem.cfg
@@ -1,10 +1,13 @@
# -*- mode:conf -*-
[general]
-outPath = newmark-1e4 # output written to ./output/outPath
+outPath = newmark # output written to ./output/outPath
gravity = 9.81 # [m/s^2]
[body]
+length = 5.0 # [m]
+height = 1.0 # [m]
+heightScaling = 0.3
bulkModulus = 4.12e7 # [Pa]
poissonRatio = 0.3 # [1]
[body.elastic]
@@ -49,11 +52,11 @@ relativeTau = 1e-4 # 1e-6
[timeSteps]
scheme = newmark
-timeSteps = 1e4
+timeSteps = 1e6
[problem]
-finalTime = 50 # [s] #1000
+finalTime = 100 # [s] #1000
bodyCount = 4
diff --git a/src/strikeslip/strikeslip.cfg b/src/strikeslip/strikeslip.cfg
index f373d429..3ae01416 100644
--- a/src/strikeslip/strikeslip.cfg
+++ b/src/strikeslip/strikeslip.cfg
@@ -1,5 +1,5 @@
[general]
-outPath = friction # output written to ./output/strikeslip/outPath
+outPath = newmark-uniform-1e4 # output written to ./output/strikeslip/outPath
gravity = 0.0 # [m/s^2]
[body0]
@@ -35,7 +35,7 @@ shearViscosity = 0.0 # [Pas]
bulkViscosity = 0.0 # [Pas]
[boundary.friction]
-C = 6 # [Pa]
+C = 0 # [Pa]
mu0 = 0.48 # [ ]
V0 = 1e-3 # [m/s]
L = 1e-6 # [m]
@@ -64,7 +64,7 @@ restarts.write = true #true
vtk.write = true
[problem]
-finalTime = 15 # [s] #1000
+finalTime = 50 # [s] #1000
bodyCount = 2
[initialTime]
--
GitLab