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
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.
@optimux - Can you explain a bit further why particles would not be practical in this case? Is it simply the number of new particle properties you would need to track?
That point aside, the new particle interpolation algorithms (linear least squares, quadratic least squares) allow a high degree of accuracy with a relatively low number of particles per cell (relative to cell average) and also have the option to use limiters. Unless there really is a major barrier to using particles, I would strongly suggest this route.
Hello Mingming,
While the DG limiter has not been expressly tested for spherical geometry and so you receive a message that it is not available in spherical. However, we have been using this in my group (2D and 3D chunk geometry) with no problems. For some projects we started with the DG limiter but then shifted to particles because it was faster, however the dynamics was the same. This confirmed that the DG limiter is working fine in spherical coordinates.
I recommend commenting out the assert that prevents you from using the limiter in spherical coordinates and then doing some tests to demonstrate to yourself that it works.
Magali
Hi @jbnaliboff , thanks for your quick response. I have 5 compositional fields, one is peridotite which will experience melting, so phase change is related to this field. The rest are heat producing elements that will undergo radioactive decay. In my plugin, the amount of change of these fields during dt will be calculated and returned. So the number of particles will decrease accordingly. What l’m concerned is how this decrease is performed in ASPECT and do I need to do this manually (maybe not, I’m on a bus now, cannot check the manual, and i have never used particles in ASPECT before)?
Cheers,
MM
Dear @mibillen,
Sorry to interrupt and add a quick comment here.
However, we have been using this in my group (2D and 3D chunk geometry) with no problems.
I have indeed been using the DG limiter in my 3D spherical models, but so far I have not noticed a clear effect. Perhaps I need to do more tests on simpler models to better understand its behavior.
For some projects we started with the DG limiter but then shifted to particles because it was faster, however the dynamics was the same.
Also, regarding the point about switching from the DG limiter to particles for some projects, I am curious. For my previous viscoelastic benchmark, I couldn’t get the particle method to reach the same accuracy and speed as the field-based approach. Do you have any insights or suggestions on this?
Thank you very much for your time and help!
Best,
Ninghui
Hi Magali,
Thanks for your reply! I remembered that I ignored the warning regarding the availability of DG limiter in spherical shell, by just setting options Global composition maximum and Global composition minimum accordingly. I did not see any effect then, so I was suspecting by default the limiter was not applied or in my case the limiter was not called at all. It was a long time ago so I’m not 100% sure about this, thanks @tiannh7 by the way.
I’m also interested in the accuracy and speed that the particle method may deliver. Could you please share more details? Really appreciate!
Cheers,
MM
Hi @jbnaliboff ,
I skimmed briefly the particle section in the manual, and found that the source term can be handled by particle addition algorithm and particle removal algorithm. However, the energy associated with the source term (e.g., latent heat during phase change) cannot be returned to energy equation, esp. when the amount of phase change is iteratively solved (that is, for each iteration, the energy from phase change should be returned). So I’m not aware of any examples that have implemented similar algorithm. Any ideas?
Best,
MM
Hi Magali,
Looks like the major part has been eaten by the browser. Could you please show us more details?
Cheers,
MM
Here is the response that I tried to send before. My response to Nighui are in blue - magali
Sorry to interrupt and add a quick comment here.
However, we have been using this in my group (2D and 3D chunk geometry) with no problems.
I have indeed been using the DG limiter in my 3D spherical models, but so far I have not noticed a clear effect. Perhaps I need to do more tests on simpler models to better understand its behavior.
Mingming was seeing overshoots and undershoots without the limiter, if you turn is on these should go away and you see sharp boundaries maintained between compositions. So, that’s the effect he should see. You could then do some test to track total mass or something like that to show the that the compositional field isn’t growing or shrinking with time.
For some projects we started with the DG limiter but then shifted to particles because it was faster, however the dynamics was the same.
Also, regarding the point about switching from the DG limiter to particles for some projects, I am curious. For my previous viscoelastic benchmark, I couldn’t get the particle method to reach the same accuracy and speed as the field-based approach. Do you have any insights or suggestions on this?
I’m not the expert on this, but rather my former PhD student did these tests. My understanding is that you don’t actually need a lot of particles to get the same accuracy. The speed of the field approach versus particles also depends on which solver your are using. My understanding is that the GMG solver is faster with particles, because assembling the equations to solve for multiple fields ends up being more expensive than just updating particles. This may only be noticeable if you have a lot of fields (like when you are track stress components for visco-elastic rheology models.
Thank you very much for your time and help!
1 Like
Hi Magali,
Thanks for your update. I’d like to do some tests on DG limiter in next few days and will talk to you later.
Cheers,
MM
Dear @mibillen,
Thanks a lot for your helpful explanation and experience sharing. I’ll try some tests with the DG limiter on simple models and track the total mass as you suggested. Your insights on particles versus field methods, especially regarding solver efficiency, clarified many points for me.
Appreciate your time and support!
Best,
Ninghui
Hello @mibillen ,
Some updates on DG+limiter.
1.DG+limter in spherical shell. The prm looks like the following. This is DG1for all compositional fields.
subsection Discretization
set Composition polynomial degree = 1
set Temperature polynomial degree = 2
set Use discontinuous composition discretization = true
subsection Stabilization parameters
set Global composition maximum = 1.0, 1.0, 1.0, 1.0, 1.0
set Global composition minimum = 0.0, 0.0, 0.0, 0.0, 0.0
set Use limiter for discontinuous composition solution = true
set Stabilization method = entropy viscosity # SUPG
set beta = 0.117 #0.078 for 2D, 0.117 for 3D
set cR = 0.33
end
end
However, the option set Use limiter for discontinuous composition solution = true will cause the following error and the code will never run:
Exception 'ExcMessage("The limiter for the discontinuous temperature and composi
tion solutions has not been tested in non-Cartesian geometries and currently
requires the use of a Cartesian geometry model.")’ on rank xxx on processing:
You mentioned you can bypass the error by commenting out the assert, do you mean commenting the source code accordingly and recompile?
2.DG without limiter. Then I’m running models with set Use limiter for discontinuous composition solution = false, but the global maximum and minimum are still present. I was hoping these extrema can set up the constraints, but the heat-producing elements could still be negative, probably due to their low concentrations (the main composition field peridotite is in [0, 1] and never becomes negative using DG). So the conclusion is that the global max and min do not work.
3.DG0 works so far. Only DG0 preserves all concentrations positive. In my first post, I used DG1 and never tried DG0, so that conclusion is now completed.
Thanks for joining the discussion @jbnaliboff @tiannh7 , really appreciate your time and help!
Cheers,
MM