Comprehensive set of tests involving terms with divergence

[amartinhuertas]

@davelee2804 @santiagobadia @cbladwell

IMPORTANT NOTE: Please, this PR is highly important, let us try to process it with care and solve it to the best of our abilities. (no rush, let us give it a profound thinking). If we suceed I think we are going to avoid a lot of pain down the road

In this PR, I have written code to perform systematic tests of variational formulation terms involving the divergence operator. I have written them, among others (although not restricted to), to shed some light to the discussion that we have had in PR https://github.com/gridapapps/GridapGeosciences.jl/pull/10 regarding this particular point:

Solve the following question: Can we actually exploit the DIV+reference integration for this term: (∇⋅(v))(0.5ui⋅ui) ?

You can see the code I have written here: https://github.com/gridapapps/GridapGeosciences.jl/blob/79363b97fb704592a9bc4a193f500d017ea388cc/test/DivTermsTests.jl

As always, number one task is to confirm that the code I have written is correct, that I did not fall into confusion/semantic errors, etc. In my view, I think it is correct, but it is always best to have more eyes to double check I did not do anything silly.

According to this code, the following conclusions can be extracted (you can run the code yourself and play around, but you can find the results of this script on my machine below):

1. In the discrete cubed sphere manifold., ∫(∇⋅(v)*p)dΩ and ∫(DIV(v)*p)dω are equivalent to machine precision, no matter mesh resolution, nor polynomial degree. This is something I have observed before (I have told it to you in the past, indeed, not sure if I stressed that enough to transfer the message). This is not new for me, but it is still puzzling me. How come this can this be true if we are computing ∇⋅(v) assuming that the Jacobian is constant within each cell? Any ideas? May it be related with the particular properties of the discrete functions v and p? (hypothesis, no idea)
2. ∫(DIV(v)*p)dω is within 7-10 times faster than ∫(∇⋅(v)*p)dΩ for the tests I have performed (not surprising).
3. ∫(∇⋅(v)*(u⋅u))dΩ and ∫(∇⋅(v)*(u_dot_u_projection))dΩ are equivalent to machine precision in a 2D problem, but not in the discrete cubed sphere manifold. Why????!!!!!
4. ∫(∇⋅(v)*(u_dot_u_projection))dΩ and ∫(DIV(v)*(u_dot_u_projection))dω are equivalent to machine precision (both in a 2D problem and the discrete cubed sphere manifold). This is reasonable and consistent with point 1. above
5. ∫(∇⋅(v)*(u⋅u))dΩ and ∫(DIV(v)*(u⋅u))dω are also equivalent to machine precision both in a 2D problem and the discrete cubed sphere manifold. I do not understand how can this be (taking into account how u is built in Gridap, with the Piola map, etc.)

Thanks for your help!

setup_FE_spaces (generic function with 1 method)

compute_relative_error (generic function with 1 method)

compare_phys_vs_ref_div_v_p_term (generic function with 1 method)

compare_div_v_u_dot_u_versus_div_v_projected_u_dot_u (generic function with 1 method)

compare_div_v_projected_u_dot_u_DIV_v_projected_u_dot_u (generic function with 1 method)

compare_div_v_u_dot_u_DIV_v_u_dot_u (generic function with 1 method)

1.0e-15

5.710127560402143
8.132490888663972
8.105644019093509
9.605480986715655
10.45951241083886
8.680363007795998
11.065744086462233
9.750486202515491
11.232162339079204
11.005677659262158
17.075798750138524
 51.839175 seconds (102.05 M allocations: 5.242 GiB, 6.70% gc time)
([1.0, 0.6666666666666666, 0.5, 0.4, 0.3333333333333333, 0.2857142857142857, 0.25, 0.2222222222222222, 0.2, 0.2, 0.1], [8.171020218150498e-17, 1.0894693624200663e-16, 1.387778780781446e-16, 8.973825869789967e-17, 1.1179063206240823e-16, 1.1139928929424335e-16, 9.41237363483903e-17, 1.1970611027533865e-16, 1.1018648905114065e-16, 1.1018648905114065e-16, 1.0300771835293419e-16], -0.05249283547206647)

Test Passed

7.905290073347969
7.675951461665748
4.659955384726491
5.826073435239282
7.781720451676158
7.613958125824359
8.832022609239253
8.38297670371396
8.72660379056484
8.991955579911698
12.055192135708632
  3.579519 seconds (8.37 M allocations: 583.297 MiB, 54.82% gc time)
