[Meta] Improve accuracy of ASPECT solutions

[bangerth]

We've had a number of discussions over the years on how to improve the accuracy of the solutions ASPECT computes. Some of these have been implemented, for example @ian-r-rose's implementation of the consistent boundary flux method for the stress in #1572, or @gassmoeller's implementation for the heat flux in #2675.

Other ideas are proposed in #2711 and #2712. Let us collect what else we can come up here or in separate issues linked here.

[maxrudolph]

I had some discussions with @tjhei, @Shangxin-Liu, and Scott King at AGU about how to address these issues. To summarize, Scott and Shangxin have found through their efforts to reproduce the benchmarks in Zhong et al. 2008 that they can only reproduce the CitcomS benchmark values for temperature and heat flux if the stabilization parameters for the entropy viscosity scheme are set to zero. However, the benchmarks in Zhong et al. 2008 are run at low to moderate Rayleigh number. At higher Ra, turning off stabilization will lead to numerical instabilities. I have basically five ideas about how to improve the accuracy of the solution to the advection diffusion equation.

1. Attempt to find a more suitable choice of the stabilization parameters that would make the existing entropy viscosity implementation viable.
2. Explore the use of the DG discretization for temperature. My understanding is that the DG discretization does not require stabilization. Correct? We tried this already and encountered failures to converge as well as inconsistency of boundary values. Someone with deeper knowledge of DG would have to help with this.
3. Re-implement in ASPECT SUPG the way that it is implemented in CitcomS. I have read the implementation in CitcomS and would be happy to write out a detailed note on how it is implemented but I will need help implementing it in aspect/deal.ii. This may be more difficult than it appears at first glance because the elements used in ASPECT and CitcomS are different.
4. Advect temperature on particles following Gerya and Yuen (2007, 2009), also now done in StagYY apparently. Gerya and Yuen add an additional subgrid scale temperature relaxation scheme which is physically motivated but potentially unsatisfactory to some as it is somewhat ad-hoc.
5. Implement a different, yet-to-be-identified scheme for the advection-diffusion equation. In choosing a scheme, we should consider not just the convergence properties but also the accuracy at the levels of resolution that are accessible with current computing resources.

I have an additional suggestion, which is that at this time, there should be a warning in the ASPECT documentation about potential inaccuracies for high Ra calculations. I think that we would all agree that the temperature solution in 3D spherical geometry models (and possibly cartesian models as well) at any attainable resolution is non-physical and not subtly so. I would hate to see people invest large amounts of human or CPU time only to find out later that there are issues, and that these issues were known but not fully documented.

[tjhei]

To summarize, Scott and Shangxin have found through their efforts to reproduce the benchmarks in Zhong et al. 2008 that they can only reproduce the CitcomS benchmark values for temperature and heat flux if the stabilization parameters for the entropy viscosity scheme are set to zero

I talked to Grant last week and he said that at some point the results without stabilization where more accurate, so that is why the runs were done without. I don't think we know if some of the fixes that already went into ASPECT recently make a difference here. This needs to be explored of course. A good start would be to compare without/with EV/with SUPG for this benchmark.

• Re-implement in ASPECT SUPG the way that it is implemented in CitcomS. I have read the implementation in CitcomS and would be happy to write out a detailed note on how it is implemented but I will need help implementing it in aspect/deal.ii

Oh, we can handle the ASPECT side. Rene and I already have a preliminary implementation and we can finish that after we agree on the correct mathematical formulation (and comparing them to your notes). It would be very helpful to see what Citcom does. Should we schedule a video meeting to discuss? Are there code snippets worth looking at? Can you write down the equation being solved?

4. Advect temperature on particles

Might be interesting to try, but I don't think that is the solution we are looking for.

5. Implement a different, yet-to-be-identified scheme

Nonlinear SUPG/shock capturing comes to mind. This would be a neat student project I could tackle. Not something we can do quickly, though.

[maxrudolph]

I will write a note on exactly what is done in CitcomS but it will not happen until after the 27th as I’m on vacation with no laptop :). If you want to look at it in the meantime the code is pretty straightforward. It is in lib/Advection_diffusion.c if memory serves.

