Compositional reaction using the particle method

[lhy11009]

Previously, I have been using the out.reaction_terms rate to reset the composition field in the ViscoPlastic Composition model. An example would be to include: for (unsigned int i=0; i < in.n_evaluation_points(); ++i) { ReactionRateOutputs<dim> *reaction_rate_out = out.template get_additional_output<ReactionRateOutputs<dim> >(); reaction_mor_compositions(i, out.reaction_terms, reaction_rate_out, in); } In the Evaluate function in the visco_plastic.cc file. But recently, I have had an issue that this only works with the field method but would be deprecated for the particle method.

So the implementation I want to have is to reset the composition at the trailing edge of the plate in 3D so that the subduction could continue without being restricted by the initial length of the plate. I have also made the choice of using the particle method to have a sharper boundary in the slab with the same mesh refinement. Here is a comparison of the different results: Upper: case with the field method, t = 14.0 Ma in model run time. Lower: a case with the particle method and t = 13.5 Ma in model run time. The left figures show the isovolume of the lower crust composition. Middle figures show another angle. The right figures show the viscosity. By looking at the volume of the lower-plate composition, the composition is reset at the trailing edge and the plates start from x = 600 km. On the other hand, the particle method fails to reset the composition and the plates start from x = 1500 km from the advance of the plate.

I would like to learn if there is a similar interface to apply composition reaction for the particle method as well? @gassmoeller @jdannberg, do you have time to look at this? I remembered the last time I adapted the out.reaction_terms term, it was from one of the plugins you included in the cookbook. Tell me what you think.

[lhy11009]

Okay I think it goes into this function in the advection.cc file: AdvectionSystem<dim>::execute, and here is where it's being parsed in: const double reaction_term = ((advection_field_is_temperature) ? 0.0 : scratch.material_model_outputs.reaction_terms[q][advection_field.compositional_variable]); Based on the name of it, it has to be related to an advection_field, does it mean it only works with a field method? I suspect the particle is not related to an advection field.

[gassmoeller]

Hi Haoyuan, Yes the reaction_terms are specifically used to compute the right-hand side of the compositional field advection equation. Therefore, they do not feed back into particle properties at all.

Your intend sounds like you want a particle property plugin that uses the reaction terms (or some other quantity you control) to update a particle property. Take a look at source/particle/property/elastic_stress.cc:124 and how it uses a material model to update the elastic stresses stored on the particles. This should give you an idea for how to do something similar for your purpose.

[lhy11009]

Hi @gassmoeller，Thanks for the information. I'll check that out first and keep this issue open for further conversation.