Within AdvectionField structs there’s the handy function is_temperature() to tell whether the struct refers to a temperature field or not. Is there any current method to access the same information in a plugin’s evaluate() function and tell whether the temperature field is currently being solved or not?
Ryan: What kind of plugin are you considering? If it is, for example, an InitialTemperature plugin, then yes: the field being considered is the temperature. If it is InitialComposition, then it isn’t. So the context often helps. If you tell us which plugin system you are considering, we might be able to help 
Best
W.
Thank you for the response!
I was considering a MaterialModel plugin. Specifically, is it possible to tell during the preconditioner building phase in the advection system whether a temperature system is about to be solved or a composition system within evaluate() for a material model?
My main motivation was looking at assembly.cc (lines 892-893) where the material model is called for local_assemble_advection_system
material_model->evaluate(scratch.material_model_inputs,
scratch.material_model_outputs);
Here advection_field (that has the is_temperature function) is accessible inside local_assemble_advection_system() but not inside the invoked evaluate().
template <int dim>
void Simulator<dim>::
local_assemble_advection_system (const AdvectionField &advection_field,
const Vector<double> &viscosity_per_cell,
const typename DoFHandler<dim>::active_cell_iterator &cell,
internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
internal::Assembly::CopyData::AdvectionSystem<dim> &data)
{
Many thanks again for the response.
Ryan: No, I don’t think there is a way to determine inside the MaterialModel::evaluate() function where it was invoked from. The function takes certain inputs and is asked to create certain outputs; who that is for shouldn’t matter to the function.
Why do you need to know this? What would you do differently inside the function if you knew?
Best
W.
1 Like
That makes sense.
It follows up a bit on our discussion during the user meeting about my attempts to advect a field representing free fluid upwards at a constant rate. I had some instabilities show up when I tried to modify the advection equation by changing the reaction_term.
I had tested your suggestion of capping the time step duration to see if that would eliminate the instabilities, but the instabilities persisted independent of time step duration. However, the instabilities do not show up if the background temperature is homogeneous, so I was curious if turning off the reaction_term assignment in the material model when the temperature is being solved for would improve things. I also changed the adaptive mesh refinement parameters; no change there.
Anyhow, the whole point is a bit moot now because today I was able to obtain my desired behavior using the darcy field setup, but I was curious about the cause of the instabilities all the same.
Ryan,
Great to hear the Darcy Field setup worked for you!
Just for testing the effect of different terms, it might be easiest to just comment these out in the assembler itself (i.e. within this function: aspect/advection.cc at ef63de22597567192111921060aa2cce452847d0 · geodynamics/aspect · GitHub). This also shows that for the assembly of the temperature field, the reaction_term is always set to zero anyway (lines 135-140), so it does not affect what happens in the assembly of the temperature equation.
(Also note that the Darcy field uses a different assembler: aspect/advection.cc at ef63de22597567192111921060aa2cce452847d0 · geodynamics/aspect · GitHub)
Cheers,
Juliane
1 Like
Yes, very glad it worked. It was a modification of the other assembler that ended up producing the correct results (lines 561, 621).
Thank you for pointing out the zeroing of the reaction_term in the temperature equation. That addresses my original underlying motivation for trying to access the is_temperature() information from evaluate().
I’m not sure yet if it’s possible to identify the cause of the instabilities but if so it won’t be from the reaction_term in the temperature equation. Given that there’s a viable workaround, I’ll be pivoting in that direction.
Many thanks to you both for the explanations and tips!