How to set limiter for CG compositional fields?

ASPECT
1

Hi folks,

After improving the convergence of Stokes subsystem, I got into another trouble. The problem is some compositional fields have negative values, which is not physically meaningful. I’m using P2CG for all compositional fields.
1.Why not using particles? Because the compositional fields have source terms, it wouldn’t be easy to update and project the required number of particles back and forth.
2.Why not using DG? DG will also suffer oscillations and if no limiter applied, the solution will gradually be contaminated. I’m using a 3D spherical shell, no limiters available yet.
3.Does P1CG help? yes, a little bit.
4.Is oscillation/undershoot/overshoot responsible for negative composition? yes, I double checked this by setting source term 0, negative values still appear.
5.Does smaller dt help? Yes, a lot. But I need dt as large as possible because this is a big simulation that needs a lot of wall time. A large dt saves wall time.

So I need to set up a limiter for compositional fields manually, but how? What I’m planning to do is setting a limiter after the solve but at the beginning of postprocess so that the postprocess and next step will have [0,1] values. I understand that the solution is protected (keyword const) for safety and consistency, it’s not easy to change the solution. I searched the forum and found prescribed velocity plugin is used to demonstrate how modifying solution is performed (Simulating multi-stage magma intrusion in ASPECT: How to dynamically update T & composition in specific spacetime domains?). But honestly, I don’t understand that at all.

So my question is, is there any ways to remove negative values (or decrease oscillations)? Do I have to use limiters? If so, where to start?

Deeply appreciate your time and help,
MM

2

An example of setting limiters on compositional fields:

namespace aspect
{
  namespace Postprocess
  {   
    template <int dim>
    std::pair<std::string,std::string>
    Laneuville<dim>::execute (TableHandler &statistics)
    {

      /* Concentration limiter: due to diffusion/oscillations, some negative concentrations must be corrected. */
      auto &solution = const_cast< LinearAlgebra::BlockVector &>(this->get_solution());
      const unsigned int n_comp = this->n_compositional_fields();
      const unsigned int peridotite_idx = this->introspection().compositional_index_for_name("peridotite");
      const auto &introspection = this->introspection();

      for(unsigned int nc = 0; nc < n_comp; ++nc)
      {
        const unsigned int idx = introspection.component_indices.compositional_fields[nc];
        auto &block = solution.block(idx);

        const auto &owned = block.locally_owned_elements();
        for (auto it = owned.begin(); it != owned.end(); ++it)
        {
          const auto i = *it;
          double value = block[I];

          if (nc == peridotite_idx)
          {
            std::cout<<"before = "<<value<<std::endl;
            value = std::max(0.0, std::min(1.0, value));
            std::cout<<"after = "<<value<<std::endl;
          }
          else
          {
            value = std::max(0.0, value);
          }

          block[i] = value;
        }
      }
      solution.compress(VectorOperation::insert);
      solution.update_ghost_values();
      /* Concentration limiter */

    /* Some other operation */
}}}

Hopefully this is useful.