([1.0, 0.6666666666666666, 0.5, 0.4, 0.3333333333333333, 0.2857142857142857, 0.25, 0.2222222222222222, 0.2, 0.2, 0.1], [2.50629346356623e-16, 2.7397467455750404e-16, 3.6398400866098096e-16, 2.6124188738696184e-16, 2.788804597125794e-16, 2.845262667442368e-16, 2.711415138006175e-16, 2.8219344409845557e-16, 2.7507525863244907e-16, 2.7507525863244907e-16, 2.749828919728138e-16], 1.554410272138009e-5)

Test Passed

6.743851830756222
7.632346243410186
6.926109458881564
5.470265339261303
3.9082712600892737
8.17032365423259
7.851674052373194
7.448766248970142
7.2457589988461075
7.656558387991294
7.641623593559819
  5.299253 seconds (21.48 M allocations: 1.186 GiB, 39.14% gc time)
([1.0, 0.6666666666666666, 0.5, 0.4, 0.3333333333333333, 0.2857142857142857, 0.25, 0.2222222222222222, 0.2, 0.2, 0.1], [2.848049683802928e-16, 2.783244677145701e-16, 2.8568720858194903e-16, 2.6409906663473615e-16, 2.704611318747015e-16, 2.773188933027618e-16, 2.8107561984968204e-16, 2.6924151168487263e-16, 2.748389085112028e-16, 2.748389085112028e-16, 2.711581925031757e-16], 0.018321910982848792)

CartesianDiscreteModel()

3.4410560032879544e-16

Test Passed

AnalyticalMapCubedSphereDiscreteModel()

0.0008252644349794302

Test Broken
  Expression: rel_error < tol

CartesianDiscreteModel()

3.261094076995979e-16

Test Passed

AnalyticalMapCubedSphereDiscreteModel()

3.861507338905354e-16

Test Passed

CartesianDiscreteModel()

2.878278973033339e-16

Test Passed

AnalyticalMapCubedSphereDiscreteModel()

3.2341303102026523e-16

Test Passed

julia> 

[codecov-commenter]

Codecov Report

Merging #12 (7ec0a99) into master (3091346) will not change coverage. The diff coverage is n/a.

@@           Coverage Diff           @@
##           master      #12   +/-   ##
=======================================
  Coverage   62.72%   62.72%           
=======================================
  Files          11       11           
  Lines         440      440           
=======================================
  Hits          276      276           
  Misses        164      164           

---

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 3091346...7ec0a99. Read the comment docs.

[santiagobadia]

Question 1

Let us start with question 1. I would say that the approximation you make is

div(Ju) = div(J)u + J div(u)

It seems that

d_i J_ij = 0

Now, let us take a look at your cubed sphere mapping. The map is orthogonal, it preserves inner products.

It transforms the dX1, dX2 1-forms into dx1, dx2, where dx1 \perp dx2. There is a scaling factor, but again it is constant. The map is homogeneous, we do the same for all points, transform segments into geodesics.

j = [dx1, dx2] expressed in R3 Euclidean space.

Now, d_x1(dx1) = 0 (note that d_x1 is a covariant derivative). As we move along the geodesic, the only change is in the normal component, no derivative on the tangent space.

As a result, I think that you can prove that the term we are assuming to be zero is in fact zero for this very specific map.

I wonder, what happens for other maps, like a map into HEX that are not affine. Do we still have the same value for div and DIV? I would say no. You can even check in 2D using a non-affine HEX mesh.

Question 3

They are equivalent in 2D, it is obvious, I would say. div(RT) subset DG0, so by the def'on of the projection, they must be equal. In any case, which order are you using? For DG0, be careful. u^2 in DG0 for u in DG0.

I am not sure why they are not the same for the sphere. Errors related to the geometrical map? I have to think about it.

Question 5

Idem question 1

[amartinhuertas]

OK, I have added some tests in order to shed some more light into some of your questions. Let me answer them separately.

Question 1 I wonder, what happens for other maps, like a map into HEX that are not affine. Do we still have the same value for div and DIV? I would say no. You can even check in 2D using a non-affine HEX mesh.

I have tried with a distorted mesh, see: here and here. I visualized the mesh in paraview to double check that the mapping is not affine (see below). The results are available here:

https://github.com/gridapapps/GridapGeosciences.jl/pull/12/checks?check_run_id=3469551972#step:6:328

