New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Time stepping frame work #1803

Merged
merged 28 commits into from Jun 30, 2017

Conversation

Projects
None yet
3 participants
@wenqing
Member

wenqing commented May 11, 2017

This PR presents a frame work for time stepping and a first application of adaptive time stepping, EvolutionaryPIDcontroller, an automatic time step control.

TODO:
1 make the example compared with synchronized fixing time stepping.
2 change all project files.

{
OGS_FATAL(
"Unknown time stepping type: `%s'. "
"The available types are \n\tSingleStep, \n\tFixedTimeStepping"

This comment has been minimized.

@wenqing

wenqing May 11, 2017

Member

FixedTimeStepping --> EvolutionaryPIDcontroller

@endJunction

This comment has been minimized.

Member

endJunction commented May 11, 2017

xmlstarlet it was. They do have move transformation.

@wenqing wenqing requested review from endJunction and chleh and removed request for endJunction May 23, 2017

@wenqing wenqing force-pushed the wenqing:dt0 branch from 31d939d to daa11ca May 23, 2017

@chleh

Many comments, but no fundamental issues.
Most importantly I think:

  • MathLib::VecNormType::INVALID should be avoided as a marker for a special case.
  • Maybe we should talk (in a greater round) about how to treat norms close to zero in general.
MathLib::VecNormType norm_type)
{
if (norm_type == MathLib::VecNormType::INVALID)
return 0.;

This comment has been minimized.

@chleh

chleh May 31, 2017

Contributor

Here there has to be OGS_FATAL() IMHO.

This comment has been minimized.

@wenqing

wenqing Jun 23, 2017

Member

Added

 OGS_FATAL("An invalid norm type has been passed");
return 0.;
const double norm2_x = MathLib::LinAlg::norm(x, norm_type);
assert(norm2_x > std::numeric_limits<double>::epsilon());

This comment has been minimized.

@chleh

chleh May 31, 2017

Contributor

This assertion might also be interesting in release builds.
On the other hand a solution x might be plain 0.

This comment has been minimized.

@wenqing

wenqing Jun 23, 2017

Member

removed this assertion.

if (norm_type == MathLib::VecNormType::INVALID)
return 0.;
const double norm2_x = MathLib::LinAlg::norm(x, norm_type);

This comment has been minimized.

@chleh

chleh May 31, 2017

Contributor

Misnomer? It's not always the 2-norm, is it?

This comment has been minimized.

@wenqing

wenqing Jun 23, 2017

Member

Changed to norm_x

// dx = x - x_old --> x - dx --> dx
MathLib::LinAlg::axpy(dx, -1.0, x);
return MathLib::LinAlg::norm(dx, norm_type) / norm2_x;

This comment has been minimized.

@chleh

chleh May 31, 2017

Contributor

Maybe add: // TODO implement per-component norms.

This comment has been minimized.

@wenqing

wenqing Jun 23, 2017

Member

The specific norm type of norm2 is removed.

return MathLib::LinAlg::norm(dx, norm_type) / norm2_x;
}
double BackwardEuler::getRelativeError(GlobalVector const& x,

This comment has been minimized.

@chleh

chleh May 31, 2017

Contributor

Maybe rename computeRelativeError to getRelativeError and remove getRelativeError from the subclasses.

Maybe the method name should be something like getRelativeChangeFromPreviousTimestep.

This comment has been minimized.

@wenqing

wenqing Jun 23, 2017

Member

Changed to getRelativeChangeFromPreviousTimestep.

local_assembler.assembleWithCoupledTerm(t, local_x, _local_M_data,
_local_K_data, _local_b_data,
local_coupling_term);
if (local_coupled_xs0.size() == 0 || local_coupled_xs.size() == 0)

This comment has been minimized.

@chleh

chleh May 31, 2017

Contributor

Minor: size() == 0empty()

This comment has been minimized.

@wenqing

wenqing Jun 23, 2017

Member

Replaced.

coupling_term.dt, coupling_term.coupled_processes,
std::move(local_coupled_xs0), std::move(local_coupled_xs));
_jacobian_assembler->assembleWithJacobianAndCouping(

This comment has been minimized.

@chleh

chleh May 31, 2017

Contributor

Btw.: There is a typo in the method name. assembleWithJacobianAndCoupling

This comment has been minimized.

@wenqing

wenqing Jun 23, 2017

Member

Corrected.

local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
_local_M_data, _local_K_data, _local_b_data, _local_Jac_data,
local_coupling_term);
if (local_coupled_xs0.size() == 0 || local_coupled_xs.size() == 0)

This comment has been minimized.

@chleh

chleh May 31, 2017

Contributor

Do the if and else branches have to be swapped, here?

This comment has been minimized.

@wenqing

wenqing Jun 23, 2017

Member

Gosh, it is wrong.

Corrected.

" <tol> 1.e-3 </tol>"
" <norm_type> NORM2 </norm_type>"
"</time_stepping>";
auto const PIDSteper = createTestTimeStepper(xml);

This comment has been minimized.

@chleh

chleh May 31, 2017

Contributor

Minor: typo.

This comment has been minimized.

@wenqing

wenqing Jun 23, 2017

Member

Corrected.

ts = PIDSteper->getTimeStep();
h_new = ts.dt();
ASSERT_EQ(2u, ts.steps());
ASSERT_NEAR(t_previous, ts.previous(), 1.e-10);

This comment has been minimized.

@chleh

chleh May 31, 2017

Contributor

Is no stricter tolerance than 1e-10 possible here?

This comment has been minimized.

@wenqing

wenqing Jun 23, 2017

Member

Reduced the tolerance.

@endJunction

This comment has been minimized.

Member

endJunction commented Jun 14, 2017

@wenqing How is the progress?

@wenqing

This comment has been minimized.

Member

wenqing commented Jun 15, 2017

@endJunction Just started revising.

@wenqing wenqing force-pushed the wenqing:dt0 branch from daa11ca to 9c82903 Jun 16, 2017

@wenqing wenqing force-pushed the wenqing:dt0 branch 7 times, most recently from 57dd81c to b84be8f Jun 16, 2017

@wenqing wenqing force-pushed the wenqing:dt0 branch from b84be8f to 22d91c2 Jun 26, 2017

*
* \param x The solution at the current timestep.
* \param norm_type The type of global vector norm.
*/

This comment has been minimized.

@endJunction

endJunction Jun 26, 2017

Member

Please add special behaviour (if some norms are < eps) to the documentation.

This comment has been minimized.

@wenqing

wenqing Jun 27, 2017

Member

Different algorithms have different method to calculate the error ( of the change ) of solutions. This special behaviour can be documented in the class of algorithm.

@wenqing wenqing force-pushed the wenqing:dt0 branch from 22d91c2 to 424e77c Jun 26, 2017

*
* @warning the value of x_old is changed to x - x_old after this computation.
*/
double computeRelativeChangeFromPreviousTimestep(

This comment has been minimized.

@endJunction

endJunction Jun 26, 2017

Member

Is it necessary to put this function into the interface; So far I have seen it is only used in the .cpp file, in the implementation.

This comment has been minimized.

@wenqing

wenqing Jun 27, 2017

Member

Since I want _dx, vector for the absolute solution increment, dx, to be created only once, I need this as a member of class. It can be move to cpp as a static function but vector dx has to be created and released in it.

This comment has been minimized.

@wenqing

wenqing Jun 27, 2017

Member

My initial implementation of it was in cpp.

_e_n_minus1 * _e_n_minus1 / (e_n * _e_n_minus2),
_kD) *
h_n;
else

This comment has been minimized.

@endJunction

endJunction Jun 26, 2017

Member

style: braces would improve readability

This comment has been minimized.

@wenqing

wenqing Jun 27, 2017

Member

Added.

}
h_new = limitStepSize(h_new, h_n);
const double checked_h_new = checkSpecificTimeReached(h_new);