On Dec 17, 2018, at 10:43 PM, Timo Heister notifications@github.com wrote:

To summarize, Scott and Shangxin have found through their efforts to reproduce the benchmarks in Zhong et al. 2008 that they can only reproduce the CitcomS benchmark values for temperature and heat flux if the stabilization parameters for the entropy viscosity scheme are set to zero

I talked to Grant last week and he said that at some point the results without stabilization where more accurate, so that is why the runs were done without. I don't think we know if some of the fixes that already went into ASPECT recently make a difference here. This needs to be explored of course. A good start would be to compare without/with EV/with SUPG for this benchmark.

Re-implement in ASPECT SUPG the way that it is implemented in CitcomS. I have read the implementation in CitcomS and would be happy to write out a detailed note on how it is implemented but I will need help implementing it in aspect/deal.ii Oh, we can handle the ASPECT side. Rene and I already have a preliminary implementation and we can finish that after we agree on the correct mathematical formulation (and comparing them to your notes). It would be very helpful to see what Citcom does. Should we schedule a video meeting to discuss? Are there code snippets worth looking at? Can you write down the equation being solved?

Advect temperature on particles Might be interesting to try, but I don't think that is the solution we are looking for.

Implement a different, yet-to-be-identified scheme Nonlinear SUPG/shock capturing comes to mind. This would be a neat student project I could tackle. Not something we can do quickly, though.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

[gassmoeller]

I will not have time to look into this until the new year (technically I am on vacation), but I wanted to chime in, because I have something related to 1. that might be worth exploring: Currently the entropy viscosity scales with the residual of the advection equation, i.e. the less the solution changes over time, the less stabilization is needed. However, in the computation of the residual the entropy viscosity term itself is not included, meaning whenever we apply entropy viscosity we create a residual that leads to entropy viscosity for the next timestep. It would introduce another nonlinearity if we made the entropy viscosity dependent on the current solution, but we could store the entropy viscosity of past timesteps, and include them in the computation for the next timestep (essentially improving the explicit residual that we compute). This might already solve the weird situation that even for steady-state benchmarks we have a noticeable influence of the stabilization term. Does this make sense, or is there a reason we need this stabilization even for steady-state solutions?

[tjhei]

the code is pretty straightforward. It is in lib/Advection_diffusion.c if memory serves.

I browsed the code, but I don't understand the meaning of all the fields involved. It seems the parameter design for the SUPG parameter is

    adiff = (unorm>0.000001)?( (uxse*xse+ueta*eta+ufai*fai)/(2.0*unorm) ):0.0;

where unorm is the velocity norm and the uxse...fai is something like the size of the advection in each coordinate direction? There is a check where advection is compared to diffusion (twodiff), which I don't understand.

The second piece is the residual

Eres[j] -=
        PG->vpt[GNVINDEX(j,i)] * dOmega->vpt[i]
              * ((dT[i] + v1[i]*tx1[i] + v2[i]*tx2[i] + v3[i]*tx3[i])*rho*cp
                 - heating )
              + diff * dOmega->vpt[i] * E->heating_latent[m][el]
              * (GNx->vpt[GNVXINDEX(0,j,i)]*tx1[i]*rtf[3][i] +
                 GNx->vpt[GNVXINDEX(1,j,i)]*tx2[i]*sint[i] +
                 GNx->vpt[GNVXINDEX(2,j,i)]*tx3[i] );

which looks like (dT/dt + u.\nabla T)*rho*Cp + diffusion*nabla^2 T. I don't understand the details or how dT is computed, though.

Also interesting:

/* This function filters the temperature field. The temperature above   */
/* Tmax0(==1.0) and Tmin0(==0.0) is removed, while conserving the total */
/* energy. See Lenardic and Kaula, JGR, 1993.                           */

Finally:

/*   Functions which solve the heat transport equations using Petrov-Galerkin
     streamline-upwind methods. The process is basically as described in Alex
     Brooks PhD thesis (Caltech) which refers back to Hughes, Liu and Brooks.  */

