HDG for advection-diffusion on the sphere

[amartinhuertas]

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

[codecov-commenter]

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]

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]

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]

(Note that D = 2).

Did you try with D=3?

[davelee2804]

...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]

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]

err.log

[amartinhuertas]

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]

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

[amartinhuertas]

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

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

The rest untouched.

[davelee2804]

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

[davelee2804]

...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!?!

[amartinhuertas]

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 https://github.com/gridapapps/GridapGeosciences.jl/pull/34/commits/1586cf7729879573c9598fe9dcf56acfad16db25

[amartinhuertas]

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]

Thank you @amartinhuertas , that was quick!

I'm not surprised there are still bugs in HGD advection the code. I update to your current version of GridapHybrid.jl and take a look (hopefully tomorrow).

Let us stick (by now) to a 2nd-order polynomial approximation of the sphere

Great idea, yes, let's do this

Do you have a reference where the formulation you are aiming to implement is written?

It is solid body rotation of the flux form advection equation, so the (constant) velocity field is 0 in the north-south direction and cos(\phi) in the east-west direction (with \phi being the latitude). This is the 'Williamson 1' test case ('Advection of a Cosine Bell) from this paper here: https://www.sciencedirect.com/science/article/pii/S0021999105800166 (with \alpha=0). However I am using a Gaussian not a cosine as this is infinitely differentiable, so a little nicer to work with. The weak HDG form of the flux form advection equation is based on the formulation in Kang et. al. JCP 2020

Given that n is a unit vector, isnt n \cdot n always 1

Yes, as you say. I have kept the n\cdot n as I wanted to be consistent with respect to the n\cdot n_o terms in the jumps, just until I am comfortable with the method, then I will remove

Note that no and n are defined differently

Yes, I noticed this, thanks. So n\cdot n_o represents a jump, ie a [[ . ]] operator yes?

it does not produce the correct result, in particular, it does not seem to time evolve

Yes, it seems to be coming out as 0 everywhere. I will look into this tomorrow...

[amartinhuertas]

I'm not surprised there are still bugs in HGD advection the code. I update to your current version of GridapHybrid.jl and take a look (hopefully tomorrow).

Please note that the bugs might not come only from your side. We are testing GridapHybrid with manifolds for the first time and there might be still pitfalls/bugs to be determined/fixed. I think you are aware of this, but just to be 100% sure. I would be extremely careful, and may be try to think in simpler tests we can write to determine whether integration of fields over \partialK is actually working as expected with manifolds.

The weak HDG form of the flux form advection equation is based on the formulation in Kang et. al. JCP 2020

Can you send me the particular Eq. number in that paper?

Yes, I noticed this, thanks. So n\cdot n_o represents a jump, ie a [[ . ]] operator yes?

Please note that, in general, the jump operator may come in two different flavours:

• u^+ - u^-, where u can be a scalar or vector double-valued discontinuous field. If u is scalar, then the jump is scalar, if it is a vector, then it is a vector.
• u^+ n^+ + u^- n^-. Where n is the outward unit normal as seen from the perspective of the two cells around. If u is an scalar, the jump is a vector. If u is a vector, the jump is an scalar. Note that for manifolds n^+ is not equal to -n^- in general.

See def 14.4 and remark 14.5 in https://www.math.tamu.edu/~guermond/M610_SPRING_2016/chap14.pdf

n \cdot n0 seems actually another thing (like trying to give a global orientation to the normal, and then adjusting if the normal from the current cell does not match the one given in the global orientation), but perhaps I am missing something. I am not saying it is incorrect, I am bit confused actually.

[davelee2804]

Hey @amartinhuertas , the relevant equations in Kang et. al. are (3a), (16) and (17).

With respect to the normals we also have that n^+ = -n^-, so assuming K^+ is the cell owning the normal, then n^+ = n_o and n^- = -n_o so u\cdot n\_o integrated over \partial K from both sides would be (n^+ - n^-)\cdot n_o (ie: a scalar jump). At least that is my understanding. I also get very confused with these binary choices (that is why I left the n\cdot n terms in the formulation :) ), there is a very strong chance that I am wrong...

