[aspect-devel] backward advection

Ian Rose ian.rose at berkeley.edu
Mon Jan 18 08:52:54 PST 2016


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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20160118/b723dfa3/attachment-0001.html>


More information about the Aspect-devel mailing list