Enable iterative viscosity dampening

[naliboff]

This [WIP] PR adds an option to average (dampen) the viscosity (or other properties) between nonlinear iterations, which has been reported to help improve convergence behavior.

At present, the PR is only partially complete as the method for storing (compositional field) and updating (reaction terms) the viscosity between nonlinear iterations does not work with standard compositional fields.

Here is a summary of what I think the issue is, borrowing from a comment by @jdannberg in #4923.

Ignoring advection, the reaction term update is C_new = C_old + R. Here, R is the reaction term and C_old is the old solution, which is not updated with nonlinear iterations.

The equation I would like to implement is
C_new = C_old + R = C_old + (-C_old + new_viscosity), where C_old is the viscosity is the previous nonlinear iteration and new_viscosity is the viscosity calculated in the most recent nonlinear iteration.

So, the first issue is that C_old is not updated during the nonlinear iterations. Can this be circumvented in some capacity by implementing a new custom field method?

The second issue is that the C_new values here are jumping between the same two values on subsequent nonlinear iterations, which is an issue @jdannberg also pointed out in the same comment.. However, I'm not sure why that is happening here as R is not equal to C_old

Hopefully there is path forward here.

[jdannberg]

I had a brief discussion with @anne-glerum and @bobmyhill about this, so just to document this here: I would use the PrescribedFieldOutputs here instead of reaction outputs. Make your viscosity field a prescribed field instead of solving an equation for it in every time step (that is cheaper anyway). You can then compute your update from the viscosity in the previous iteration stored in in.composition[i][this->introspection().compositional_index_for_name("viscosity_field")] and the new viscosity in out.viscosities[i].

Instead of solving, you would then just copy your prescribed field outputs into your compositional field (no reaction terms required). That happens in the interpolate_material_output_into_advection_field function. I just checked again and in.composition in this case is based on the values in the solution vector, and the solution vector for your prescribed field would be updated in this function as well (directly after you have computed your update). So that would solve your issues with the updates.

For an example, see this file:

aspect/tests/multiple_named_additional_outputs.cc

Lines 51 to 67 in 735ef46

	template <int dim>
	void
	PrescribedFieldMaterial<dim>::
	evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
	MaterialModel::MaterialModelOutputs<dim> &out) const
	{
	Simple<dim>::evaluate(in, out);
	// set up variable to interpolate prescribed field outputs onto compositional fields
	PrescribedFieldOutputs<dim> *prescribed_field_out = out.template get_additional_output<PrescribedFieldOutputs<dim>>();
	if (prescribed_field_out != NULL)
	for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
	{
	const double y = in.position[i](1);
	prescribed_field_out->prescribed_field_outputs[i][1] = std::exp(-y*y/2.0);
	}

[naliboff]

@jdannberg - fantastic, thank you! The approach you outlined should do the trick.