diff --git a/dune/dorie/interface/transport_simulation.cc b/dune/dorie/interface/transport_simulation.cc index 0008e8dbc79f1b85b658483a9516c862e5ff8e34..1f9b241e8d84c8e002c2f9a28b61481b5b8b842e 100644 --- a/dune/dorie/interface/transport_simulation.cc +++ b/dune/dorie/interface/transport_simulation.cc @@ -144,9 +144,12 @@ void TransportSimulation::step () auto t = controller->getTime(); suggest_timestep(controller->getDT()); - // for some reason dt!=this->dt - const auto _dt = this->dt; + double courant_number = 0.5; //FIXME: add keyword in inifile + if (not ts_param->implicit()) + this->suggest_timestep(courant_number*Dune::Dorie::cfl_condition(*gf_water_flux)); + // for some reason dt!=this->dt, any idea why? + // Check whether time intervals of grid functions are valid if (Dune::FloatCmp::gt(water_flux_interval.begin,t+this->dt)) DUNE_THROW(Dune::Exception," "); diff --git a/dune/dorie/interface/transport_simulation.hh b/dune/dorie/interface/transport_simulation.hh index 5e0d9ba26adf302353d77ef85f52367cb0640dcd..fb7c7671252340f681457d05d982ab2aca313215 100644 --- a/dune/dorie/interface/transport_simulation.hh +++ b/dune/dorie/interface/transport_simulation.hh @@ -298,11 +298,7 @@ public: void suggest_timestep(TimeField _dt) override { - auto cfl_dt = 0.3*Dune::Dorie::cfl_condition(*gf_water_flux); // FIXME - std::cout << "cfl_dt: " << cfl_dt << std::endl; - this->dt = std::min(_dt,controller->getDT()); - this->dt = std::min(dt,cfl_dt); - std::cout << "dt: " << cfl_dt << std::endl; + dt = std::min(_dt,controller->getDT()); } // TODO: Define data exchange with richards diff --git a/dune/dorie/solver/util_cfl_condition.hh b/dune/dorie/solver/util_cfl_condition.hh index 473444b53379667e1ed1b05c30f2265a58dfecce..429494037628d2529b250c6590cbf17c8baa76c1 100644 --- a/dune/dorie/solver/util_cfl_condition.hh +++ b/dune/dorie/solver/util_cfl_condition.hh @@ -55,7 +55,7 @@ namespace Dorie{ for (int i = 1; i < geo.corners(); ++i) for (int j = 0; j < i-1; ++j) h_T = std::max(h_T,(geo.corner(i)-geo.corner(j)).two_norm()); - std::cout << "h_T: " << h_T << std::endl; + // Check that h_T is not zero in debug mode, who knows. assert(Dune::FloatCmp::gt(h_T,0.)); @@ -69,15 +69,12 @@ namespace Dorie{ gf.evaluate(cell,x,water_flux); max_water_flux = std::max(max_water_flux,water_flux.two_norm()); } - std::cout << "max_water_flux: " << max_water_flux << std::endl; + // Calculate the cfl for this cell if (Dune::FloatCmp::gt(max_water_flux,0.)) suggested_dt = std::min(suggested_dt,h_T/max_water_flux); - - std::cout << "suggested_dt: " << suggested_dt << std::endl; } - std::cout << "FINAL suggested_dt: " << suggested_dt << std::endl; // Communicate the lowest cfl number to all procesors suggested_dt = gv.comm().min(suggested_dt); return suggested_dt;