[jaustermann]

Hi all - I just wanted to add that I'm interested in getting this fixed as well and would be happy to help running tests once that's useful. Max asked me at AGU about who else might have run into this issue. The only person I could think of was maybe @SiqiZhang, so looping him in here in case he's interested in this or has any insight.

[maxrudolph]

Also this may be relevant to @heronphi.

[sdkatvt]

the code is pretty straightforward. It is in lib/Advection_diffusion.c if memory serves.

I browsed the code, but I don't understand the meaning of all the fields involved. It seems the parameter design for the SUPG parameter is

    adiff = (unorm>0.000001)?( (uxse*xse+ueta*eta+ufai*fai)/(2.0*unorm) ):0.0;

where unorm is the velocity norm and the uxse...fai is something like the size of the advection in each coordinate direction? There is a check where advection is compared to diffusion (twodiff), which I don't understand.

You essentially have it right. There is a terse description of the method with no derivation of proof (referring to Hughes and Brooks, 1979) in my 1990 ConMan paper (3.2). I've never tried the drag and drop feature on github so I hope this works. If not I can post the conman paper somewhere or e-mail it. There is a lot of similarity between Citcom and ConMan at the FE level.

1990KingEtAl.pdf

The second piece is the residual

Eres[j] -=
      PG->vpt[GNVINDEX(j,i)] * dOmega->vpt[i]
              * ((dT[i] + v1[i]*tx1[i] + v2[i]*tx2[i] + v3[i]*tx3[i])*rho*cp
                 - heating )
              + diff * dOmega->vpt[i] * E->heating_latent[m][el]
              * (GNx->vpt[GNVXINDEX(0,j,i)]*tx1[i]*rtf[3][i] +
                 GNx->vpt[GNVXINDEX(1,j,i)]*tx2[i]*sint[i] +
                 GNx->vpt[GNVXINDEX(2,j,i)]*tx3[i] );

which looks like (dT/dt + u.\nabla T)*rho*Cp + diffusion*nabla^2 T. I don't understand the details or how dT is computed, though.

This is part of a predictor corrector method, and dT is the time derivative of the temperature field. Again a terse description is in the ConMan paper (section 3.2.1). The key point is that the PG Shape functions are only used on the advection terms, for the diffusion terms, the standard shape functions are used. If you look at equation 29 in the ConMan paper that's Eres[i].

We use a lumped mass matrix so there is no matrix solve and everything is explicit. There is an implicit version of all this too but I don't think I have a useful writeup of it.

Also interesting:

/* This function filters the temperature field. The temperature above   */
/* Tmax0(==1.0) and Tmin0(==0.0) is removed, while conserving the total */
/* energy. See Lenardic and Kaula, JGR, 1993.                           */

Yes, this came about in the days when we were using fields rather than tracers for composition so basically 0 or 1, essentially a shock. We realized we could do layered mantle convection by copying the advection-diffusion routine a second time with a second field. Myself and others looked at a lot of "shock capturing methods" to try and figure out how to deal with the problem rather than try to use something like particles and I found many of them didn't work as well as we might have hoped. Adrian designed a filter so that the over/undershoot mass was calculated and that mass was added back into the gradient region. It was a brute force method that we (I) pretty much only used for composition advection. It's not too far off from the idea of level sets (Arthur Raefsky and I almost got level sets worked out back in 1988, recognizing that the 0.5 contour was really accurately located and you really just wanted a way not to have to deal with the gradient). It should not be used for temperature, or at least I would not use it for temperature.

Finally:

/*   Functions which solve the heat transport equations using Petrov-Galerkin
     streamline-upwind methods. The process is basically as described in Alex
     Brooks PhD thesis (Caltech) which refers back to Hughes, Liu and Brooks.  */

