Skip to content
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

HDG for advection-diffusion on the sphere #34

Draft
wants to merge 27 commits into
base: master
Choose a base branch
from
Draft

Conversation

amartinhuertas
Copy link
Collaborator

Super WIP. Opening PR in draft mode ... just for discussion purposes

@codecov-commenter
Copy link

codecov-commenter commented Jun 29, 2022

Codecov Report

Merging #34 (d96d72a) into master (288cf35) will not change coverage.
The diff coverage is 0.00%.

@@           Coverage Diff           @@
##           master     #34    +/-   ##
=======================================
  Coverage    0.00%   0.00%            
=======================================
  Files          18      19     +1     
  Lines        1246    1346   +100     
=======================================
- Misses       1246    1346   +100     
Impacted Files Coverage Δ
src/AdvectionHDG.jl 0.00% <0.00%> (ø)
src/mpi/CubedSphereDistributedDiscreteModels.jl 0.00% <0.00%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 288cf35...d96d72a. Read the comment docs.

@davelee2804
Copy link
Collaborator

Hey @amartinhuertas , I'm having some basic trouble with the HDG advection on the sphere. I am getting the error:

It is not possible to perform the operation "dot" on the given cell fields.

here:
https://github.com/gridapapps/GridapGeosciences.jl/blob/hdg_adv/src/AdvectionHDG.jl#L19
Currently I am passing in un directly from the function:
https://github.com/gridapapps/GridapGeosciences.jl/blob/hdg_adv/driver/sequential/SBRAdvectionHDG.jl#L19
However I have also tried projecting un onto a FESpace via
https://github.com/gridapapps/GridapGeosciences.jl/blob/hdg_adv/src/AdvectionHDG.jl#L57
https://github.com/gridapapps/GridapGeosciences.jl/blob/hdg_adv/src/AdvectionHDG.jl#L90
I'm wondering if the vector field on the manifold and the L2 reference element vector field (as I've defined it) are somehow incompatible? (Note that D = 2)...

@amartinhuertas
Copy link
Collaborator Author

Can you please separate the terms in the bilinear form one by one as in e.g., https://github.com/gridap/GridapHybrid.jl/blob/main/test/DarcyHDGTests.jl#L71
and let me know which is the one(s) that triggers the error? Thanks!

@amartinhuertas
Copy link
Collaborator Author

(Note that D = 2).

Did you try with D=3?

@davelee2804
Copy link
Collaborator

davelee2804 commented Jun 30, 2022

