Interface triangulations and CellFields

[JanSuntajs]

circle_n_2_size_2.0_80.0_x0_0.0_0.0_0.0_end_0.msh.zip

Hello,

I am trying out different ways to simulate diffusion in a heterogenous medium using Gridap. When I try to model the problem using separate (independent) FESpaces on two different domains, I have observed some strange behavior, possibly related to how CellFields behave on interface triangulations.

Unfortunately, as my project contains some glueing code with convenience routines (particulary for easier meshing and interfacing with solvers from DifferentialEquations.jl), I cannot produce a MWE quickly, but I am appending the relevant code bits here.

I believe the crux of the problem lies here:

# set the diffusion constants (trivial test case, everything should be homogenous here).
# For simplicity, I am modelling only two domains.
D1 = 1.
D2 = 1.

# I am appending the mesh used at the top of this message. Physical domains are denoted as "face<i>" in the
# .msh mesh, which can then be converted into domain tags 1 and 2 inside Gridap.
tags_field = CellField(domain_tags, Ω)
D_from_tag(tag) = tag == 1 ? D1 : D2    # very simple mapping function in this case
D_field = D_from_tag ∘ tags_field

# some more code ...
# The heart of the problem: mean and jump operators
@inline _jump(w1, w2) = w1.⁺ -  w2.⁻
@inline _mean(w1, w2) = w1.⁺ + w2.⁻

# this fails - namely, the calculated concentration field is discontinous
# even though this is used in terms used to enforce continuity
@inline _mean(u1, u2, n) = (D_field.⁺ ⋅ ∇(u1.⁺) ⋅ n.⁺ +  D_field.⁻ ⋅ ∇(u2.⁻) ⋅ n.⁻)
# strangely enough, this works as expected:
@inline _mean(u1, u2, n) = (D1 ⋅ ∇(u1.⁺) ⋅ n.⁺ +  D2 ⋅ ∇(u2.⁻) ⋅ n.⁻)

As implied above, if I define the mean operator on the interface as in the first case, I obtain discontinuities in the concentration field, which is not the case using the second approach. For clarity, I now also append the preparation of FESpaces, weak forms and relevant operators.

# preparation of FESpaces
order = 2
degree = 2*order

Ω = Triangulation(model)
Γ = BoundaryTriangulation(model, tags="Casing")
dΓ = Measure(Γ, degree)
n_Γ = get_normal_vector(Γ)
dΩ = Measure(Ω, degree)

Ω1 = Interior(model, tags="face1")
Ω2 = Interior(model, tags="face2")
dΩ1 = Measure(Ω1, degree)
dΩ2 = Measure(Ω2, degree)

Γ_if = InterfaceTriangulation(Ω1, Ω2)
d_Γ_ifs = Measure(Γ_if, degree)
n_Γ_ifs = get_normal_vector(Γ_if)

reffe = ReferenceFE(lagrangian, Float64, order)
V1 = TestFESpace(Ω1, reffe, conformity=:H1)
V2 = TestFESpace(Ω2, reffe, conformity=:H1)
U1 = TransientTrialFESpace(V1)
U2 = TransientTrialFESpace(V2)

V = MultiFieldFESpace([V1, V2])
U = TransientMultiFieldFESpace([U1, U2])


D1 = 1.
D2 = 1.

C1 = 1.
C2 = 1.

flux(t) = t <= 1 ? sin(π * t) : 0.

# Diffusion constants

tags_field = CellField(domain_tags, Ω)
D_from_tag(tag) = tag == 1 ? D1 : D2    # very simple mapping function in this case
D_field = D_from_tag ∘ tags_field

# weak form

function _prefactors(triangulation, measure, param)

    dim = num_cell_dims(triangulation)
    _hmat = get_array(∫(1)measure)
    _α = CellField(lazy_map(h -> param/h.^dim, _hmat), triangulation)

    return _α
end
_α = _prefactors(Γ_if, d_Γ_ifs, 0.1)

# Note that, since we no longer have a single global FE space, we need to provide
# left and right functions to define jumps and averages at the interface.
@inline _jump(w1, w2) = w1.⁺ -  w2.⁻
@inline _mean(w1, w2) = w1.⁺ + w2.⁻

# this works
@inline _mean(u1, u2, n) = (D1 ⋅ ∇(u1.⁺) ⋅ n.⁺ +  D2 ⋅ ∇(u2.⁻) ⋅ n.⁻)
# this doesn't - a discontinuity in concentration is produced; 
# @inline _mean(u1, u2, n) = (D_field.⁺ ⋅ ∇(u1.⁺) ⋅ n.⁺ +  D_field.⁻ ⋅ ∇(u2.⁻) ⋅ n.⁻)
# note that in regular volume and surface integrations, I can work with CellFields just fine. 

# for some reason, this works better with no symmetrization term in the Nitsche scheme
res(t, (u1, u2), (v1, v2)) = begin

    ∫(∂t(u1) ⋅ v1)dΩ1 + ∫(∂t(u2) ⋅ v2)dΩ2 +
    ∫(-flux(t) ⋅ v1)dΓ + # only on one side this time
    ∫((D_field ⋅ ∇(u1)) ⋅ ∇(v1))dΩ1 + ∫((D_field ⋅  ∇(u2)) ⋅ ∇(v2))dΩ2 - 
    ∫(_jump(v1, v2) ⋅ _mean(u1, u2, n_Γ_ifs))d_Γ_ifs +
    ∫(_α ⋅ _jump(u1, u2) ⋅ _jump(v1, v2))d_Γ_ifs 

end

jac(t, (u1, u2), (du1, du2), (v1, v2)) = begin

    ∫((D_field ⋅ ∇(du1)) ⋅ ∇(v1))dΩ1 + ∫((D_field ⋅ ∇(du2)) ⋅ ∇(v2))dΩ2 -
    ∫(_jump(v1, v2) ⋅ _mean(du1, du2, n_Γ_ifs))d_Γ_ifs +
    ∫(_α ⋅ _jump(du1, du2) ⋅ _jump(v1, v2))d_Γ_ifs 

end

jac_t(t, (u1, u2), (dtu1, dtu2), (v1, v2)) = begin

    ∫(dtu1 ⋅ v1)dΩ1 + ∫(dtu2 ⋅ v2)dΩ2

end

# Operators

init_guess = interpolate_everywhere([c1, c2], U(0.))
op = TransientFEOperator(res, jac, jac_t, U, V)

# As said above, I then use some convenience routines to solve this via DifferentialEquations.jl routines 

Can someone reproduce this behavior or would be willing to look into it? I would be happy to cooperate on this issue. Interestingly, this same issue does not appear if I use global FESpaces U and V (hence if I do not break the problem into subdomains), there the CellFields work nicely on the interfaces. The reason I am bringing this up is that I would eventually want to treat more complex diffusion cases with possible discontinuities in certain quantities (think of generalized diffusion with continous chemical potentials across interfaces and consequently discontinous concentration fields).

Also, the issue with interfaces and AD still persists (it probably boils down to some missing signature) - namely, one needs to manually add the Jacobians if working with interfaces in transient problems, as otherwise AD will complain about some function belonging to an interface definition.

Best,
Jan

[JordiManyer]

I think your signs are wrong in the mean and jump functions. The normals are not absolute, they are relative to the cell they are attached to, so n.⁺ and n.⁻ have opposite signs. Therefore on a SkeletonTriangulation (such as the one one returned by InterfaceTriangulation(Ω1, Ω2)) the jump and mean of quantities containing the normals go with a + and - respectively (not the other way around).

[JanSuntajs]

Hello,
thanks for your reply, this is indeed a good catch! Unfortunately, and perhaps rather surprisingly, changing the sign from plus to minus in the definition of the mean doesn't seem to have much effect on the concentration field. I am judging this from the comparison of the L2 error norm of the difference between the solution of the above problem and the solution of the homogenous diffusion problem (with $D1 = D2$).

Also, for some reason, if I define the mean operators via CellFields rather than with scalar diffusion coefficients, the solver performance (I am typically casting the equations in mass matrix form and using FBDF or QNDF solvers) worsens by orders of magnitude (a quick estimate would be from approx. 10 seconds to tens of minutes). Peculiarly, this only affects the interface terms, in the bulk and surface terms everything works fine.

[JanSuntajs]

@JordiManyer After some further study, I think this pretty much boils down to Issue #982. As I've noticed, multiplying the negative side quantities by a scalar indeed evaluates to zero, which is also why no difference was observed while changing the sign.

Is there a simple way around this or an easy to implement fix?

Best,
Jan

[JordiManyer]

I am going to give a detailed explanation that you don't need to fully understand, I just want to leave it on record:

Consider the following MWE:

using Gridap
using Gridap.Geometry

model = CartesianDiscreteModel((0,1,0,1),(2,2))

Ω = Triangulation(model)
Ω1 = Triangulation(model,[1,2])
Ω2 = Triangulation(model,[3,4])
Λ = InterfaceTriangulation(Ω1,Ω2)

reffe = ReferenceFE(lagrangian,Float64,1)
V1 = FESpace(Ω1,reffe)
V2 = FESpace(Ω2,reffe)

u1 = interpolate(1.0,V1)
u2 = interpolate(1.0,V2)

dΛ = Measure(Λ,2)
f(u) = sum(∫(u)dΛ)

f(u1.plus)             # OK
f(u2.minus)         # OK
f(1.0*u1.plus)     # OK
f(1.0*u2.minus) # BAD

o = CellField(1.0,Ω)
o1 = CellField(1.0,Ω1)
o2 = CellField(1.0,Ω2)
f(o*u2.minus)     # OK
f(o1*u1.plus)      # OK
f(o2*u2.minus)  # BAD

o2bis = CellField(fill(1.0,num_cells(Ω2)),Ω2)
f(o2bis*u2.minus) # BREAKS, which is expected

# FIX A: 
o2 = CellField(1.0,Ω2)
f(o2.minus*u2.minus)  # OK

# FIX B: 
o = CellField(1.0,Ω)
f(o*u2.minus)  # OK

# FIX C: 
dΛm = Measure(Λ.minus,2)
sum(∫(1.0*u2.minus)dΛm)

The underlying issue is that the constant gets mapped to a constant CellField defined on the triangulation where the operand (here u2) is defined. Therefore you end up with the equivalent of o2. The underlying array for this CellField is a Fill, since it is a constant field over the triangulation.
Then: when performing the change of domain, o2 gets passed onto this function which by default takes the .plus side. In general, if the constant CellField is defined over both sides of the interface the comment in the function is correct and we can take the plus side without any issues. However in this particular case things do not work because our constant is not defined on the plus side of the interface.

The correct way of doing things, in my view, is to

• FIX A: Ensure the correct side is taken during the intergrals
• FIX B: Ensure the constant is defined everywhere
• FIX C: define the integral on a single side of the interface

Overall, I would say this is a case of Gridap trying to be too flexible. We should probably not allow these ambiguous cases to happen, and just enforce that integrals over SkeletonTriangulations should ONLY involve a) CellField objects defined on the skeleton itself and b) CellFieldAt objects. Otherwise, we leave the door open for errors like these which are hard to understand for users that do not know how things happen under the hood.

[JanSuntajs]

Thank you for your effort and a comprehensive response!

Taking all this into account, I was able to simulate a generalized example of heterogenous diffusion, in which
the gradients of the chemical potential drive the diffusion and the concentration field can have jumps at the interfaces (shown are two snapshots at two different times):

[JordiManyer]

Glad to hear! If you think it's different enough from what we already have and you have time, consider making a tutorial!

Cheers

[JanSuntajs]

Thanks for your prompt responses and help! Will surely consider it once I get some pending projects done. Though the solutions were obtained using OrdinaryDiffEq.jl solvers, which on its own could be a subject for a tutorial.

Best,
Jan