[davelee2804]

Hey @amartinhuertas , just a heads up, when I check the sum() of the input fields un, pn, these are non-zero, however the half step field, ph has zero sum(). I've also looked at individual matrix terms as vectors (assuming pn for the trial function), and these also have non-zero sum()

[amartinhuertas]

u\cdot n_o integrated over \partial K from both sides would be (n^+ - n^-)\cdot n_o

Why u\cdot n_o integrated over both sides and then summed is (n^+ - n^-)\cdot n_o? In the latter expression, I do not see the u. ??? What am I missing ...

With respect to the normals we also have that n^+ = -n^-

I have doubts in regards to this equivalence for discretized manifolds. I think this equivalence is only approximate in this case, as the normals might be contained in different planes (the ones defined by the cells). At least the code does not assume this, and compute n^+ and n^- from each cell without relying on this equivalence. I guess that the finer the mesh, the more accurate the approximation, but not totally sure. I can write an small numerical experiment to test whether this is true or not.

[amartinhuertas]

Hey @amartinhuertas , just a heads up, when I check the sum() of the input fields un, pn, these are non-zero, however the half step field, ph has zero sum(). I've also looked at individual matrix terms as vectors (assuming pn for the trial function), and these also have non-zero sum()

Thanks for reporting this! I found another BUG for GridapHybrid.jl + manifolds. I have fixed it in this PR https://github.com/gridap/GridapHybrid.jl/pull/24/files. I have merged it into main. Please execute up in the package manager to get this BUG fix on your env.

I have executed the HDG advection driver after the bug fix, and it now time evolves, but it does seem to generate just rubbish. We still have to find extra bug fixes, formulation pitfalls etc.

[amartinhuertas]

@davelee2804 ... for convenience, I wrote a MWE (Minimum Working Example) to find the previous BUG. I am copying it here so that you can also leverage it if require. I think it can be highly helpful as a template to write MWEs for other yet-to-spot BUGs.

using LinearAlgebra
using SparseArrays
using Gridap
using GridapGeosciences
using GridapHybrid

function Gridap.Geometry.push_normal(invJt,n)
  v = invJt⋅n
  m = sqrt(inner(v,v))
  if m < eps()
    return zero(v)
  else
    return v/m
  end
end

function project_initial_conditions(dΩ, P, Q, p₀, U, V, u₀, mass_matrix_solver)
  # the tracer
  b₁(q)    = ∫(q*p₀)dΩ
  a₁(p,q)  = ∫(q*p)dΩ
  rhs₁     = assemble_vector(b₁, Q)
  L2MM     = assemble_matrix(a₁, P, Q)
  L2MMchol = numerical_setup(symbolic_setup(mass_matrix_solver,L2MM),L2MM)
  pn       = FEFunction(Q, copy(rhs₁))
  pnv      = get_free_dof_values(pn)

  solve!(pnv, L2MMchol, pnv)

  # the velocity field
  b₂(v)    = ∫(v⋅u₀)dΩ
  a₂(u,v)  = ∫(v⋅u)dΩ
  rhs₂     = assemble_vector(b₂, V)
  MM       = assemble_matrix(a₂, U, V)
  MMchol   = numerical_setup(symbolic_setup(mass_matrix_solver,MM),MM)
  un       = FEFunction(V, copy(rhs₂))
  unv      = get_free_dof_values(un)

  solve!(unv, MMchol, unv)

  pn, pnv, L2MM, un
end

order  = 1
degree = 3
n      = 16
dt     = 0.25*2.0*π/((order+1)*4*n)
# single rotation about the sphere
nstep  = Int(2.0*π/dt)

# solid body rotation velocity field
function u₀(xyz)
  θϕr = xyz2θϕr(xyz)
  u   = cos(θϕr[2])
  spherical_to_cartesian_matrix(θϕr)⋅VectorValue(u,0,0)
end