...I tried D=3, but it failed even earlier (I can't remember where exactly sorry). Thanks for the suggestion, the failure is at:
_op2 = (un⋅n). It also manifests here: (∇(q)⋅un)

@amartinhuertas
Copy link
Collaborator Author

can you attach a file with the whole backtrace? I would like to understand the cause behind the inability to perform the dot product operator.

@davelee2804
Copy link
Collaborator

err.log

@amartinhuertas
Copy link
Collaborator Author

Ok, digging into the error log, I saw this line of code:

setindex!(A::Vector{Gridap.TensorValues.VectorValue{2, Float64}}, x::Gridap.TensorValues.VectorValue{3, Float64}, i1::Int64)

There we are trying to set a 3-components field vector (Gridap.TensorValues.VectorValue{3, Float64}) as an entry of a Vector of 2-component field vectors (Vector{Gridap.TensorValues.VectorValue{2, Float64}}). Thus, as you say, it seems we are messing things up with the number of spatial dimensions.

What is the error that you get with D=3?

@davelee2804
Copy link
Collaborator

davelee2804 commented Jun 30, 2022

Well spotted! Here is the D=3 log file:
err-2.log

@amartinhuertas
Copy link
Collaborator Author

ok. Use D=3 ONLY for this line, i.e.,

reffeᵤ = ReferenceFE(lagrangian,VectorValue{3,Float64},order;space=:P)

The rest untouched.

@davelee2804
Copy link
Collaborator

...Same error as the original, as far as I can see: :(
err-3.log

@davelee2804
Copy link
Collaborator

davelee2804 commented Jun 30, 2022

...I just got this error:
err-4.log
...with this source:
AdvectionHDG.log
(so un is here projected onto U, but is maybe not being interpreted as a vector field)?

There is this:
Vector{Gridap.ReferenceFEs.GenericRefFE{Gridap.ReferenceFEs.RaviartThomas, 2}}
But we don't want an RT ReferenceFE, we want an L2 vector field!?!

not crash). The current results it is generating are wrong.
@amartinhuertas
Copy link
Collaborator Author

Hey @amartinhuertas , I'm having some basic trouble with the HDG advection on the sphere. I am getting the error:

For the records, this was addressed in 1586cf7

@amartinhuertas
Copy link
Collaborator Author

Hi @davelee2804,

please read the below carefully, let me know if you have questions, and please answer my questions (no hurries).

Current status of this development:

  • I have solved several minor issues in GridapHybrid.jl when dealing with manifolds. (see https://github.com/gridap/GridapHybrid.jl/pull/23/files). These have been already merged to GridapHybrid.jl#main. As we are pointing there in the GridapGeosciences.jl#hdg_adv branch, after pulling from remote hdg_adv, you just have to write up in the Julia REPL package manager to update GridapHybrid.jl.
  • Now the HDG driver for advection on the sphere runs bottom-top (it does not crash). HOWEVER it does not produce the correct result, in particular, it does not seem to time evolve.

Important remarks:

  • Let us stick (by now) to a 2nd-order polynomial approximation of the sphere (note that in the latest version of the driver I am passing polynomial degree 2 in this call https://github.com/gridapapps/GridapGeosciences.jl/pull/34/files#diff-d6ace757a93640319c464479d34d4617d5d53c419c630a9afee6b9c8db9b8d46R41 ). I still have to go through the analytical mapping cubed sphere mesh to see whether it properly represents the 1D skeleton in 3D ambient space.
  • I had to change the order of the operands in terms like this m*p*((un⋅nₒ) + τ*(n⋅nₒ)) as ((un⋅nₒ) + τ*(n⋅nₒ))*m*p. In the former form, GridapHybrid.jl generates an internal error (I will solve it when I have time, I think there are some functions missing, but I still have to fully understand this). However, in the latter form, it does not.

Questions:

  • Do you have a reference where the formulation you are aiming to implement is written? If not, can you write it? This may help to (1) try to understand what you are trying to do. (2) try to determine if the code is doing what you are trying to do or not.
  • I have seen you have terms like n \cdot n. Given that n is a unit vector, isnt n \cdot n always 1?
  • Note that no and n are defined differently. n is the tangent outward unit normal (it always points outwards). no points outwards if the corresponding cell is the OWNER of the edge, and points inwards if the corresponding cell is not the owner. Thus, for a given cell K, and edge E in the bounday of K, n \cdot no is 1 if K owns E, and it is not 1 (its value depends on the cosine of the angle among the two normals) if it does not own it. We assign the owner cell of an edge arbitrarily to be the cell with the minimum global identifier among the two cells surrounding the edge.

@davelee2804
Copy link
Collaborator

Thank you @amartinhuertas ! I had not appreciated the subtlety that in for example Kang et al. eqns (14) and (15) the first skeleton integral is over \partial K while the second is over e, this makes much more sense now! Looking at (7) in Nguyen et al I see that all the skeleton integrals are over \partial T, which is consistent with

integrand restricted to each edge seen from the perspective of K

as you say.

I am now implementing a convergence study for the solid body rotation test (second order in space + time), to verify that everything is OK.

A few days ago I got some nonsensical results with order 4, so I will also double check that as well...

@amartinhuertas
Copy link
Collaborator Author

Thank you @amartinhuertas ! I had not appreciated the subtlety that in for example Kang et al. eqns (14) and (15) the first skeleton integral is over \partial K while the second is over e, this makes much more sense now! Looking at (7) in Nguyen et al I see that all the skeleton integrals are over \partial T, which is consistent with

Exactly. You can write hybridizable methods either by edges or grouping together edges into cell boundaries \partial K, for each K. The second approach is the one required to implement static condensation locally at each cell, towards assembly of the global problem on the skeleton, which is the main motivation for hybridizable methods.

I am now implementing a convergence study for the solid body rotation test (second order in space + time), to verify that everything is OK.
A few days ago I got some nonsensical results with order 4, so I will also double check that as well...

Ok. Let us see.

@davelee2804
Copy link
Collaborator

Hey @amartinhuertas , just a heads up that the convergence is sub-optimal. There is too much dispersion, perhaps as a result of how I have formulated the time stepping. I'll take a closer look at this...

@amartinhuertas
Copy link
Collaborator Author

Hey @amartinhuertas , just a heads up that the convergence is sub-optimal. There is too much dispersion, perhaps as a result of how I have formulated the time stepping. I'll take a closer look at this...

I would try to first solve a steady advection problem with manufactured PDE solution. This way we eliminate noise related to time discretization.

also removing the (n.n) terms. results greatly improved
@davelee2804
Copy link
Collaborator

davelee2804 commented Jul 10, 2022

Thanks @amartinhuertas , yes great point, this is the correct approach, I will do this.

Just a heads up, I changed the time integration to a stiffly-stable second order implicit runge-kutta scheme, and now the solid body rotation test is converging at the correct rate (second order convergence for second order in space+time discretisation) after a single revolution round the sphere (at least for the low resolutions that I have checked).

The entropy conservation error is also converging quadratically, so that is also a good sign...

@davelee2804
Copy link
Collaborator

hdg_convergence

@davelee2804
Copy link
Collaborator

Hi @amartinhuertas , I've added a convergence test for the solid body rotation configuration of the HDG advection solver. I'm wondering if now is a good time to do a code review and merge this branch, before starting work on the linear wave equation with jump penalisation on the pressure gradient tangent component - what do you think?

@davelee2804
Copy link
Collaborator

...SBRAdvectionHDG tests are failing in github-ci (but passing on my laptop) :
https://github.com/gridapapps/GridapGeosciences.jl/runs/7274708478?check_suite_focus=true#step:9:1037

Should I also commit the Manifest.toml?

@amartinhuertas
Copy link
Collaborator Author

Should I also commit the Manifest.toml?

Yes, because we are pointing to versions of some of the packages which are not yet registered in the Julia registries.

@amartinhuertas
Copy link
Collaborator Author

Also, please, pull from master if you didnt lately. I see in this PR the following changes that are already in master, so that they should not be in this PR.

https://github.com/gridapapps/GridapGeosciences.jl/pull/34/files#diff-8ad2a593c8bedda9e0bea392ea69c0bbbe41609839d4d3b33bc4a81155a208dcL48

@amartinhuertas
Copy link
Collaborator Author

Hi @amartinhuertas , I've added a convergence test for the solid body rotation configuration of the HDG advection solver.

I see in the attached plot that there is a mild degradation in slope for the highest resolution that you tested. Have you tried even higher resolutions to be 100% sure that the convergence is actually quadratic?

@davelee2804
Copy link
Collaborator

I see in the attached plot that there is a mild degradation in slope for the highest resolution that you tested

Hi @amartinhuertas , higher resolutions are to the left on the x-axis, so the convergence rate is actually increasing as resolution increases. Does that make sense? Have I missed your point somehow?

@davelee2804
Copy link
Collaborator

...Ie: convergence rate from dx=0.12 to 0.06 is higher than from dx=0.24 to 0.12. Does that make sense?

@amartinhuertas
Copy link
Collaborator Author

...Ie: convergence rate from dx=0.12 to 0.06 is higher than from dx=0.24 to 0.12. Does that make sense?

Ok, you are right, I was confused. The convergence curve looks like a piece-wise linear function with 2 pieces, where the right piece seems to be parallel to the theoretical convergence, and the left piece seems to bend down a little bit, with an slope which is no longer parallel to the theoretical convergence. But, as you say, bending down means even faster convergence, and not the opposite. (Here it was my confusion).

@amartinhuertas
Copy link
Collaborator Author

I'm wondering if now is a good time to do a code review and merge this branch

Yes, it might be a good time to do this. The only thing is that I have some other code reviews on the queue, I wont be able to do this immediately.

Apart from this code review, we need:

  1. To investigate issue Fix failing MPI tests #36
  2. To fix (if actually required) and test the geometrical discretization of the sphere using an analytical mapping.

@davelee2804
Copy link
Collaborator

No problem. Note that for the SBRAdvectionHDGTests file that I committed, convergence between the lowest two resolutions is sub-optimal, but quadratic convergence is achieved between the next two.

@davelee2804
Copy link
Collaborator

No problem @amartinhuertas , there is no rush from me. Please let me know if there is anything I can do to help with any of the issues you mention. We are using an isoparametric mapping right now, is that correct?

@amartinhuertas
Copy link
Collaborator Author

We are using an isoparametric mapping right now, is that correct?

We are using biquadratic elements for the geometry, and bilinear elements for the unknown. I think that the term "isoparametric" is used to refer to those scenarios in which the order of the polynomial to approximate the geometry matches the order of the polynomial to approximate the unknown.

@davelee2804
Copy link
Collaborator

Aaah, yes, that is what the 2 is for in the CubedSphereDiscreteModel :) OK. Is there a particular reason we need to support an analytical mapping as well?

@amartinhuertas
Copy link
Collaborator Author

Is there a particular reason we need to support an analytical mapping as well?

We have been using an analytical mapping in all the results we obtained with GridapGeosciences + compatible FEs so far.

If we can represent the geometry exactly, better to stick to it, instead of approximating it, so that we eliminate any noise caused by geometrical errors (leaving integration errors apart). I remember e.g., that, with a bilinear approximation of the cubed sphere we were not able to obtained the theoretical convergence rates of RT+DG FEs for Darcy on the sphere.

@davelee2804
Copy link
Collaborator

OK, fair enough. To my mind supporting iso or super parametric mappings (as you are already doing) is better, as if we ever want to add mountains when we will have to ditch the analytic mapping anyway, so ensuring everything works for a numerical mapping makes it easier to transition to such a case if we ever want to. But I take your point, analytic is nicer...

@amartinhuertas
Copy link
Collaborator Author

Yes, sure, better to have both.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants