[aspect-devel] backward advection
Austermann, Jacqueline
jaustermann at fas.harvard.edu
Mon Jan 18 08:30:00 PST 2016
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<http://core.cc> given some boolean switch was set to true. I declared and parsed this switch parameter in parameters.cc<http://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<http://assembly.cc> b/source/simulator/assembly.cc<http://assembly.cc>
index d2c3dc1..0179ed1 100644
--- a/source/simulator/assembly.cc<http://assembly.cc>
+++ b/source/simulator/assembly.cc<http://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<http://core.cc> b/source/simulator/core.cc<http://core.cc>
index 0ea06a8..eea66b6 100644
--- a/source/simulator/core.cc<http://core.cc>
+++ b/source/simulator/core.cc<http://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<http://parameters.cc> b/source/simulator/parameters.cc<http://parameters.cc>
index 13cb69f..b26ac42 100644
--- a/source/simulator/parameters.cc<http://parameters.cc>
+++ b/source/simulator/parameters.cc<http://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<http://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!
----------------------------------------------------
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20160118/a4721f3c/attachment-0001.html>
More information about the Aspect-devel
mailing list