Yes, I think I have copies of the Hughes, Liu and Brooks papers still somewhere. Maybe only hardcopy as pdf did not exist then. I also have a copy of Brooks' thesis. That is where the heavy lifting was gone. Hughes went on to recognize that SUPG was basically a special case of Galerkin Least Squares. Whether that still resonates in the FE world I don't know. It's about at that point that I decided I needed a thesis and could not keep playing around in finite element land, much as I would have liked to.

[sdkatvt]

Grant should be back next week. I can see with him exactly what versions he has tried. I know that with the latest heat flux and the artificial diffusion turned off he gets excellent agreement with Zhong et al 2008. We were not sure if this is just because it's a steady state solution or because even C1 (Ra=1.0e5) is such a low Rayleigh number that you don't need stabilization. I asked him to try getting the latest stabilization edits and try that. I also asked him to try discontinuous Galerkin. I don't know if he got to that or not. If someone has a working example of a discontinuous Galerkin problem for the temperature, we could adapt it and try.

[tjhei]

Thanks for your comments!

in my 1990 ConMan paper (3.2)

Great. I will take a look to see if this helps me in understanding.

We use a lumped mass matrix so there is no matrix solve and everything is explicit.

Just to be sure: this is also what is done in CitcomS still?

[maxrudolph]

Thanks for your comments!

in my 1990 ConMan paper (3.2)

Great. I will take a look to see if this helps me in understanding.

We use a lumped mass matrix so there is no matrix solve and everything is explicit.

Just to be sure: this is also what is done in CitcomS still?

It is explicit (no matrix solve) but the expression starting on line 634 does not look like it uses nodal quadrature. (dT+u.gradT)*(w+p) gets evaluated at the gauss quadrature points (E->N.vpt) from what I can tell. @sdkatvt can correct me if I'm wrong on this.

edit: I think that nodal quadrature must be implied by the omission of what is labeled M* in equation 28 of Scott's paper.

[sdkatvt]

Yeah. CitcomS also use a lumped mass matrix for the energy equation. It’s interesting how similar they are. Also interesting that we’re basically working with 1970’s and 80’s algorithms.

Scott

Sent from my iPhone

On Jan 9, 2019, at 5:41 PM, Timo Heister notifications@github.com wrote:

Thanks for your comments!

in my 1990 ConMan paper (3.2)

Great. I will take a look to see if this helps me in understanding.

We use a lumped mass matrix so there is no matrix solve and everything is explicit.

Just to be sure: this is also what is done in CitcomS still?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

[sdkatvt]

Indeed Eres[i] is computed at quadrature points but the code within "static void pg_solver" suggests that a lumped mass matrix is used and if you struggle around and find the right places where the variable are actually computed, indeed it is mass lumping.
The coding is pretty obscure but Josh and I spent a lot of this fall going through this section pretty carefully because we were trying to take melt calculated at integration points and project that to nodes. This is one area where CitcomS suffers somewhat from having many coders and the absence of a clear lead and/or clearly defined nomenclature standard within the code.

I’m not sure that the issue of a lumped mass matrix matters from the POV of SUPG elements, as I have a version of ConMan that uses an implicit solver with SUPG elements. In terms of actually trying to figure out what the routine means I can imagine for someone who didn’t “grow up with it” it is fairly challenging to parse what is what. LOL.

On Jan 9, 2019, at 6:34 PM, Max Rudolph notifications@github.com wrote:

Thanks for your comments!

in my 1990 ConMan paper (3.2)

Great. I will take a look to see if this helps me in understanding.

We use a lumped mass matrix so there is no matrix solve and everything is explicit.

Just to be sure: this is also what is done in CitcomS still?

It is explicit (no matrix solve) but the expression starting on line 634 does not look like it uses nodal quadrature. (dT+u.gradT)*(w+p) gets evaluated at the gauss quadrature points (E->n.vpt) from what I can tell. @sdkatvt https://github.com/sdkatvt can correct me if I'm wrong on this.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/geodynamics/aspect/issues/2713#issuecomment-452910090, or mute the thread https://github.com/notifications/unsubscribe-auth/AYbmOHjsKe6PdVN5EeMLzcHrs1HcV2btks5vBnyZgaJpZM4Y2ZWg.