# Gaussian tracer initial condition
function p₀(xyz)
  rsq = (xyz[1] - 1.0)*(xyz[1] - 1.0) + xyz[2]*xyz[2] + xyz[3]*xyz[3]
  exp(-4.0*rsq)
end

model = CubedSphereDiscreteModel(1,2; radius=1)

# Forward integration of the advection equation
D = num_cell_dims(model)
Ω = Triangulation(ReferenceFE{D},model)
Γ = Triangulation(ReferenceFE{D-1},model)
∂K = GridapHybrid.Skeleton(model)

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

# Define test FESpaces
V = TestFESpace(Ω, reffeᵤ; conformity=:L2)
Q = TestFESpace(Ω, reffeₚ; conformity=:L2)
M = TestFESpace(Γ, reffeₗ; conformity=:L2)
Y = MultiFieldFESpace([Q,M])

U = TrialFESpace(V)
P = TrialFESpace(Q)
L = TrialFESpace(M)
X = MultiFieldFESpace([P,L])

dΩ  = Measure(Ω,degree)
d∂K = Measure(∂K,degree)

# HDG stabilisation parameter
τ = 1.0

# Project the initial conditions onto the trial spaces
pn, pnv, L2MM, un = project_initial_conditions(dΩ, P, Q, p₀, U, V, u₀, Gridap.Algebra.BackslashSolver())

dth = 0.5*dt

n  = get_cell_normal_vector(∂K)
nₒ = get_cell_owner_normal_vector(∂K)

# First stage
b₁((q,m)) = ∫(q*pn)dΩ - ∫(m*0.0)d∂K
a₁((p,l),(q,m)) = ∫(q*p - dt*(∇(q)⋅un)*p)dΩ + ∫(((un⋅n) + τ*(n⋅n))*dt*q*p)d∂K -   # [q,p] block
              ∫(dt*q*l*τ*(n⋅n))d∂K +                                              # [q,l] block
              ∫(((un⋅nₒ) + τ*(n⋅nₒ))*m*p)d∂K -                                    # [m,p] block
              ∫(m*l*τ*(n⋅nₒ))d∂K                                                  # [m,l] block

a11(p,q)=∫(q*p - dt*(∇(q)⋅un)*p)dΩ + ∫(((un⋅n) + τ*(n⋅n))*dt*q*p)d∂K
A11=assemble_matrix(a11,X[1],Y[1])
a12(l,q)=∫(-dt*q*l*τ*(n⋅n))d∂K
A12=assemble_matrix(a12,X[2],Y[1])
a21(p,m)=∫(((un⋅nₒ) + τ*(n⋅nₒ))*m*p)d∂K
A21=assemble_matrix(a21,X[1],Y[2])
a22(l,m)=∫(-1.0*m*l*τ*(n⋅nₒ))d∂K
A22=assemble_matrix(a22,X[2],Y[2])

b1(q)= ∫(q*pn)dΩ
b2(m)=∫(-m*0.0)d∂K

A=assemble_matrix(a₁,X,Y)

B1=assemble_vector(b1,Y[1])
B2=assemble_vector(b2,Y[2])

BΓ  = B2  - Matrix(A21)*(A11\B1)
BΓΓ = A22 - Matrix(A21)*inv(Matrix(A11))*Matrix(A12)

op₁  = HybridAffineFEOperator((u,v)->(a₁(u,v),b₁(v)), X, Y, [1], [2])
Xh   = solve(op₁)
ph,_ = Xh

op_BΓΓ = op₁.skeleton_op.op.matrix
op_BΓ = op₁.skeleton_op.op.vector

@assert norm(op_BΓΓ-BΓΓ) < 1.0e-12
norm(BΓ-op_BΓ) < 1.0e-12

[davelee2804]

Hey @amartinhuertas

Why u\cdot n_o integrated over both sides and then summed is (n^+ - n^-)\cdot n_o? In the latter expression, I do not see the u. ??? What am I missing ...

Sorry, silly typo, that should read (u^+ - u^-)\cdot n_o, silly typo on my part