This comment has been minimized.

@endJunction

endJunction Jun 26, 2017

Member

minor: one could reuse h_new...

This comment has been minimized.

@wenqing

wenqing Jun 27, 2017

Member

reused.

return h_new;
const double specific_time = _specific_times.back();
const double zero_threshlod = std::numeric_limits<double>::epsilon();

This comment has been minimized.

@endJunction

endJunction Jun 26, 2017

Member

typo: threshold.

This comment has been minimized.

@wenqing

wenqing Jun 27, 2017

Member

Corrected.

/// Create an EvolutionaryPIDcontroller time stepper from the given
/// configuration
std::unique_ptr<TimeStepAlgorithm> createEvolutionaryPIDcontroller(

This comment has been minimized.

@endJunction

endJunction Jun 26, 2017

Member

Would suggest to move this function definition into own file. Also create a separate header for the declaration.

This comment has been minimized.

@wenqing

wenqing Jun 27, 2017

Member

Created h/cpp files for this creator.

wenqing added some commits Mar 22, 2017

[dt] Added a member to classes of time discretization to get relative…
… solution error

and moved some member definitions to cpp file.
[dt] Introduced a member _dx to store the solution increment
as a buffer for computing the relative error.

@wenqing wenqing force-pushed the wenqing:dt0 branch from 96e8087 to ea57db1 Jun 29, 2017

[dt] Use the norm type of ConvergenceCriterion in the error
caculation for the time steppers.

@wenqing wenqing force-pushed the wenqing:dt0 branch from ea57db1 to 364bd9d Jun 29, 2017

@wenqing

This comment has been minimized.

Member

wenqing commented Jun 29, 2017

All tests are green now.
@chleh @endJunction Thanks for the comments.

@wenqing wenqing referenced this pull request Jun 30, 2017

Merged

Adaptive timestepping #60

@endJunction endJunction merged commit 80f8821 into ufz:master Jun 30, 2017

3 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/jenkins/pr-merge This commit looks good
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment