Check that we do not use a discontinuous pressure when we compute pressure gradients

[bangerth]

Pressure gradients are not useful quantities, yet they are computed all over the place in the code:

rprojects/aspect> grep -r pressure_gradient source/
source/simulator/assembly.cc:          std::vector<Tensor<1,dim> > old_pressure_gradients;
source/simulator/assembly.cc:          std::vector<Tensor<1,dim> > old_old_pressure_gradients;
source/simulator/assembly.cc:          old_pressure_gradients (quadrature.size(), Utilities::signaling_nan<Tensor<1,dim> >()),
source/simulator/assembly.cc:          old_old_pressure_gradients (quadrature.size(), Utilities::signaling_nan<Tensor<1,dim> >()),
source/simulator/assembly.cc:          old_pressure_gradients (scratch.old_pressure_gradients),
source/simulator/assembly.cc:          old_old_pressure_gradients (scratch.old_old_pressure_gradients),
source/simulator/assembly.cc:            scratch.old_pressure_gradients);
source/simulator/assembly.cc:            scratch.old_old_pressure_gradients);
source/simulator/assembly.cc:            scratch.explicit_material_model_inputs.pressure_gradient[q] = (scratch.old_pressure_gradients[q] + scratch.old_old_pressure_gradients[q]) / 2;
source/simulator/assembly.cc:        material_model_inputs.pressure_gradient);
source/simulator/assembly.cc:        scratch.old_pressure_gradients);
source/simulator/assembly.cc:        scratch.old_old_pressure_gradients);
source/simulator/assembly.cc:          scratch.explicit_material_model_inputs.pressure_gradient[q] = (scratch.old_pressure_gradients[q] + scratch.old_old_pressure_gradients[q]) / 2;
source/simulator/lateral_averaging.cc:                                                                                           in.pressure_gradient);
source/simulator/helper_functions.cc:    std::vector<Tensor<1,dim> > pressure_gradients(n_q_points);
source/simulator/helper_functions.cc:                                                                                   pressure_gradients);
source/simulator/helper_functions.cc:                  in.pressure_gradient[q] = pressure_gradients[q];
source/simulator/nullspace.cc:              fe[introspection.extractors.pressure].get_function_gradients (relevant_dst, in.pressure_gradient);
source/simulator/nullspace.cc:              fe[introspection.extractors.pressure].get_function_gradients (relevant_dst, in.pressure_gradient);
source/mesh_refinement/thermal_energy_density.cc:                                                                                         in.pressure_gradient);
source/mesh_refinement/viscosity.cc:                                                                                         in.pressure_gradient);
source/mesh_refinement/density.cc:                                                                                         in.pressure_gradient);
source/heating_model/adiabatic_heating.cc:            heating_model_outputs.heating_source_terms[q] = (material_model_inputs.velocity[q] * material_model_inputs.pressure_gradient[q])
source/heating_model/latent_heat.cc:                                                          * (material_model_inputs.velocity[q] * material_model_inputs.pressure_gradient[q]);
source/material_model/interface.cc:      pressure_gradient(n_points, aspect::Utilities::signaling_nan<Tensor<1,dim> >()),
source/postprocess/visualization/dynamic_topography.cc:                .get_function_gradients (this->get_solution(), in.pressure_gradient);
source/postprocess/visualization/density.cc:                in.pressure_gradient[q][d] = duh[q][this->introspection().component_indices.pressure][d];
source/postprocess/visualization/thermal_expansivity.cc:                in.pressure_gradient[q][d] = duh[q][this->introspection().component_indices.pressure][d];
source/postprocess/visualization/shear_stress.cc:                in.pressure_gradient[q][d] = duh[q][this->introspection().component_indices.pressure][d];
source/postprocess/visualization/friction_heating.cc:                in.pressure_gradient[q][d] = duh[q][this->introspection().component_indices.pressure][d];
source/postprocess/visualization/viscosity.cc:                in.pressure_gradient[q][d] = duh[q][this->introspection().component_indices.pressure][d];
source/postprocess/visualization/stress.cc:                in.pressure_gradient[q][d] = duh[q][this->introspection().component_indices.pressure][d];
source/postprocess/visualization/material_properties.cc:                in.pressure_gradient[q][d] = duh[q][this->introspection().component_indices.pressure][d];
source/postprocess/visualization/heating.cc:                in.pressure_gradient[q][d] = duh[q][this->introspection().component_indices.pressure][d];
source/postprocess/visualization/specific_heat.cc:                in.pressure_gradient[q][d] = duh[q][this->introspection().component_indices.pressure][d];
source/postprocess/visualization/vertical_heat_flux.cc:                in.pressure_gradient[q][d] = duh[q][this->introspection().component_indices.pressure][d];
source/postprocess/mass_flux_statistics.cc:                    in.pressure_gradient);
source/postprocess/heat_flux_statistics.cc:                    in.pressure_gradient);
source/postprocess/dynamic_topography.cc:              .get_function_gradients (this->get_solution(), in.pressure_gradient);
source/postprocess/viscous_dissipation_statistics.cc:                                                                                         in.pressure_gradient);
source/postprocess/heating_statistics.cc:                                     in.pressure_gradient);
source/adiabatic_conditions/initial_profile.cc:            in.pressure_gradient[0] = g/(g.norm() != 0.0 ? g.norm() : 1.0)
source/adiabatic_conditions/initial_profile.cc:            in.pressure_gradient[0] = Tensor <1,dim> ();

I think we should look at which of these terms are in fact used and ensure that in all of the other places, we don't compute them at all.

[bangerth]

Specifically, MaterialModelInputs also has them. I don't think it should be in there.

[bangerth]

And for reference, @tjhei 's #710 already takes care of some of these locations.

[jdannberg]

We need the pressure gradients in some of the heating plugins, for example in the adiabatic heating. I think that's why we added them to the MaterialModelInputs.

[bangerth]

@jdannberg -- that's awkward. How are you going to deal with the case of someone using a discontinuous pressure element (the "conservative" option)? Can you state what term it is you need this in?

[jdannberg]

There is an option to use density * gravity instead of the pressure gradient (it is called Use simplified adiabatic heating), but at the moment we do not make sure that if the discontinuous pressure element is used one can only use the simplified adiabatic heating. So that is something to improve.

The term we need that in is the RHS term of the adiabatic heating: _velocity * pressure_gradient * thermal_expansioncoefficient * temperature

and the latent heat: _temperature * density * entropy_derivative_pressure * velocity * pressuregradient.

So we could use density * gravity everywhere, but the pressure gradient is more exact.

[bangerth]

I don't know if the numerical gradient of the pressure is actually more exact. But in either case, we should probably add an assertion about which element is used in all of the places where the gradient is used.