I think this equivalence is only approximate in this case, as the normals might be contained in different planes

Great point, yes, I think as one refines the mesh n^+ converges to -n^-, but it will probably not be exact (the H(div) Piola map gives this exactly though yes?)

[davelee2804]

I have executed the HDG advection driver after the bug fix, and it now time evolves, but it does seem to generate just rubbish. We still have to find extra bug fixes, formulation pitfalls etc.

Thank you for identifying the error @amartinhuertas ! The subsequent errors are probably on my end, I will do some digging...

..and thanks for the example, yes it does look super helpful. Now I will understand how to interface with and replicate the hybrid solver a little better...

[amartinhuertas]

hi @davelee2804 ... I am trying to make sense out of the HDG formulation that you have written, but I am struggling quite a bit (I cannot say whether it is correct or not, just that I dont understand).

Looking at https://arxiv.org/abs/1605.03228, Eqs. (3.1)-(3.3), they briefly present the HDG upwind formulation for the (steady-state) advection equation (originally proposed in ref [5] https://users.oden.utexas.edu/~tanbui/PublishedPapers/phoga.pdf). Is this the formulation that are you trying to implement? (see the snapshot below).

For convenience, the correspondence among the article notation and the code notation is:

• \beta -> un
• u -> p
• v -> q
• \hat{u} -> l
• \mu -> m

By now, I have sevaral concerns:

1. The article mentions that beta is assumed to be a continuous function. However, un being used in the code is the projection of a continuous function into a discontinuous FE space. Correct? I would like to understand why beta has to be continuous.
2. The term which is integrated over the bulk (i.e., K) in the article does not match the one in the code.
3. I do not know how to proof the equivalence among the terms in the code and the terms in the article over \partial K.

Perhaps you are trying to implement a different formulation ... not sure ....

[davelee2804]

Hey @amartinhuertas ,

Is this the formulation that are you trying to implement?

