[aspect-devel] backward advection
Timo Heister
timo.heister at gmail.com
Mon Jan 18 08:54:58 PST 2016
The problem you are running into is that a negative delta_t will make
the linear systems negative definite. This breaks linear solvers like
CG. It might be easier to work with a positive delta_t than trying to
fix all places that will be messed up because of this.
On Mon, Jan 18, 2016 at 10:52 AM, Ian Rose <ian.rose at berkeley.edu> wrote:
> Hi Jacky,
>
> Since you are turning off diffusion, I wonder whether it would be enough to
> just reverse the direction of gravity. The timestepping would then not have
> to change at all, you would just interpret the times as negative rather than
> positive.
>
> Ian
>
> On Mon, Jan 18, 2016 at 8:30 AM, Austermann, Jacqueline
> <jaustermann at fas.harvard.edu> wrote:
>>
>> Hi everyone,
>>
>> I was wondering if anyone had thought about / implemented backward
>> advection in ASPECT. The way I would think to implement this is to let time
>> run in the negative direction (while setting diffusion to zero). I’ve
>> started looking through places where to implement this and was naively
>> trying out to just set time_step to - time_step in core.cc given some
>> boolean switch was set to true. I declared and parsed this switch parameter
>> in parameters.cc and adjusted termination_criteria/end_time.cc in which the
>> current time is queried. See the git diff below. I would expect that there
>> are other places that this might run into problems (namely where the time or
>> the timestep are queried), and I was hoping to fix these things one at a
>> time when they come up. However, I didn’t quite get there, because the
>> advection solver fails. With the changes below I tried running a cookbook
>> backward in time and it failed with the error message at the end of this
>> email. The advection solver did not converge with a residual of 1.70158e+18
>> (when I tried running the convection box cookbook backward the residual was
>> ~ 1e-5). I was going through the tutorial on how timestepping is done is
>> ASPECT and would think that it can handle negative time, but something is
>> clearly going wrong (maybe the solver can handle it and I’m just doing
>> something wrong).
>>
>> tl;dr I’m trying to implemented backward advection in ASPECT and am
>> wondering whether the advection solver can handle negative time or if there
>> is a better way to approach this than what is shown in the git diff.
>>
>> I’m happy to just keep thinking about this on my own but any advice or
>> insight would be appreciated.
>>
>> Best,
>> Jacky
>>
>>
>>
>> [jaustermann at rclogin06 simulator]$ git diff
>> diff --git a/include/aspect/parameters.h b/include/aspect/parameters.h
>> index 49ea611..8f3327c 100644
>> --- a/include/aspect/parameters.h
>> +++ b/include/aspect/parameters.h
>> @@ -161,6 +161,7 @@ namespace aspect
>> */
>> typename NonlinearSolver::Kind nonlinear_solver;
>>
>> + bool back_in_time;
>> double nonlinear_tolerance;
>> bool resume_computation;
>> double start_time;
>> diff --git a/source/simulator/assembly.cc b/source/simulator/assembly.cc
>> index d2c3dc1..0179ed1 100644
>> --- a/source/simulator/assembly.cc
>> +++ b/source/simulator/assembly.cc
>> @@ -724,7 +724,7 @@ namespace aspect
>> return max_viscosity;
>> else
>> {
>> - Assert (old_time_step > 0, ExcInternalError());
>> + Assert (old_time_step != 0, ExcInternalError());
>>
>> double entropy_viscosity;
>> if (parameters.stabilization_alpha == 2)
>> diff --git a/source/simulator/core.cc b/source/simulator/core.cc
>> index 0ea06a8..eea66b6 100644
>> --- a/source/simulator/core.cc
>> +++ b/source/simulator/core.cc
>> @@ -2071,6 +2071,10 @@ namespace aspect
>> parameters.maximum_time_step);
>> time_step =
>> termination_manager.check_for_last_time_step(time_step);
>>
>> + // switch time to go backward in time if switched on
>> + if (parameters.back_in_time == true)
>> + time_step = -time_step;
>> +
>> if (parameters.convert_to_years == true)
>> statistics.add_value("Time step size (years)", time_step /
>> year_in_seconds);
>> else
>> diff --git a/source/simulator/parameters.cc
>> b/source/simulator/parameters.cc
>> index 13cb69f..b26ac42 100644
>> --- a/source/simulator/parameters.cc
>> +++ b/source/simulator/parameters.cc
>> @@ -44,6 +44,9 @@ namespace aspect
>> Parameters<dim>::
>> declare_parameters (ParameterHandler &prm)
>> {
>> + prm.declare_entry ("Run backward in time", "false",
>> + Patterns::Bool (),
>> + "Option to let the calculation to go backward in
>> time.");
>> prm.declare_entry ("Dimension", "2",
>> Patterns::Integer (2,4),
>> "The number of space dimensions you want to run
>> this program in. "
>> @@ -708,6 +711,8 @@ namespace aspect
>> parse_parameters (ParameterHandler &prm,
>> const MPI_Comm mpi_communicator)
>> {
>> + back_in_time = prm.get_bool ("Run backward in time");
>> +
>> // first, make sure that the ParameterHandler parser agrees
>> // with the code in main() about the meaning of the "Dimension"
>> // parameter
>> diff --git a/source/termination_criteria/end_time.cc
>> b/source/termination_criteria/end_time.cc
>> index 6ca50f2..f182ac1 100644
>> --- a/source/termination_criteria/end_time.cc
>> +++ b/source/termination_criteria/end_time.cc
>> @@ -30,16 +30,16 @@ namespace aspect
>> bool
>> EndTime<dim>::execute()
>> {
>> - return (this->get_time() > end_time);
>> + return (std::abs(this->get_time()) > end_time);
>> }
>>
>> template <int dim>
>> double EndTime<dim>::check_for_last_time_step (const double
>> time_step) const
>> {
>> - if ((this->get_time()<end_time)
>> + if ((std::abs(this->get_time())<end_time)
>> &&
>> - (this->get_time()+time_step > end_time))
>> - return end_time - this->get_time();
>> + (std::abs(this->get_time())+time_step > end_time))
>> + return end_time - std::abs(this->get_time());
>> else
>> return time_step;
>> }
>>
>>
>>
>>
>>
>>
>> [jaustermann at rclogin12 run]$ ../aspect-build/aspect shell_simple_2d.prm
>>
>>
>> -----------------------------------------------------------------------------
>> -- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion.
>> -- . version 1.4.0-pre
>> -- . running in DEBUG mode
>> -- . running with 1 MPI process
>> -- . using Trilinos
>>
>> -----------------------------------------------------------------------------
>>
>> Number of active cells: 3,072 (on 6 levels)
>> Number of degrees of freedom: 40,836 (25,090+3,201+12,545)
>>
>> *** Timestep 0: t=0 years
>> Solving temperature system... 0 iterations.
>> Rebuilding Stokes preconditioner...
>> Solving Stokes system... 30+7 iterations.
>>
>> Number of active cells: 5,298 (on 7 levels)
>> Number of degrees of freedom: 70,924 (43,594+5,533+21,797)
>>
>> *** Timestep 0: t=0 years
>> Solving temperature system... 0 iterations.
>> Rebuilding Stokes preconditioner...
>> Solving Stokes system... 30+8 iterations.
>>
>> Number of active cells: 7,938 (on 7 levels)
>> Number of degrees of freedom: 105,212 (64,686+8,183+32,343)
>>
>> *** Timestep 0: t=0 years
>> Solving temperature system... 0 iterations.
>> Rebuilding Stokes preconditioner...
>> Solving Stokes system... 30+7 iterations.
>>
>> Number of active cells: 10,962 (on 7 levels)
>> Number of degrees of freedom: 144,929 (89,118+11,252+44,559)
>>
>> *** Timestep 0: t=0 years
>> Solving temperature system... 0 iterations.
>> Rebuilding Stokes preconditioner...
>> Solving Stokes system... 30+7 iterations.
>>
>> Number of active cells: 17,538 (on 8 levels)
>> Number of degrees of freedom: 232,798 (143,162+18,055+71,581)
>>
>> *** Timestep 0: t=0 years
>> Solving temperature system... 0 iterations.
>> Rebuilding Stokes preconditioner...
>> Solving Stokes system... 30+8 iterations.
>>
>> Postprocessing:
>> Writing graphical output: output/solution-00000
>> RMS, max velocity: 0.0664 m/year, 0.133 m/year
>> Temperature min/avg/max: 973 K, 2486 K, 4273 K
>> Heat fluxes through boundary parts: -2.845e+04 W, 5.635e+04 W, 1830
>> W, 1830 W
>> Writing depth average output/depth_average.gnuplot
>>
>> *** Timestep 1: t=-83550.1 years
>> Solving temperature system...
>>
>> ----------------------------------------------------
>> Exception on processing:
>>
>> --------------------------------------------------------
>> An error occurred in line <446> of file
>> </n/home09/jaustermann/aspect/aspect/source/simulator/solver.cc> in function
>> double aspect::Simulator<dim>::solve_advection(const
>> aspect::Simulator<dim>::AdvectionField &) [with dim = 2]
>> The violated condition was:
>> false
>> The name and call sequence of the exception was:
>> ExcMessage (std::string("The iterative advection solver " "did not
>> converge. It reported the following error:\n\n") + exc.what())
>> Additional Information:
>> The iterative advection solver did not converge. It reported the following
>> error:
>>
>>
>> --------------------------------------------------------
>> An error occurred in line <895> of file
>> </n/home09/jaustermann/bin_new/dealii-install/include/deal.II/lac/solver_gmres.h>
>> in function
>> void dealii::SolverGMRES<VECTOR>::solve(const MATRIX &, VECTOR &,
>> const VECTOR &, const PRECONDITIONER &) [with MATRIX =
>> dealii::TrilinosWrappers::SparseMatrix, PRECONDITIONER =
>> dealii::TrilinosWrappers::PreconditionILU, VECTOR =
>> dealii::TrilinosWrappers::MPI::Vector]
>> The violated condition was:
>> iteration_state == SolverControl::success
>> The name and call sequence of the exception was:
>> SolverControl::NoConvergence (accumulated_iterations, last_res)
>> Additional Information:
>> Iterative method reported convergence failure in step 1000. The residual
>> in the last step was 1.70158e+18.
>>
>> This error message can indicate that you have simply not allowed a
>> sufficiently large number of iterations for your iterative solver to
>> converge. This often happens when you increase the size of your problem. In
>> such cases, the last residual will likely still be very small, and you can
>> make the error go away by increasing the allowed number of iterations when
>> setting up the SolverControl object that determines the maximal number of
>> iterations you allow.
>>
>> The other situation where this error may occur is when your matrix is not
>> invertible (e.g., your matrix has a null-space), or if you try to apply the
>> wrong solver to a matrix (e.g., using CG for a matrix that is not symmetric
>> or not positive definite). In these cases, the residual in the last
>> iteration is likely going to be large.
>> --------------------------------------------------------
>>
>> --------------------------------------------------------
>>
>> Aborting!
>> ----------------------------------------------------
>>
>> _______________________________________________
>> Aspect-devel mailing list
>> Aspect-devel at geodynamics.org
>> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
>
>
>
> _______________________________________________
> Aspect-devel mailing list
> Aspect-devel at geodynamics.org
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
More information about the Aspect-devel
mailing list