-
Notifications
You must be signed in to change notification settings - Fork 235
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
[Meta] Improve accuracy of ASPECT solutions #2713
Comments
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.
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. |
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.
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?
Might be interesting to try, but I don't think that is the solution we are looking for.
Nonlinear SUPG/shock capturing comes to mind. This would be a neat student project I could tackle. Not something we can do quickly, though. |
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 ***@***.***> 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.
|
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: |
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
where unorm is the velocity norm and the The second piece is the residual
which looks like Also interesting:
Finally:
|
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. |
Also this may be relevant to @heronphi. |
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.
This is part of a predictor corrector method, and dT is the time derivative of the temperature field. 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.
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.
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. |
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. |
Thanks for your comments!
Great. I will take a look to see if this helps me in understanding.
Just to be sure: this is also what is done in CitcomS still? |
It is explicit (no matrix solve) edit: I think that nodal quadrature must be implied by the omission of what is labeled M* in equation 28 of Scott's paper. |
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 ***@***.***> 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.
|
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 ***@***.***> 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 <#2713 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AYbmOHjsKe6PdVN5EeMLzcHrs1HcV2btks5vBnyZgaJpZM4Y2ZWg>.
|
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.
The text was updated successfully, but these errors were encountered: