SUPG Stabilization

#1

Hi,
I’m trying to use the SUPG advection stabilization in visco-plastic convection models. I’ve run into a problem where when I use supg stabilization instead of the default, the output heat flux is roughly twice the order of magnitude of the default (10^7 default vs 10^12 supg). I’ve tried running the supg stabilization on other material models with the same result except for using the simple material model. Any insight into this would be appreciated.
Thanks,
Erin Heilman

#2

Hi Erin,

I am not too surprised as the SUPG scheme has not been tested extensively, especially with all other features that might interact with it. SUPG is based on computing a residual of the temperature equation. My assumption is something that you are doing makes this residual very large.
What are your boundary conditions? Do you have a free surface?

#3

It might also be related to the compositional fields of the strain components. Here we use reaction terms to modify them. Maybe we don’t take these into account correctly.
Can you set up a simple example that you can run with and without SUPG and check if the compositional fields look similar?

#4

I have also noticed this issue where heat fluxes are computed many orders of magnitude too large when using SUPG. I have also been having problems with SUPG in 3D (both cartesian and spherical geometries, both with and without AMR). In every model that I have run, there are large (30%) temperature overshoots and undershoots and eventually the temperature solver will fail to converge, terminating the simulation. Increasing the ‘GMRES solver restart length’ will delay but not prevent crashing. I do have a relatively simple 3D test case that will produce an error but it takes about a day to run on 256 processors. Timo - would this be helpful?

#5

I wonder if this is the same problem that we had with the EV method that at Dirichlet boundaries we still compute the residual as if there were no boundary condition (making the residual extremely large). For EV we have code in source/entropy_viscosity.cc:343-412 that sets the artificial viscosity (=stabilization) to 0 at dirichlet boundaries. Maybe this should be moved outside the switch between EV and SUPG and just always apply?
Max do you mostly/only see these oscillations close to dirichlet boundaries, or also if you advect a temperature gradient in the interior of the domain or towards an insulating boundary?

#6

The over/under shoots seem to be confined to convective upwellings and downwellings. I have not run models with an insulating boundary. I could but it will take some time.

#7

I went back and checked everything running models using all the default parameters for the visco plastic material model. For the boundary conditions, I prescribe the outer spherical boundary at 273 K and the inner at 2973 K with no free surface. I looked and the residuals are small and the compositional fields are similar between the models.

#8

Rene,
I wonder if the unreasonably high heat flux values could be due to the inclusion of the stabilization term in the computation of heat flow? I remember some time ago that a decision was made to calculate heat fluxes using the product of the temperature gradient and the max(artificial viscosity,kThermal) where kThermal is the physical thermal conductivity. Is this being done also for SUPG, and does the SUPG stabilization parameter have the same physical units as thermal conductivity?

Max

#9

I made a one line change to entropy_viscosity.cc to disable stabilization near the boundary for SUPG in addition to EV and the heat fluxes are closer to what I would expect for the first couple of timesteps. This suggests that yes, the issue was related to the inclusion of the SUPG parameter when calculating heat fluxes. I am testing this branch now to see whether this also solves the failure to converge issue but it will take a day or so to run.

I think that if the non-physical contribution to thermal conductivity is going to be included in the heat transport, it should only be applied in the direction of the streamlines, so effectively (kI + supg_parameterv/norm(v)) * grad(T) where k is the thermal conductivity, I is the identity matrix, supg_parameter is the stabilization parameter and v/norm(v) is a unit vector along a streamline.

#10

The solution to Erin’s initial problem - unreasonably large reported values of heat flux - lies in the postprocessor for heat flux, which calculates heat flow as k_eff * grad(T). k_eff is defined as the maximum of the actual thermal conductivity and the artificial viscosity. For entropy viscosity stabilization, this kind of makes sense (though I would argue is not a useful output and should really come with a warning). The SUPG parameter has from what I can tell the MKS units of m^3 K/W, i.e. not a thermal conductivity, so it should not be included. The following change to source/postprocess/heat_flux_map.cc will disable this behavior for supg. I will submit a pull request for this issue in a few minutes and start a separate thread related to the convergence problems.

diff --git a/source/postprocess/heat_flux_map.cc b/source/postprocess/heat_flux_map.cc
index aabd84587..594b55f63 100644
--- a/source/postprocess/heat_flux_map.cc
+++ b/source/postprocess/heat_flux_map.cc
@@ -192,9 +192,14 @@ namespace aspect
                   const double material_prefactor = density_c_P + latent_heat_LHS;

                   const double artificial_viscosity_cell = static_cast<double>(artificial_viscosity(cell->active_cell_index()));
-                  const double diffusion_constant = std::max(out.thermal_conductivities[q],
-                                                             artificial_viscosity_cell);
-
+                 // The SUPG parameter tau does not have the physical dimensions of a diffusivity and as such should not be included in heat flux calculations
+                 const double diffusion_constant = (simulator_access.get_parameters().advection_stabilization_method == Parameters<dim>::AdvectionStabilizationMethod::supg) ?
+                   out.thermal_conductivities[q]
+                   :
+                   std::max(out.thermal_conductivities[q],
+                            artificial_viscosity_cell)
+                   ;
+
                   for (unsigned int i = 0; i<dofs_per_cell; ++i)
#11

Hi Erin,
your original problem should have been fixed by Max in the following pull request which was merged this morning: https://github.com/geodynamics/aspect/pull/3330

Could you check if your issue is solved in the current ASPECT development version?

Thanks for reporting the issue!
Rene

#12

Hi Rene,
I just tried it out and I’m getting the same values as the entropy viscosity stabilizer now. Thanks everyone!