SUPG Convergence Issues

maxrudolph commented 4 years ago

I have been encountering failures to converge when using SUPG in 3D. We actually haven't been successful in running a single 3D model to completion with SUPG. This includes both 3D spherical geometry models and 3D cartesian models, both with and without AMR. I have tried to come up with the simplest possible testcase, which I have attached. This crashes after about 40 timesteps due to failure to converge. It takes less than a minute to run on 32 cores.

solver_history.txt supg_testcase.prm.txt log.txt

I haven't tried using SUPG in 2D so I'm not sure if these issues are 3D-specific. In 2D it is generally cheap enough to run at a high enough refinement level to reduce the entropy viscosity to values below the thermal conductivity.

tjhei commented 4 years ago

Thanks for this test case. :+1: I have not used SUPG for any larger or more realistic computations, so I am not too surprised that you are encountering problems. I have some ideas on how to fix this. Will report back.

tjhei commented 4 years ago

This is the error I get in debug mode after 62 timesteps:


*** Timestep 62:  t=2.3924e-05 seconds
   Solving temperature system... 30 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 83+0 iterations.

   Postprocessing:
     RMS, max velocity:                  9.45e+04 m/s, 5.35e+05 m/s
     Temperature min/avg/max:            -1.013 K, 0.4936 K, 1.603 K
     Heat fluxes through boundary parts: 6906 W, -6906 W, -3992 W, 3992 W, -2883 W, 1645 W
     Writing graphical output:           output_supg_testcase/solution/solution-00062

--------------------------------------------------------
An error occurred in line <605> of file <../source/simulator/helper_functions.cc> in function
    double aspect::Simulator<3>::compute_time_step() const [dim = 3]
The violated condition was:
    rho * c_p > 0
Additional information:
    The product of density and c_P needs to be a non-negative quantity.

A couple of things about this:

if rho*cp is negative, I am not surprised that the temperature solver will fail.
We should rethink if this should also be detected in release mode @gassmoeller ?
Should the material model maybe make sure that the density stays positive?
Your setup is severely underresolved, which probably explains the situation. Maybe we need to tune the SUPG parameters (in 3d)?

maxrudolph commented 4 years ago

It shouldn't be possible for any of the density terms in the energy equation to become negative. This calculation uses Boussinesq approximation so density is constant and equal to the reference density in the energy equation. Is rho*Cp being used somewhere in the stabilization terms in such a way that it's inconsistent with the "Formulation" selected? The momentum equation could contain a negative density if temperature undershoots are severe.

tjhei commented 4 years ago

negative density if temperature undershoots are severe.

Yes, we have about 20% undershoots, so T=-0.2, which makes the density negative locally.

maxrudolph commented 4 years ago

Yes, I agree that the undershoots are severe. The magnitude of the density over/undershoots should decrease at higher resolution and this test case was meant to fail quickly - hence it's severely under-resolved as you say. However, for a boussinesq problem I think that the undershoots still should not matter. The density used in the energy equation is the reference density and the only place where density should contribute to the momentum equation is the buoyancy term, which should not care about whether density is positive or negative. Now, if density is being used in calculating stabilization parameters or scaling factors (e.g. scaling the pressure to improve condition number) I could see how this becomes a problem.

maxrudolph commented 4 years ago

Following up on this, I think that this code in helper_functions.cc should likely be modified for Boussinesq models (and maybe others too?) to use the reference values of k,rho, and Cp.

