search

geodynamics / aspect

A parallel, extensible finite element code to simulate convection in both 2D and 3D models.
https://aspect.geodynamics.org/
Other
224 stars 235 forks source link

Models with constant composition do not converge #282

Closed gassmoeller closed 8 years ago

gassmoeller commented 9 years ago

The composition solver produces problems for constant compositional fields, because we use the right hand side of the equation for the tolerance computation, which is approximately zero in many of these cases. For now I always use slightly perturbed fields, but for some cases this might not be possible (e.g. if you track the sum of something), and it is nowhere documented that compositional fields should have at least some perturbation, so new users fall into this pit as well.

Any ideas for a better tolerance computation, or a criterion to skip the solution right away if the new solution would be equal to the old?

bangerth commented 9 years ago

Could you concoct a simple input file that shows the problem and attach it here?

gassmoeller commented 9 years ago

I have uploaded two parameter files here: https://gist.github.com/gassmoeller/a13190166616cd86c08a They describe the same model, but in one setup the density contrast is created by the temperature, in the other by the composition. There is no boundary condition providing additional force to drive the flow.

The file with a constant composition (where the temperature drives the flow) crashes with an up-to-date master in timestep 1 in debug mode, because the composition solver does not converge. The file with constant temperature (composition driving the flow) completes 27 timesteps before also crashing in the composition solver. This is around the time, when the velocities slowed down a lot and the fields are more or less constant. In the constant composition parameter file you see 3 commented lines, each of them separately changes the behaviour. Significantly decreasing the CFL number, increasing the Composition solver tolerance or adding a very small perturbation to the compositional field lets the model behave similar to the constant temperature case. The interesting thing is that it only happens in the composition solver, a constant temperature does not seem to create a problem.

ian-r-rose commented 9 years ago

I recall running into something similar to this with Timo a while back, and the fix for it was a bit of a hack, and may need revisiting. I'm not sure if this accounts for what you are currently seeing, but it could cut down on debug time.

Computing the entropy viscosity for the advection equation involves a division by the global variation in the advected field. If that is zero, there will clearly be problems. The current fix is to simply skip this step for a field with variation below some tolerance, but the fix/tolerance may not be chosen in the most sensible way. The relevant function for this is Simulator::compute_viscosity().

Anyways, this may or may not help. Best, Ian

On Tue, Feb 24, 2015 at 7:22 AM, Rene Gassmöller notifications@github.com wrote:

I have uploaded two parameter files here: https://gist.github.com/gassmoeller/a13190166616cd86c08a They describe the same model, but in one setup the density contrast is created by the temperature, in the other by the composition. There is no boundary condition providing additional force to drive the flow.

The file with a constant composition (where the temperature drives the flow) crashes with an up-to-date master in timestep 1 in debug mode, because the composition solver does not converge. The file with constant temperature (composition driving the flow) completes 27 timesteps before also crashing in the composition solver. This is around the time, when the velocities slowed down a lot and the fields are more or less constant. In the constant composition parameter file you see 3 commented lines, each of them separately changes the behaviour. Significantly decreasing the CFL number, increasing the Composition solver tolerance or adding a very small perturbation to the compositional field lets the model behave similar to the constant temperature case. The interesting thing is that it only happens in the composition solver, a constant temperature does not seem to create a problem.

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/aspect/issues/282#issuecomment-75776214.

gassmoeller commented 9 years ago

Hi Ian, thanks for the hint. I went through the assembly of the advection blocks for these parameter files line by line (debugging the values in the time step that does not converge), and I do think we handle the global variation right. In both parameter files timestep 1 uses the constant viscosity (assembly.cc:662 gets triggered), and the rest of the assembly is very similar. In fact the only difference I can see at the moment is that the temperature is multiplied by density and specific heat and the composition is not. I played around a bit, and interestingly the models behave more similar, when I do the multiplication also for the composition. They still do not converge at timestep 27 / 32 respectively, but the constant composition model now converges in timestep 1.

Therefore I wondered, if the constant viscosity value of 5e-3 * cell_diameter in assembly.cc:665 is only a reasonable value for the temperature, but not for the composition. The entries of the right-hand-side in the composition block are generally smaller by a factor of around 1e9 (rho_cp_T) compared to the temperature. A value of 5e-9 * cell_diameter also leads to convergence for the constant_composition model in timestep 1. This however, would not yet explain the non-convergence in the later timesteps.

bangerth commented 8 years ago

Can someone check whether this is still a problem?

gassmoeller commented 8 years ago

I think we can close this issue. The initial problem was solved by #321. The tests I posted above still crash after 26 timesteps, but this is likely related to #293, which we can address separately.