diff --git a/src/Makefile.am b/src/Makefile.am index ac90649d2259744dce386774ee8dbd0726891303..7649f800aa2510de0e32509c4a3bace528975758 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -25,6 +25,7 @@ run-one-body-sample-gdb: one-body-sample one_body_sample_SOURCES = \ assemblers.cc \ compute_state.cc \ + compute_state_ruina.cc \ mysolver.cc \ one-body-sample.cc \ print.cc \ diff --git a/src/compute_state_ruina.cc b/src/compute_state_ruina.cc new file mode 100644 index 0000000000000000000000000000000000000000..ddba439cd1f68fcb36a951cba3be6424c35fff92 --- /dev/null +++ b/src/compute_state_ruina.cc @@ -0,0 +1,61 @@ +#include <cassert> + +#include <dune/common/fvector.hh> +#include <dune/common/fmatrix.hh> +#include <dune/fufem/interval.hh> + +#include <dune/tnnmg/problem-classes/bisection.hh> + +#include <dune/tectonic/mydirectionalconvexfunction.hh> + +#include "compute_state_ruina.hh" + +// Ruina's nonlinearity +class Ruina { +public: + typedef Dune::FieldVector<double, 1> VectorType; + typedef Dune::FieldMatrix<double, 1, 1> MatrixType; + + Ruina(double h, double uol) : h(h), uol(uol) {} + + void directionalSubDiff(VectorType const &u, VectorType const &v, + Interval<double> &D) const { + if (uol == 0) + D[0] = D[1] = v[0] * uol * (u[0]); + else + D[0] = D[1] = v[0] * uol * (u[0] + std::log(uol / h)); + } + + void directionalDomain(VectorType const &, VectorType const &, + Interval<double> &dom) const { + dom[0] = -std::numeric_limits<double>::max(); + dom[1] = std::numeric_limits<double>::max(); + } + +private: + double const h; + double const uol; +}; + +double state_update_ruina_bisection(double h, double uol, double old_state) { + Ruina::VectorType const start(0); + Ruina::VectorType const direction(1); + Ruina const phi(h, uol); + MyDirectionalConvexFunction<Ruina> const J(1.0, old_state, phi, start, + direction); + + int bisectionsteps = 0; + Bisection const bisection(0.0, 1.0, 1e-12, true, 0); // TODO + return bisection.minimize(J, 0.0, 0.0, bisectionsteps); // TODO +} + +double state_update_ruina(double h, double uol, // unorm/L + double old_state) { + double const uol_log_uol = (uol == 0) ? 0.0 : uol * std::log(uol); + + double ret1 = state_update_ruina_bisection(h, uol, old_state); + assert(std::abs(ret1 + uol * ret1 - old_state + uol_log_uol - + uol * std::log(h)) < 1e-14); + + return ret1; +} diff --git a/src/compute_state_ruina.hh b/src/compute_state_ruina.hh new file mode 100644 index 0000000000000000000000000000000000000000..43e8eba961ff564de9f27bf297cdbab7b2eaecdb --- /dev/null +++ b/src/compute_state_ruina.hh @@ -0,0 +1,6 @@ +#ifndef COMPUTE_STATE_RUINA_HH +#define COMPUTE_STATE_RUINA_HH + +double state_update_ruina(double h, double uol, double old_state); + +#endif diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index 991213dc062ada0d49e1ba96952a8243dd968458..306156131f388bcc90b848b408f2865a322200cd 100644 --- a/src/one-body-sample.cc +++ b/src/one-body-sample.cc @@ -63,6 +63,7 @@ #include "assemblers.hh" #include "compute_state.hh" +#include "compute_state_ruina.hh" #include "mysolver.hh" #include "print.hh" #include "vtk.hh" @@ -277,7 +278,14 @@ int main(int argc, char *argv[]) { // velocity // std::cout << std::log(L/unorm * h) << std::endl; - (*s4_new)[i] = state_update_dieterich(h, unorm / L, s4_old[i]); + auto const model = + parset.get<std::string>("boundary.friction.state.model"); + if (model == std::string("Dieterich")) + (*s4_new)[i] = state_update_dieterich(h, unorm / L, s4_old[i]); + else if (model == std::string("Ruina")) + (*s4_new)[i] = state_update_ruina(h, unorm / L, s4_old[i]); + else + assert(false); } } if (parset.get<bool>("printProgress")) { diff --git a/src/one-body-sample.parset b/src/one-body-sample.parset index c65b72443f813b42e7863038957146fd4f0362f4..9e2a31d9ecaa391d5166eee5bd34e4c924516d09 100644 --- a/src/one-body-sample.parset +++ b/src/one-body-sample.parset @@ -73,6 +73,7 @@ model = Ruina evolve = true # log(alpha(0)) initial = -2.42037 # Dirichlet: -4.72295 (1e-4), -2.42037 (1e-3), -0.117783 (1e-2) +model = Dieterich # Ruina [boundary.friction.ruina] # "For rocks, typical values of A and B range from 0.005 to 0.015"