for different mixed RT-DG FE pair orders. Assuming that these results are correct, it seems that we still have the same values for div and DIV on a 2D mesh with a non-affine map. How can this be? Any ideas? (I am still digesting your reasoning to question 1). Perhaps we should start looking at the details of how the divergence is actually being computed under the hood. Might it be end being computed using AD by chance? (no idea, just a crazy hypothesis).

Question 3 They are equivalent in 2D, it is obvious, I would say. div(RT) subset DG0, so by the def'on of the projection, they must be equal. In any case, which order are you using? For DG0, be careful. u^2 in DG0 for u in DG0.

First, here "For DG0, be careful. u^2 in DG0 for u in DG0." the last DG0 is RT0, right? Why u^2 in DG0? if I take a generic polynomial vector-valued function in Q1,0xQ0,1, and take the dot product by itself, I get a polynomial with monomials in {y^2 , y , x^2, x, 1}.

I was solely using order 0. You can see the results for different orders here:

https://github.com/gridapapps/GridapGeosciences.jl/pull/12/checks?check_run_id=3469551972#step:6:344

It seems that the larger the order, the larger the relative error (???), although it keeps being quite small even for order 5. Perhaps a conditioning issue of the mass matrix? (no idea, crazy hyphotesis)

I am not sure why they are not the same for the sphere. Errors related to the geometrical map? I have to think about it.

I have also run this case for increasing orders here:

https://github.com/gridapapps/GridapGeosciences.jl/pull/12/checks?check_run_id=3469551972#step:6:350

The higher the order, the higher the error.

Question 5 Idem question 1

Why? For me, it is not obvious that q1 implies q5 (still thinking on why q1 implies q5).

[santiagobadia]

With the proof I have sent you @amartinhuertas all points are solved now but 4.

4 seems to be due to the fact that div(v) is in the space of the L2 projection for u.^2

When you have a non-affine map, the div does not belong to this space anymore and thus this is not the same.

[davelee2804]

Hi @amartinhuertas , @santiagobadia , So by points (3-5) ∫(DIV(v)*(u⋅u))d\omega and ∫(DIV(v)*(u_dot_u_projection))d\omega are also NOT equivalent on a non-affine geometry!! I'm guessing this is due to inexact integration?? So I suppose this re-enforces the importance of projecting the potential (\phi = u.u/2 + gh) onto L2 before taking the gradient. This is probably the right thing to do from a Hamiltonian perspective, since \phi is the variational derivative of the energy w.r.t. h, and we want to compute this uniquely (or as close to as possible) to balance out the corresponding energy exchange from the divergence of the mass flux (which we are also projecting)

[davelee2804]

Why u^2 in DG0?

Interesting question @amartinhuertas . From a fluids perspective, u.u is a potential (it is associated with divergent motions, not rotation), and as such belongs in the same space as the other potential terms (the DG space), as opposed to the rotational terms, which are in the CG space. Also there is no assumption of continuity on u.u.

HOWEVER, as you point out, u.u as a higher polynomial degree, so by projecting into the DG space we loose information right?

I don't have an answer to this, except that in differential forms the equivalent expression is u^(u), so if u is a k-form (in a space of dimension n), u (the Hodge-* of u) is a n-k form, so the wedge product of these two things is an n-form (ie: a volume form, the DG space).

[davelee2804]

...Note however that the choice to put the potential terms in the DG space and the rotational terms in the CG space is somewhat arbitrary. In practice we could put the rotational terms in the DG space, the potential terms in the CG space (and the velocity in H(curl) :) ). Then the u.u term would be (u^(u)), a 0-form (in the CG space)

[amartinhuertas]

With the proof I have sent you @amartinhuertas all points are solved now but 4. 4 seems to be due to the fact that div(v) is in the space of the L2 projection for u.^2 When you have a non-affine map, the div does not belong to this space anymore and thus this is not the same.

@santiagobadia ... here do you actually mean 3 instead of 4?

[amartinhuertas]

This PR can be already closed. The main conclusion (proven mathematically) is that ∫(∇⋅(v)*(x))dΩ is equivalent to ∫(DIV(v)*(x))dω no matter what x actually is, e.g.,

• p in P \subset L2(Omega)
• g*h + 0.5*u⋅u, with h in P \subset L2(Omega) and u in U \subset H(div,Omega)

On the other hand, on a non-affine map, ∫(DIV(v)*(u⋅u))dω, with u,v in U \subset H(div,Omega) is NOT equivalent to ∫(DIV(v)*(uu_L2))dω, with uu_L2 being the L2 projection of u⋅u into P \subset L2(Omega). The latter should be this avoided in order to introduce unnecessary extra error.