diff --git a/src/compute_state.cc b/src/compute_state.cc index 4580c45903eec3c26e5f6842d4cca2417f141102..415ed209de60fccb5c904ead060b8570596ef9ad 100644 --- a/src/compute_state.cc +++ b/src/compute_state.cc @@ -10,7 +10,8 @@ #include "compute_state.hh" -// The nonlinearity exp(-x) +// The nonlinearity exp(-x) [ Note: the Dieterich law has an additional linear +// term! ] class DecayingExponential { public: typedef Dune::FieldVector<double, 1> VectorType; @@ -33,22 +34,23 @@ class DecayingExponential { double const h; }; -double compute_state_update_bisection(double h, double unorm, double L, - double old_state) { +double state_update_dieterich_bisection(double h, double uol, + double old_state) { DecayingExponential::VectorType const start(0); DecayingExponential::VectorType const direction(1); DecayingExponential const phi(h); MyDirectionalConvexFunction<DecayingExponential> const J( - 1.0, old_state - unorm / L, phi, start, direction); + 1.0, old_state - uol, 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 compute_state_update(double h, double unorm, double L, - double old_state) { - double ret1 = compute_state_update_bisection(h, unorm, L, old_state); - assert(std::abs(ret1 - old_state + unorm / L - h * std::exp(-ret1)) < 1e-12); +double state_update_dieterich(double h, double uol, // unorm / L + double old_state) { + double const ret1 = state_update_dieterich_bisection(h, uol, old_state); + assert(std::abs(ret1 - old_state + uol - h * std::exp(-ret1)) < 1e-12); return ret1; } diff --git a/src/compute_state.hh b/src/compute_state.hh index 6b4e8778affcf364b5a3812b8087800c1ff02ffc..b1190cc071ee577f425ffcb4ab87a08241fb45fc 100644 --- a/src/compute_state.hh +++ b/src/compute_state.hh @@ -1,6 +1,6 @@ #ifndef COMPUTE_STATE_HH #define COMPUTE_STATE_HH -double compute_state_update(double h, double unorm, double L, double old_state); +double state_update_dieterich(double h, double uol, double old_state); #endif diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index 15a0abf5340f8b088c19d6bc18e64817f02c3ecb..991213dc062ada0d49e1ba96952a8243dd968458 100644 --- a/src/one-body-sample.cc +++ b/src/one-body-sample.cc @@ -277,7 +277,7 @@ int main(int argc, char *argv[]) { // velocity // std::cout << std::log(L/unorm * h) << std::endl; - (*s4_new)[i] = compute_state_update(h, unorm, L, s4_old[i]); + (*s4_new)[i] = state_update_dieterich(h, unorm / L, s4_old[i]); } } if (parset.get<bool>("printProgress")) {