They are doing material form advection \beta.\nabla u, whereas I am doing flux form \nabla\cdot(\beta u() For a divergence free vector field, \nabla\cdot\beta = 0, the two are the same. The velocity field we are using (\cos(\phi),0) is divergence free in the continuous form, but for \beta\in H1(\Omega I would imagine there will be some numerical divergence, so they won't be exactly the same

However, un being used in the code is the projection of a continuous function into a discontinuous FE space. Correct? I would like to understand why beta has to be continuous.

We could do un continuous or discrete. We are only solving a single equation u -> p so there is the velocity field \beta->un will always be prescribed. When we extend to the shallow water equations then we solve for \beta->un

The term which is integrated over the bulk (i.e., K) in the article does not match the one in the code.

Indeed, that is a result of the marerial vs. flux form advection.

I don't understand the differences in the mesh skeleton term w.r.t. Kang et al 2020. (ie the use of |\bea\cdot n|) - do you have any thoughts on this??

[amartinhuertas]

They are doing material form advection \beta.\nabla u, whereas I am doing flux form \nabla\cdot(\beta u) For a divergence free vector field, \nabla\cdot\beta = 0, the two are the same.

Ok, I now understand that if the fluid flow is incompressible, then both the term in the code and the one in the article are equivalent, as I have double checked in the following note. Agreed?

[amartinhuertas]

We could do un continuous or discrete.

I would do it continuous (as far as possible). Is there any reason you need to interpolate it into a FE space? The only concern is that if we have a differential operator applied to the continuous solution we have to ensure that we apply the spherical variant of the differential operator, as we do, e.g., in the Darcy test with the laplacian applied to the analytical solution.

We are only solving a single equation u -> p so there is the velocity field \beta->un will always be prescribed. When we extend to the shallow water equations then we solve for \beta->un

Yes, I was aware of this. \beta is input data to the advection problem.

[amartinhuertas]

I don't understand the differences in the mesh skeleton term w.r.t. Kang et al 2020. (ie the use of |\bea\cdot n|) - do you have any thoughts on this??

I think that Kang et al and the iHDG paper (the one from which I extracted the capture above) use the same abstract family of HDG formulations for hyperbolic problems of the general form

For the advection problem (3.1), we have that

F_k(u)=\beta_k u 
C=0

Thus, A_k is a 3x3 tensor with entries defined as [A_k]_kk=\beta_k and 0 otherwise.

Agreed? (please let me know if everything clear)

[amartinhuertas]

The general form of the upwind HDG discretization for problem (2.1) reads as:

In the case of problem (3.1), A is diagonal matrix, with \beta_k n_k in the diagonal, for k=1..3. Thus, after computing the eigendecomposition of A, R is the identity, and S is equal to A, and |A| is |\beta \cdot n|, which is precisely what they use in (3.2) and (3.3). This choice of the upwind flux is the solution of a Riemman problem on the interface of the cells (TBH, I still do not fully understand this part, I need to devote more time to understand where the definition of the \hat{F} flux comes from and why it provides the desired convergence, stability, etc. properties.).

I guess that in Kang et. al. they are leveraging the eigenvalues of the A matrix that appears for the shallow water equations, and thus this does not apply verbatim to the (much simpler) advection equation.

How did you derive the formulation in the code for the advection equation from the one of the shallow water equations in Kang et. al.? Could this be what is wrong?

[davelee2804]

Hey @amartinhuertas , I think I agree with what you say. If I unroll the definition of \hat{F}\cdot n from eqn (3) in Muralikrishnan et al as used in section (4) of that paper I get:

\hat{F}\cdot n = \beta\cdot n u + |\beta\cdot n|(u - \hat{u}), as in eqn (19), where \beta is the trial eigenvalue for the advection problem

if I unroll (17) from Kang et al (which is what I was trying to do) I get (using the same nomenclature) the same (using n\cdot n = 1)

Based on this I think my current implementation is wrong in two ways:

1. I am not using an absolute value, and
2. I am not using the correct \tau (I am just using \tau=1 for convenience, based on the max value of the absolute value of |\beta\cdot n| for this simple problem.

Is there a way to apply an absolute value in Gridap? I know there is mean() and jump(), is there also an abs() ?

[amartinhuertas]

I know there is mean() and jump(), is there also an abs() ?

I think that is possible. AFAIK, the high level API supports an invokation to an arbitrary Julia function passing a CellField as an argument (\beta \cdot n is a CellField). Anyway, I have never tried that with GridapHybrid.jl ... let me know if that does not work.

[amartinhuertas]

BTW, did you notice that in (2.2a) they are using exactly the same bulk term as the one you are using in your code? :-)

[davelee2804]

BTW, did you notice that in (2.2a) they are using exactly the same bulk term as the one you are using in your code? :-)

I haven't checked it experimentally, but for a divergence free \beta I think it should be OK.

I've just made a change that leads to a much more coherent solution (and exact mass conservation), however, there is no jump() in there now as this leads to a code error as: It is not possible to perform the operation "*" on the given cell fields.

Do you know if the jump() operator is supported for integration with respect to test functions on the Skeleton using the measure d\partial K ?

[amartinhuertas]

Do you know if the jump() operator is supported for integration with respect to test functions on the Skeleton using the measure d\partial K ?

Yes. The jump() operator is not supported, on purpose. The measure d\partial K represents the boundary of the cells, on a cell by cell basis, such that each cell is independent of each other. The jump() operator can be used only with facet-wise triangulations that glue together the two cells around. I am not sure why you need the jump() operator explicity.

[davelee2804]

OK, great to know, thanks for the confirmation @amartinhuertas . Well, how do we represent [[.]] then? I guess this is what n\_o is for... I'll keep trying to find a way to replicate the jump() with this. Right now the results look coherant and conservative, but I think over long times they may go oscillatory without some sort of jump penalisation...

[amartinhuertas]

Well, how do we represent [[.]] then?

Let us assume that we want to compute, e.g.,

< [[u]], v>_e, for all edges e in the mesh, where [[u]] = u^+ - u^-

Then you have to write

< u*(n \cdot n_o), v>_d\partial K, for all cells K in the mesh.

So I think that you already got that from what I see in your comments.

[amartinhuertas]

Right now the results look coherant and conservative, but I think over long times they may go oscillatory without some sort of jump penalisation...

This is great! I am surprised that we only needed a moderate amount of BUG fixes. :-)