starting line 595 in helper_functions.cc:


              // calculate the corresponding conduction timestep, if applicable
              for (unsigned int q=0; q<n_q_points; ++q)
                {
                  const double k = out.thermal_conductivities[q];
                  const double rho = out.densities[q];
                  const double c_p = out.specific_heat[q];

                  Assert(rho * c_p > 0,
                         ExcMessage ("The product of density and c_P needs to be a "
                                     "non-negative quantity."));

                  const double thermal_diffusivity = k/(rho*c_p);

                  if (thermal_diffusivity > 0)
                    {
                      min_local_conduction_timestep = std::min(min_local_conduction_timestep,
                                                               parameters.CFL_number*pow(cell->minimum_vertex_distance(),2)```

tjhei commented 4 years ago

Following up on this, I think that this code in helper_functions.cc should likely be modified for Boussinesq models (and maybe others too?) to use the reference values of k,rho, and Cp.

yes, I agree, see #3356 I think the SUPG formulation itself also needs to be fixed. I will take a look there next.

tjhei commented 4 years ago

I think the SUPG formulation itself also needs to be fixed. I will take a look there next.

SUPG parameters are computed correctly.

maxrudolph commented 4 years ago

Thank you.

gassmoeller commented 4 years ago

That looks promising. I do not yet see a place where the SUPG residual would be computed with the real density for BA though. Before the residual computation the following code is executed in entropy_viscosity.cc:531

        if (parameters.formulation_temperature_equation
            == Parameters<dim>::Formulation::TemperatureEquation::reference_density_profile)
          {
            // Overwrite the density by the reference density coming from the
            // adiabatic conditions as required by the formulation
            for (unsigned int q=0; q<n_q_points; ++q)
              scratch.material_model_outputs.densities[q] = adiabatic_conditions->density(scratch.material_model_inputs.position[q]);
          }
        else if (parameters.formulation_temperature_equation
                 == Parameters<dim>::Formulation::TemperatureEquation::real_density)
          {
            // use real density
          }
        else
          AssertThrow(false, ExcNotImplemented());

which replaces the real density with the reference density for all the stabilization terms computations.

Edit: Timo fixed it before I could finish my post :-)

gassmoeller commented 4 years ago

In #3356 Timo posted:

well, it reduces the timestep size significantly so it is not crashing for me (at least in the first 100 steps).

But I do not understand this. If we use a reference density similar to the real density then over-undershootings and using the real density should cause timesteps that are too small, not timesteps that are too big. I think the model setup has a problem in that it uses an adiabatic surface temperature of 0.0 (default value, not prescribed in file), and a reference temperature for the simple material model of 0.5. This mean the reference density in the material model is 1.0 (at T=0.5), while the reference density that is used in the equations is actually 1.5 (at T=0.0), I am not sure this is intended?. But even considering that I can not see how this change should decrease the timestep instead of increasing it. Could you (or @maxrudolph) post density and adiabatic density paraview output for this model (by adding 'density, adiabat' to the 'List of output variables')?

maxrudolph commented 4 years ago

I can do that. I did not know that the adiabatic surface temperature would matter in this case since there is no adiabat under the Boussinesq approximation.

maxrudolph commented 4 years ago

Running with #3356 I am still seeing a convergence failure on the second timestep. I ran with the debug version.

tjhei commented 4 years ago

second timestep

Any other changes to the .prm you posted? How many cores? It is running for 100+ timesteps for me now.

maxrudolph commented 4 years ago

Oops - I accidentally left artificial viscosity smoothing turned on from a test. Re-running now.

maxrudolph commented 4 years ago

This does seem much better. Mine has made it to ~250 timesteps. I will try the higher-resolution testcase again.

Shangxin-Liu commented 4 years ago

Hi all @gassmoeller @tjhei @maxrudolph ,

I hope everyone is staying well at home during this special period. Following the last short discussion in the virtual meeting, we'd @GEuen like to also report the detail of the similar advection solver convergence failure with SUPG in our attempt to benchmark high Ra thermal convection. First of all, this issue indeed seems obscure because low Ra number (10^3 - 10^5) cases work very well now with SUPG but high Ra number cases (10^8) still lead to instability and convergence failure of the advection solver. We're able to come up with a failed high Ra case in 3D spherical shell with the public ASPECT version and would like to share this to help debug this issue. The crash happens after about 70 tilmestep. I ran it in the debug mode and it takes ~40 minutes to crash with 32 cores (1 node) on our cluster. The Release mode should be much faster. This test is ran with 16 layers. The negative temperature does not appear but the advection solver still fail to converge. We also plot the time evolution of the average temperature, RMS velocity, Nu_t, and Nu_b. We find that this instability is pretty obscure that only the Nu_t plot can reveal the over/undershoot problems before crash. I attach the prm file, error file, statistic plots here. I'm also running it in 32 layers, but higher resolution is expected to need more time to crash. I guess debugging SUPG in lower resolution is more time economic.

high_Ra_case2_CFL_1.0.prm.txt slurm-317472.out.txt log.txt solver_history.txt statistics.txt average_T_evolution.pdf average_V_RMS_evolution.pdf Nu_b.pdf Nu_t.pdf

Feel free to let me know if you need more material to help isolate this issue. @gassmoeller @tjhei

tjhei commented 4 years ago

We find that this instability is pretty obscure that only the Nu_t plot can reveal the over/undershoot problems before crash. I attach the prm file, error file, statistic plots here.

excellent, thank you for the test case. Can you share "venus_like.txt" please?

Shangxin-Liu commented 4 years ago

oops, I forgot this data file. here it is. @tjhei venus_like.txt

tjhei commented 4 years ago

I can reproduce the crash. Awesome, this is very helpful! :+1:

tjhei commented 4 years ago

I am a bit concerned that you lowered the Stokes tolerance to 1e-3. I am playing a bit with solver settings and I notice that it is quite difficult to get the solver to converge with the default residual. An inaccurate velocity can cause nonphysical behavior as well. It is worth figuring out why the Stokes solve is so difficult, I think.

Shangxin-Liu commented 4 years ago

I'm pretty sure that 1e-3 tolerance here is enough. The velocity is physically reasonable under this value at least for these kinds of cases from my previous experience. I've played with different values of Stokes linear solver tolerance. For the low Ra test cases, for example, 1e-3 and 1-7 lead to almost the same results both within 1% error. Yes, using the default solver with AMG, the Stokes solver takes a number of iterations to converge within the first 2-3 timesteps. However, after that, it behaves normally. I guess it's due to the sharp viscosity jump of the four layer reference viscosity (I set the depth range of the jump to very small) and the strong temperature dependent viscosity. Today I ran one more test using 1e-5 Stokes linear tolerance but set the GMRES solver length to 150 to accelerate the convergence. I also attach the prm file and the statistic plotting here. The advection solver still crash with failed convergence after ~60 timesteps.

high_Ra_case2_CFL_1.0.prm.txt slurm-320970.out.txt log.txt solver_history.txt statistics.txt average_T_evolution.pdf average_V_RMS_evolution.pdf Nu_b.pdf Nu_t.pdf

Let me know whether you can also reproduce this case. @tjhei We can talk more tomorrow.

Shangxin-Liu commented 4 years ago

I tested one more case to further turn off the depth-dependent viscosity based on my last one and this case also crashed after ~140 timesteps with a failed convergence of advection solver. The whole setting is still nondimensional with 1e-5 linear Stokes tolerance. Without radial jump, it only takes a number of Stokes iterations at tilmestep 0 (200 + 49) and after that the default AMG Stokes solver becomes easier (only coarse outer iterations). Because the viscosity is much smaller for this case (at constant nondimensional 1) than the above depth-dependent one, the nondimenisional velocity is much larger. The under/overshoot of temperature solution is even more obvious for this case, which indeed indicates that faster speed tends to result in instability more easily. Just before the time point of crash, negative temperature also appears. The instability can be seen in RMS velocity, Nu_b, and Nu_t (extreme oscillation). This case may be more helpful to debug SUPG because it is even simpler. I still attach everything here. @tjhei @gassmoeller @maxrudolph @bangerth

high_Ra_case2_CFL_1.0_no_depth_visc.prm.txt log.txt slurm-321747.out.txt solver_history.txt statistics.txt average_T_evolution.pdf average_V_RMS_evolution.pdf Nu_b.pdf Nu_t.pdf

maxrudolph commented 4 years ago

Interesting. It makes sense that the crash happens even more quickly with smaller viscosities (higher Ra, vrms). I think that my case without depth dependent viscosity may have omitted the lowest viscosity layer rather than reducing all of the viscosities to the minimum value.

On Wed, Mar 25, 2020 at 8:40 PM Shangxin Liu notifications@github.com wrote:

I tested one more case to further turn off the depth-dependent viscosity based on my last one and this case also crashed after ~140 timesteps with a failed convergence of advection solver. The whole setting is still nondimensional with 1e-5 linear Stokes tolerance. Without radial jump, it only takes a number of Stokes iterations at tilmestep 0 (200 + 49) and after that the default AMG Stokes solver becomes easier (only coarse outer iterations). Because the viscosity is much smaller for this case (at constant nondimensional 1) than the above depth-dependent one, the nondimenisional velocity is much larger. The under/overshoot of temperature solution is even more obvious for this case, which indeed indicates that faster speed tends to result in instability more easily. Just before the time point of crash, negative temperature also appears. The instability can be seen in RMS velocity, Nu_b, and Nu_t (extreme oscillation). This case may be more helpful to debug SUPG because it is even simpler. I still attach everything here. @tjhei https://github.com/tjhei @gassmoeller https://github.com/gassmoeller @maxrudolph https://github.com/maxrudolph @bangerth https://github.com/bangerth

high_Ra_case2_CFL_1.0_no_depth_visc.prm.txt https://github.com/geodynamics/aspect/files/4384790/high_Ra_case2_CFL_1.0_no_depth_visc.prm.txt log.txt https://github.com/geodynamics/aspect/files/4384791/log.txt slurm-321747.out.txt https://github.com/geodynamics/aspect/files/4384792/slurm-321747.out.txt solver_history.txt https://github.com/geodynamics/aspect/files/4384793/solver_history.txt statistics.txt https://github.com/geodynamics/aspect/files/4384794/statistics.txt average_T_evolution.pdf https://github.com/geodynamics/aspect/files/4384795/average_T_evolution.pdf average_V_RMS_evolution.pdf https://github.com/geodynamics/aspect/files/4384796/average_V_RMS_evolution.pdf Nu_b.pdf https://github.com/geodynamics/aspect/files/4384798/Nu_b.pdf Nu_t.pdf https://github.com/geodynamics/aspect/files/4384799/Nu_t.pdf

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/geodynamics/aspect/issues/3335#issuecomment-604210770, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6G6NXZKYZHEDK3Y6HX2WDRJLFCFANCNFSM4KH2KHCA .

tjhei commented 4 years ago

I have been working on this some more and noticed a few things:

It turns out that these examples generate negative densities at some point due to temperature oscillations (often earlier than the code actually crashing). This means it would fail early in debug mode if we were patient enough to run it that way. We should get back to validation #2025, #3358, #919. For now, I just modified the Simple model to always return a positive density. I will report back how far I can run with this.
We can also tune the preconditioner/solve for the temperature equation to be more robust and continue further. I need to experiment with this more before making a PR.

tjhei commented 4 years ago

when enforcing positive density, the last example is still running at

*** Timestep 429:  t=5.51588e-06 seconds

maxrudolph commented 4 years ago

Is this a Boussinesq model? If so I'm surprised that the density used would not be the reference density.

On Tue, Apr 7, 2020 at 12:56 PM Timo Heister notifications@github.com wrote:

when enforcing positive density, the last example is still running at

*** Timestep 429: t=5.51588e-06 seconds

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/geodynamics/aspect/issues/3335#issuecomment-610589328, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6G6NTZX6CWNAQX73ZS6CLRLOAOJANCNFSM4KH2KHCA .

tjhei commented 4 years ago

Is this a Boussinesq model? If so I'm surprised that the density used would not be the reference density.

I agree. This is using Boussinesq. Let me dig into this some more.

tjhei commented 4 years ago

Let me dig into this some more.

SUPG correctly uses the adiabatic density. Maybe negative densities are enough to cause nonsensical velocities from the buoyancy term.

maxrudolph commented 4 years ago

That makes sense.

On Tue, Apr 7, 2020 at 2:16 PM Timo Heister notifications@github.com wrote:

Let me dig into this some more.

SUPG correctly uses the adiabatic density. Maybe negative densities are enough to cause nonsensical velocities from the buoyancy term.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/geodynamics/aspect/issues/3335#issuecomment-610624811, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6G6NX7QZMZ4B5DS2FUJ53RLOJ2LANCNFSM4KH2KHCA .

Shangxin-Liu commented 4 years ago

I use "function" to set adiabatic conditions in my prm file, which prescribes the reference density in energy equation to 1.0, so this should not be a problem. The negative density may appear in the momentum equation due to the temperature overshoot from the energy equation, and then further cause the non-physical velocity solution, which goes into energy equation again. Anyway, the crash of the advection solver seems coming from both temperature over/undershoot and non-physical velocity of the energy equation. @tjhei @maxrudolph

tjhei commented 4 years ago

It turns out that this problem has nothing to do with negative densities. Instead, the iterative advection solver seems to struggle in certain situations with SUPG (maybe with large viscosity jumps). The PR #3453 works around the problem for me.

maxrudolph commented 4 years ago

@tjhei I tried to build #3453 this morning but had errors compiling stokes_matrix_free.cc. Do I need to update my deal.ii for this to work?

/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc: In instantiation of ‘std::pair<double, double> aspect::StokesMatrixFreeHandlerImplementation<dim, velocity_degree>::solve() [with int dim = 2; int velocity_degree = 2]’:
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:2323:3:   required from here
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:1552:9: error: ‘void dealii::PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::estimate_eigenvalues(const VectorType&) const [with MatrixType = aspect::MatrixFreeStokesOperators::ABlockOperator<2, 2, double>; VectorType = dealii::LinearAlgebra::distributed::Vector<double>; PreconditionerType = dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double> >]’ is private within this context
         mg_smoother_A[level].estimate_eigenvalues(temp_velocity);
         ^~~~~~~~~~~~~
In file included from /home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/generic_linear_algebra.h:23:0,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/global.h:29,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/stokes_matrix_free.h:25,
                 from /home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:22:
/home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/precondition.h:2239:1: note: declared private here
 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
 ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:1553:9: error: ‘void dealii::PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::estimate_eigenvalues(const VectorType&) const [with MatrixType = aspect::MatrixFreeStokesOperators::MassMatrixOperator<2, 1, double>; VectorType = dealii::LinearAlgebra::distributed::Vector<double>; PreconditionerType = dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double> >]’ is private within this context
         mg_smoother_Schur[level].estimate_eigenvalues(temp_pressure);
         ^~~~~~~~~~~~~~~~~
In file included from /home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/generic_linear_algebra.h:23:0,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/global.h:29,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/stokes_matrix_free.h:25,
                 from /home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:22:
/home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/precondition.h:2239:1: note: declared private here
 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
 ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[ 82%] Building CXX object CMakeFiles/aspect.dir/unit_tests/additional_outputs.cc.o
[ 82%] Building CXX object CMakeFiles/aspect.dir/unit_tests/first.cc.o
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc: In instantiation of ‘std::pair<double, double> aspect::StokesMatrixFreeHandlerImplementation<dim, velocity_degree>::solve() [with int dim = 2; int velocity_degree = 3]’:
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:2323:3:   required from here
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:1552:9: error: ‘void dealii::PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::estimate_eigenvalues(const VectorType&) const [with MatrixType = aspect::MatrixFreeStokesOperators::ABlockOperator<2, 3, double>; VectorType = dealii::LinearAlgebra::distributed::Vector<double>; PreconditionerType = dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double> >]’ is private within this context
         mg_smoother_A[level].estimate_eigenvalues(temp_velocity);
         ^~~~~~~~~~~~~
In file included from /home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/generic_linear_algebra.h:23:0,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/global.h:29,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/stokes_matrix_free.h:25,
                 from /home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:22:
/home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/precondition.h:2239:1: note: declared private here
 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
 ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:1553:9: error: ‘void dealii::PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::estimate_eigenvalues(const VectorType&) const [with MatrixType = aspect::MatrixFreeStokesOperators::MassMatrixOperator<2, 2, double>; VectorType = dealii::LinearAlgebra::distributed::Vector<double>; PreconditionerType = dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double> >]’ is private within this context
         mg_smoother_Schur[level].estimate_eigenvalues(temp_pressure);
         ^~~~~~~~~~~~~~~~~
In file included from /home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/generic_linear_algebra.h:23:0,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/global.h:29,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/stokes_matrix_free.h:25,
                 from /home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:22:
/home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/precondition.h:2239:1: note: declared private here
 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
 ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc: In instantiation of ‘std::pair<double, double> aspect::StokesMatrixFreeHandlerImplementation<dim, velocity_degree>::solve() [with int dim = 3; int velocity_degree = 2]’:
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:2323:3:   required from here
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:1552:9: error: ‘void dealii::PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::estimate_eigenvalues(const VectorType&) const [with MatrixType = aspect::MatrixFreeStokesOperators::ABlockOperator<3, 2, double>; VectorType = dealii::LinearAlgebra::distributed::Vector<double>; PreconditionerType = dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double> >]’ is private within this context
         mg_smoother_A[level].estimate_eigenvalues(temp_velocity);
         ^~~~~~~~~~~~~
In file included from /home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/generic_linear_algebra.h:23:0,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/global.h:29,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/stokes_matrix_free.h:25,
                 from /home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:22:
/home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/precondition.h:2239:1: note: declared private here
 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
 ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:1553:9: error: ‘void dealii::PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::estimate_eigenvalues(const VectorType&) const [with MatrixType = aspect::MatrixFreeStokesOperators::MassMatrixOperator<3, 1, double>; VectorType = dealii::LinearAlgebra::distributed::Vector<double>; PreconditionerType = dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double> >]’ is private within this context
         mg_smoother_Schur[level].estimate_eigenvalues(temp_pressure);
         ^~~~~~~~~~~~~~~~~
In file included from /home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/generic_linear_algebra.h:23:0,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/global.h:29,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/stokes_matrix_free.h:25,
                 from /home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:22:
/home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/precondition.h:2239:1: note: declared private here
 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
 ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[ 82%] Building CXX object CMakeFiles/aspect.dir/unit_tests/parse_map_to_double_array.cc.o
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc: In instantiation of ‘std::pair<double, double> aspect::StokesMatrixFreeHandlerImplementation<dim, velocity_degree>::solve() [with int dim = 3; int velocity_degree = 3]’:
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:2323:3:   required from here
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:1552:9: error: ‘void dealii::PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::estimate_eigenvalues(const VectorType&) const [with MatrixType = aspect::MatrixFreeStokesOperators::ABlockOperator<3, 3, double>; VectorType = dealii::LinearAlgebra::distributed::Vector<double>; PreconditionerType = dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double> >]’ is private within this context
         mg_smoother_A[level].estimate_eigenvalues(temp_velocity);
         ^~~~~~~~~~~~~
In file included from /home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/generic_linear_algebra.h:23:0,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/global.h:29,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/stokes_matrix_free.h:25,
                 from /home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:22:
/home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/precondition.h:2239:1: note: declared private here
 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
 ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:1553:9: error: ‘void dealii::PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::estimate_eigenvalues(const VectorType&) const [with MatrixType = aspect::MatrixFreeStokesOperators::MassMatrixOperator<3, 2, double>; VectorType = dealii::LinearAlgebra::distributed::Vector<double>; PreconditionerType = dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double> >]’ is private within this context
         mg_smoother_Schur[level].estimate_eigenvalues(temp_pressure);
         ^~~~~~~~~~~~~~~~~
In file included from /home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/generic_linear_algebra.h:23:0,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/global.h:29,
                 from /home/rudolph/sw/candi/aspect-tjhei/include/aspect/stokes_matrix_free.h:25,
                 from /home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:22:
/home/rudolph/sw/dealii-toolchain-openmpi-4.0.2-hdf5/deal.II-master/include/deal.II/lac/precondition.h:2239:1: note: declared private here
 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
 ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[ 82%] Building CXX object CMakeFiles/aspect.dir/unit_tests/particles.cc.o
[ 83%] Building CXX object CMakeFiles/aspect.dir/unit_tests/termination_criteria.cc.o
[ 83%] Building CXX object CMakeFiles/aspect.dir/unit_tests/utilities.cc.o
[ 83%] Building CXX object CMakeFiles/aspect.dir/contrib/world_builder/source/coordinate_systems/cartesian.cc.o
[ 84%] Building CXX object CMakeFiles/aspect.dir/contrib/world_builder/source/coordinate_systems/interface.cc.o
[ 84%] Building CXX object CMakeFiles/aspect.dir/contrib/world_builder/source/coordinate_systems/spherical.cc.o
[ 84%] Building CXX object CMakeFiles/aspect.dir/contrib/world_builder/source/features/continental_plate.cc.o
[ 84%] Building CXX object CMakeFiles/aspect.dir/contrib/world_builder/source/features/continental_plate_models/composition/interface.cc.o
CMakeFiles/aspect.dir/build.make:6686: recipe for target 'CMakeFiles/aspect.dir/source/simulator/stokes_matrix_free.cc.o' failed
make[2]: *** [CMakeFiles/aspect.dir/source/simulator/stokes_matrix_free.cc.o] Error 1
make[2]: *** Waiting for unfinished jobs....
CMakeFiles/Makefile2:67: recipe for target 'CMakeFiles/aspect.dir/all' failed
make[1]: *** [CMakeFiles/aspect.dir/all] Error 2
Makefile:129: recipe for target 'all' failed
make: *** [all] Error 2

tjhei commented 4 years ago

to update my deal.ii for this to work?

I assume you are using an older development version of deal.II? The error would go away if you a) update deal.II to the most recent master or b) comment out the call to estimate_eigenvalues in /home/rudolph/sw/candi/aspect-tjhei/source/simulator/stokes_matrix_free.cc:1552. It doesn't have any negative effect.

maxrudolph commented 4 years ago

Thanks. I think that my development deal.ii was a couple of months stale.

maxrudolph commented 4 years ago

Unfortunately, #3453 does not fix the failure to converge that I've encountered.

gassmoeller commented 4 years ago

Should we do some more tests to exclude iterative solver stability as the reason for the crashes? Some suggestions (open for opinions about what makes sense and what not):

Increasing the GMRES restart length (setting the "Solver parameters/Advection solver parameters/GMRES solver restart length to 200 instead of 50)? This was a problem for the Stokes GMRES solver in some cases, and at least Shangxins example above does not converge once the Advection iterations exceed 50 (not conclusive, but a hint?).
Using an AMG instead of ILU preconditioner? I know too little about the ILU preconditioner we use. I am not too concerned about using more expensive but more effective preconditioners, as the time spent for building the advection preconditioner and solving the advection system is currently below 1% overall runtime even if the Stokes system is simple (assembling the matrix is expensive, but solving is cheap).

tjhei commented 4 years ago

Should we do some more tests to exclude iterative solver stability as the reason for the crashes?

That makes sense, yes. If I have a prm that breaks with #3453 I can continue to investigate.

Increasing the GMRES restart length

Good idea, I haven't tried this.

Using an AMG instead of ILU preconditioner?

I thought about that before as well and I am happy to try if I have an example where things break.

tjhei commented 4 years ago

For the record: I tried the file from Shangxin from March 23 with

  subsection Advection solver parameters
    set GMRES solver restart length = 150
  end

and I still need #3453 in step 75, 76, 77, 78, 115, etc.. I stopped at 819 and here are the results (overlay is Shangxin's solution until crash): supg-flux-top

GEuen commented 4 years ago

Hi all. I’ve been trying to run a 3D model with depth dependent viscosity and #3453 (which appears to have fixed the original crashing problem for me). Somehow though, the model will start, step forward once as I would expect, then each subsequent time step is at the exact same time as the first. Many output values are quite high as well. I’ve messed with several parameters and can’t seem to find the issue. I’ve attached the files for a 2D version of the model that shows the same behavior, as well as some graphs showing different values with respect to both time and time step. Hopefully, I’ve just done something silly here.

statistics.txt venus_like_dim.txt dimensional.prm.txt GMTHeatFlux_2D.pdf GMTHeatFluxDensity_2D.pdf GMTTemp_2D.pdf GMTVrms_2D.pdf

gassmoeller commented 4 years ago

Hi Grant, I cant tell you the underlying reason, but here are some clues (you may have seen that already, but you did not explicitly state it above): Your model does a single time step, after which it develops enormously large velocities (>1e5 m/yr), and all subsequent timesteps are extremely short as shown in the statistics file (4e-3 yrs per timestep), this is why it looks like time does not progress. The viscosities in your model look somewhat large, but fine (average viscosity in your statistics file 1e24), but your integrated total mass is negative .(-3.20753364e+28), this is what I would investigate. Did you run your model in debug mode? It may crash with an assertion about negative densities. Its probably something about the settings in your material model. What happens if you use the simple material model directly without the detour through the depth dependent material model?

GEuen commented 4 years ago

Here are two cases more recent cases. The only difference between them is one uses only the simple material model, the other uses depth dependent as before. The simple version has a convergence error very quickly. Depth dependent looks like the previous example. I'll try them both in debug mode as well.

Simple version: error_file.txt log.txt original.prm.txt parameters.prm.txt solver_history.txt statistics.txt venus_like_dim.txt

Depth Dependent version: GMTHeatFlux.pdf GMTHeatFluxDensity.pdf GMTTemp.pdf GMTVrms.pdf log.txt original.prm.txt parameters.prm.txt statistics.txt venus_like_dim.txt

gassmoeller commented 4 years ago

In the simple version your viscosity is 1 Pas and your Rayleigh number is 1e27, this will break your model quite easily. The depth dependent model is the interesting case because it breaks at higher viscosities, but I could guess that it maybe uses the correct viscosity, but uses an incorrect reference viscosity (i.e. the 1 Pas from the simple model)? We use the reference viscosity and a length scale to scale velocities to be on the same order of magnitude as the pressures (what the Rayleigh number does for nondimensional models), so using a reference viscosity that is vastly different than the actual viscosities will break the solvers. Can you try increasing the viscosity for the simple material model, and then scaling the depth dependent factors in your depth dependent file accordingly to create the same final viscosities?

Also this discussion is not really related to SUPG at the moment, could you open a separate issue or a forum discussion?

geodynamics / aspect

SUPG Convergence Issues #3335