[amartinhuertas]

I haven't checked it experimentally, but for a divergence free \beta I think it should be OK.

I am surprised they do not mention that \beta is assumed to be divergence free at any point in the iHDG paper ... Perhaps it is something well-known in hyperbolic fluid flow problems or something like that ...

[davelee2804]

This is great! I am surprised that we only needed a moderate amount of BUG fixes. :-)

Doesn't mean it is working correctly yet however :). Could be that it will still go oscillatory or not advect correctly (I only checked a couple of time steps).

When I add in the jumps as n\cdot n\_o I get completely nonsense results (and a total loss of mass conservation), so I am wondering if there is some subtle interplay between the n\cdot n_o jump and the |A| absolute value that I don't properly understand....

[davelee2804]

I am surprised they do not mention that \beta is assumed to be divergence free at any point in the iHDG paper ... Perhaps it is something well-known in hyperbolic fluid flow problems or something like that ...

\beta could have a divergence. It is just that if \beta has non-zero divergence then the material and flux form advection equations are no longer equivalent. I thought I would do flux form as that is how the advection equation is done in shallow water.

[amartinhuertas]

\beta could have a divergence. It is just that if \beta has non-zero divergence then the material and flux form advection equations are no longer equivalent. I thought I would do flux form as that is how the advection equation is done in shallow water.

Ok. Thanks for the clarification. I think I 100% understand now.

[amartinhuertas]

Doesn't mean it is working correctly yet however :). Could be that it will still go oscillatory or not advect correctly (I only checked a couple of time steps). When I add in the jumps as n\cdot n_o I get completely nonsense results (and a total loss of mass conservation), so I am wondering if there is some subtle interplay between the n\cdot n_o jump and the |A| absolute value that I don't properly understand..

To which kind of oscillations are you referring to? You mean spurious modes that you usually have in the Galerkin method if you do not use upwinding? (such as SUPG?).

I am not wrong, the formulation we are leveraging already comes with upwinding "batteries-on", no need to do anything else. But perhaps I am missing something ...

[davelee2804]

Yes, the formulation we are trying to implement includes upwinding, agreed. However I am not able to get the jump penalisation working, which is why I think it will eventually go oscillatory... (ie: the n\cdot n\_o term leads to large errors...)

[davelee2804]

...Does the integral over the measure d\partial K already assume that the jump is baked in?

[amartinhuertas]

I am not able to get the jump penalisation working

Which is the evidence that supports this statement?

[amartinhuertas]

...Does the integral over the measure d\partial K already assume that the jump is baked in?

No, AFAIK.

[amartinhuertas]

Have you noticed that both Kang et al and MURALIKRISHNAN et al defined the jump as [[u]] = u^+ + u^- ??? With a + not a -

[amartinhuertas]

Which is the evidence that supports this statement?

Ok, I think you already answered to that, sorry. Namely here

ie: the n\cdot n_o term leads to large errors...

I am not sure we have to use n\cdot n_o here, given the definition of the jump operator in that papers.

[davelee2804]

Hey @amartinhuertas , as mentioned, the application of 'n\cdot n_o` leads to erroneous results for me. With the current configuration the results are correct (conservative and advecting at the correct speed), see images attached. So I'm very confused about the application of the jumps. In this article (which also discusses continuity of the tangent vectors), there are no jumps explicitly mentioned (see eqn 7 for example): https://www.sciencedirect.com/science/article/pii/S0021999111003226

:

[davelee2804]

...so I have no idea if the penalisation is properly applied or not. Given that the results are non-oscillatory and conservative, it looks like maybe this is being done correctly, but I don't know enough about how this works under the hood to say. Also I will need to implement a more sophisticated version of \tau (as |u\cdot n|, not 1).