Problem with numerical fluxes between elements interface (dS). How to define two different functions for two different subdomains using (dS)

[tfkonopka]

Hello everyone,

I need some help to define a function in two different subdomains. The function F has four parameters which two of them are different in each subdomain. For example, at the subdomain dx(0), no = 1 nw =1. At the subdomain dx(1), no = 2 nw = 2. For the integrals across the subdomain, dx(0) and dx(1), it is easy.

However, I don’t know how to deal with ds and dS . For example, at the stabilization term bellow we have (F(s, mu, no, nw)("+") - un("-") * F(s, mu, no, nw)("-") ) * dS. So, how can I calculate two different F function (no and nw are different according the subdomain) between the elements interface at dx(0) and dx(1) or as in a shared boundary like the next figure.

In this code s, un, un_h, u and r are trial or test function; mu = 1 sbar = 1 are constants.

Any suggestions would be very helpful.

Summarized code:

un = 0.5 * (inner(u0, n) + sqrt(inner(u0, n) * inner(u0, n)))
un_h = 0.5 * (inner(u0, n) - sqrt(inner(u0, n) * inner(u0, n)))

def F(s, mu_rel, no, nw):
    return s ** nw / (s ** nw + mu_rel * (1.0 - s) ** no)

    stabilisation = dt("+") * inner(
        jump(r), un("+") * F(s, mu, no, nw)("+") - un("-") * F(s, mu, no, nw)("-")
    ) * dS + dt * r * un_h * sbar * ds(1)

    L3 = (
        phi * r * (s - s0) * dx(0)
        + r * (s - s0) * dx(1)
        - dt * inner(grad(r), F(s, mu, no, nw) * u) * dx(0)
        - dt * inner(grad(r), F(s, mu, no, nw) * u) * dx(1)
        + dt * r * F(s, mu, no, nw) * un * ds
        + stabilisation
    )

    # Total L
    L = L1 + L2 + L3

Entire code:

# Total mobility
def lmbdainv(s, mu_rel):
    return 1.0 / ((1.0 / mu_rel) * s ** 2 + (1.0 - s) ** 2)


# Fractional flow function
def F(s, mu_rel, no, nw):
    return s ** nw / (s ** no + mu_rel * (1.0 - s) ** nw)


class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[1], (0.3, 0.7)) and between(x[0], (0.2, 0.6))


def tensor_jump(v, n):
    return ufl.outer(v, n)("+") + ufl.outer(v, n)("-")    

Kinv = Constant(1 / (k_matriz))    
dt = Constant(0.5)
mu = 1
mu_b = 1
phi = Constant(0.4)
t = 0.0
T = 50 * float(dt)

mesh = UnitSquareMesh(Nx, Ny, "crossed")
n = FacetNormal(mesh)
order = 1
BDM = FiniteElement("BDM", mesh.ufl_cell(), order)
DG = FiniteElement("DG", mesh.ufl_cell(), order - 1)
Element = MixedElement([BDM, DG, DG])
mixed_space = FunctionSpace(mesh, Element)  

# Function spaces and functions
V = TestFunction(mixed_space)
dU = TrialFunction(mixed_space)
U = Function(mixed_space)
U0 = Function(mixed_space)

v, q, r = split(V)
u, p, s = split(U)
u0, p0, s0 = split(U0)

s_mid = 0.5 * (s0 + s)

obstacle = Obstacle()
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
obstacle.mark(domains, 1)

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
left = AutoSubDomain(lambda x: near(x[0], 0.0))
right = AutoSubDomain(lambda x: near(x[0], 1.0))
bottom = AutoSubDomain(lambda x: near(x[1], 0.0))
top = AutoSubDomain(lambda x: near(x[1], 1.0))

# Define boundary markers
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dx = Measure("dx", domain=mesh, subdomain_data=domains)

alpha = 35
h = CellDiameter(mesh)
h2 = ufl.Min(h("+"), h("-"))

stab = (
	mu * (alpha / h2) * inner(tensor_jump(u, n), tensor_jump(v, n)) * dS
	- mu * inner(avg(grad(u)), tensor_jump(v, n)) * dS
	- mu * inner(avg(grad(v)), tensor_jump(u, n)) * dS
)

L1 = (
	mu_b * inner(grad(u), grad(v)) * dx(1)
	+ inner(v, lmbdainv(s_mid, mu) * Kinv * u) * dx(0)
	- div(v) * p * dx(1)
	- div(v) * p * dx(0)
	# + pin * dot(v, n) * ds(1)
	+ pout * dot(v, n) * ds(3)
	+ stab
)

L2 = q * div(u) * dx(1) + q * div(u) * dx(0)
un = 0.5 * (inner(u0, n) + sqrt(inner(u0, n) * inner(u0, n)))
un_h = 0.5 * (inner(u0, n) - sqrt(inner(u0, n) * inner(u0, n)))	

stabilisation = dt("+") * inner(
	jump(r), un("+") * F(s_mid, mu, no, nw)("+") - un("-") * F(s_mid, mu, no, nw)("-")
) * dS + dt * r * un_h * sbar * ds(1)

L3 = (
	phi * r * (s - s0) * dx(0)
	+ r * (s - s0) * dx(1)
	- dt * inner(grad(r), F(s_mid, mu, no, nw) * u) * dx(0)
	- dt * inner(grad(r), F(s_mid, mu, no, nw) * u) * dx(1)
	+ dt * r * F(s_mid, mu, no, nw) * un * ds
	+ stabilisation
)

# Total L
L = L1 + L2 + L3

# Jacobian
J = derivative(L, U, dU)
problem = NonlinearVariationalProblem(L, U, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)

while t < T:
	t += float(dt)
	U0.assign(U)
	solver.solve()
	u, p, s = U.split()   

[conpierce8]

You could define the piecewise-constant parameters no and nw with a UserExpression by implementing the eval_cell method, as in the simplified example below.

from fenics import *

class PiecewiseConstant(UserExpression):
    def __init__(self, values, markers, **kwargs):
        self._values = values
        self._markers = markers
        super().__init__(**kwargs)
    
    def eval_cell(self, values, x, cell):
        values[0] = self._values[self._markers[cell.index]]
    
    def value_shape(self):
        return tuple()

class Subdomains(SubDomain):
    def __init__(self, example_no, **kwargs):
        self._ex_no = example_no
        super().__init__(**kwargs)
    def inside(self, x, on_boundary):
        if self._ex_no == 1:
            return between(x[1], (0.3, 0.7)) and between(x[0], (0.2, 0.6))
        else:
            return between(x[1], (0.3, 0.7)) and between(x[0], (0.6, 1.0))

class SubdomainsBoundary(SubDomain):
    def __init__(self, example_no, **kwargs):
        self._ex_no = example_no
        super().__init__(**kwargs)
    def inside(self, x, on_boundary):
        if self._ex_no == 1:
            return (
                near(abs(x[0] - 0.4), 0.2) and abs(x[1] - 0.5) < 0.2 + DOLFIN_EPS
                or
                near(abs(x[1] - 0.5), 0.2) and abs(x[0] - 0.4) < 0.2 + DOLFIN_EPS
            )
        else:
            return (
                near(x[0], 0.6) and abs(x[1] - 0.5) < 0.2 + DOLFIN_EPS
                or
                near(abs(x[1] - 0.5), 0.2) and abs(x[0] - 0.8) < 0.2 + DOLFIN_EPS
            )

# Example usage
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "DG", 0)
s = interpolate(Expression("1", degree=1), V)

no = {1:1.0, 2:2.0}  #       { 1.0   ; x \in dx(1)
                     #  no = {
                     #       { 2.0   ; x \in dx(2)

def F(s, no):
    return s * no

all_F = []
for i in [1, 2]:
    mf_cells = MeshFunction("size_t", mesh, mesh.topology().dim(), 1)
    sub  = Subdomains(i)
    sub.mark(mf_cells, 2)
    mf_facets = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
    
    subBnd = SubdomainsBoundary(i)
    subBnd.mark(mf_facets, 1)
    dS = Measure('dS', mesh, subdomain_data=mf_facets)

    pc = PiecewiseConstant( no, mf_cells )
    Fi = project(F(s, pc), V)
    Fi.rename(f"F{i}", f"F{i}")
    all_F.append(Fi)
    
    print(f"avg(F{i})*dS(1) =", assemble(avg(Fi)*dS(1)))

with XDMFFile("example.xdmf") as xdmf:
    xdmf.parameters.update(
        {
            "functions_share_mesh": True,
            "rewrite_function_mesh": False
        }
    )
    for Fi in all_F:
        xdmf.write(Fi, 0)

producing

avg(F1)*dS(1) = 2.3999999999999995
avg(F2)*dS(1) = 1.7999999999999996

[tfkonopka]

Dear @conpierce8

I wolud like to thank for the help.

I have being working on it since last week , but still not working properly.

The code bellow is the best answer that I could reach, but the simulation does not converge for all time steps and the the physics is not correct.

The space “s” represents the saturation on this problem and must be between 0 and 1. Furthermore, for bordering cells the downstream saturation in a cell must inferior to the upstream one. Which does not happen as we can see in the figures bellow

Any suggestions would be very helpful.

Best Regards

The modified code:

from fenics import *
import numpy as np
from mshr import *
import ufl


# Total mobility
def lmbdainv(s, mu_rel, no, nw):
    return 1.0 / ((1.0 / mu_rel) * s ** nw + (1.0 - s) ** no)


# Fractional flow function
def F(s, mu_rel, no, nw):
    return s ** nw / (s ** nw + mu_rel * (1.0 - s) ** no)


class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[1], (0.3, 0.7)) and between(x[0], (0.2, 0.6))


def tensor_jump(v, n):
    return ufl.outer(v, n)("+") + ufl.outer(v, n)("-")


class PiecewiseConstant(UserExpression):
    def __init__(self, values, markers, **kwargs):
        self._values = values
        self._markers = markers
        super().__init__(**kwargs)

    def eval_cell(self, values, x, cell):
        values[0] = self._values[self._markers[cell.index]]

    def value_shape(self):
        return tuple()


class Subdomains(SubDomain):
    def __init__(self, example_no, **kwargs):
        self._ex_no = example_no
        super().__init__(**kwargs)

    def inside(self, x, on_boundary):
        if self._ex_no == 1:
            return between(x[1], (0.3, 0.7)) and between(x[0], (0.2, 0.6))
        else:
            return between(x[1], (0.3, 0.7)) and between(x[0], (0.6, 1.0))


class SubdomainsBoundary(SubDomain):
    def __init__(self, example_no, **kwargs):
        self._ex_no = example_no
        super().__init__(**kwargs)

    def inside(self, x, on_boundary):
        if self._ex_no == 1:
            return (
                near(abs(x[0] - 0.4), 0.2)
                and abs(x[1] - 0.5) < 0.2 + DOLFIN_EPS
                or near(abs(x[1] - 0.5), 0.2)
                and abs(x[0] - 0.4) < 0.2 + DOLFIN_EPS
            )
        else:
            return (
                near(x[0], 0.6)
                and abs(x[1] - 0.5) < 0.2 + DOLFIN_EPS
                or near(abs(x[1] - 0.5), 0.2)
                and abs(x[0] - 0.8) < 0.2 + DOLFIN_EPS
            )


Kinv = Constant(1 / (1e-6))
# Time step
dt = Constant(0.5)
mu = 1
mu_b = 1
phi = Constant(0.4)
pin = 6000
pout = 1000
sbar = 1

no_0 = 2
nw_0 = 2
no_1 = 1
nw_1 = 1

t = 0.0
T = 300 * float(dt)

mesh = UnitSquareMesh(20, 20, "crossed")
n = FacetNormal(mesh)

order = 1

BDM = FiniteElement("BDM", mesh.ufl_cell(), order)
DG = FiniteElement("DG", mesh.ufl_cell(), order - 1)
Element = MixedElement([BDM, DG, DG])
mixed_space = FunctionSpace(mesh, Element)

# Function spaces and functions
V = TestFunction(mixed_space)
dU = TrialFunction(mixed_space)
U = Function(mixed_space)
U0 = Function(mixed_space)

v, q, r = split(V)
u, p, s = split(U)
u0, p0, s0 = split(U0)

s_mid = 0.5 * (s0 + s)

obstacle = Obstacle()
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
obstacle.mark(domains, 1)

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
left = AutoSubDomain(lambda x: near(x[0], 0.0))
right = AutoSubDomain(lambda x: near(x[0], 1.0))
bottom = AutoSubDomain(lambda x: near(x[1], 0.0))
top = AutoSubDomain(lambda x: near(x[1], 1.0))

# Define boundary markers
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dx = Measure("dx", domain=mesh, subdomain_data=domains)

bc1 = DirichletBC(mixed_space.sub(0), Constant((6.0e-3, 0.0)), boundaries, 1)
bc2 = DirichletBC(mixed_space.sub(0), Constant((0.0, 0.0)), boundaries, 2)
# # bc3 = DirichletBC(VQ.sub(0), Constant((0.0, 0.0)), boundaries, 3)
bc4 = DirichletBC(mixed_space.sub(0), Constant((0.0, 0.0)), boundaries, 4)

bcs = [bc1, bc2, bc4]  # velocity BC

alpha = 35
h = CellDiameter(mesh)
h2 = ufl.Min(h("+"), h("-"))

# ------------------------------ Added code -------------------------------------
mf_cells = MeshFunction("size_t", mesh, mesh.topology().dim(), 1)
mf_facets = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
dS = Measure("dS", mesh, subdomain_data=mf_facets)
# --------------------------------------------------------------------------------------
stab = (
    mu * (alpha / h2) * inner(tensor_jump(u, n), tensor_jump(v, n)) * dS
    - mu * inner(avg(grad(u)), tensor_jump(v, n)) * dS
    - mu * inner(avg(grad(v)), tensor_jump(u, n)) * dS
)

L1 = (
    mu_b * inner(grad(u), grad(v)) * dx(1)
    + inner(v, lmbdainv(s_mid, mu, no_0, nw_0) * Kinv * u) * dx(0)
    - div(v) * p * dx(1)
    - div(v) * p * dx(0)
    # + pin * dot(v, n) * ds(1)
    + pout * dot(v, n) * ds(3)
    + stab
)

L2 = q * div(u) * dx(1) + q * div(u) * dx(0)
un = 0.5 * (inner(u0, n) + sqrt(inner(u0, n) * inner(u0, n)))
un_h = 0.5 * (inner(u0, n) - sqrt(inner(u0, n) * inner(u0, n)))


# ------------------------------ Added code -------------------------------------
no = {1: 1.0, 2: 2.0}

VVV = FunctionSpace(mesh, "DG", 0)

for i in [1]:
    sub = Subdomains(i)
    sub.mark(mf_cells, 2)
    subBnd = SubdomainsBoundary(i)
    subBnd.mark(mf_facets, 1)
    dS = Measure("dS", mesh, subdomain_data=mf_facets)
    noo = PiecewiseConstant(no, mf_cells)
    noo_proj = project(noo, VVV)


stabilisation = (
    dt("+")
    * inner(
        jump(r),
        un("+") * F(s_mid, mu, noo_proj, noo_proj)("+")
        - un("-") * F(s_mid, mu, noo_proj, noo_proj)("-"),
    )
    * dS
    + dt * r * un_h * sbar * ds(1)
)

# -----------------------------------------------------------------------

L3 = (
    phi * r * (s - s0) * dx(0)
    + r * (s - s0) * dx(1)
    - dt * inner(grad(r), F(s_mid, mu, no_0, nw_0) * u) * dx(0)
    - dt * inner(grad(r), F(s_mid, mu, no_1, nw_1) * u) * dx(1)
    + dt * r * F(s_mid, mu, no_0, nw_0) * un * ds
    + stabilisation
)

# Total L
L = L1 + L2 + L3

J = derivative(L, U, dU)

problem = NonlinearVariationalProblem(L, U, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)

s_file = File("teste/saturation.pvd")

while t < T:
    t += float(dt)
    U0.assign(U)
    solver.solve()
    u, p, s = U.split()
    s_file << s

[conpierce8]

Can you provide the equations you are trying to implement or a reference to the problem you want to implement? I don’t recognize the physics you are solving so I can’t provide much help. (A further disclaimer: I don’t know much about stabilization techniques or convergence, so I probably can’t provide much help anyway.)

I did notice the following:

1. Your tensor_jump function is implemented as twice the average of u \otimes n rather than the jump (difference). Changing this:

   @tfkonopka

   def tensor_jump(v, n):
       return ufl.outer(v, n)("+") + ufl.outer(v, n)("-")
   

   to this:

   @tfkonopka

   def tensor_jump(v, n):
       return jump(ufl.outer(v, n))
   

   improves convergence for me.

2. The combination of the following code snippets

   @tfkonopka

   obstacle = Obstacle()
   domains = MeshFunction("size_t", mesh, mesh.topology().dim())
   domains.set_all(0)
   obstacle.mark(domains, 1)
   

   @tfkonopka

   L1 = (
       mu_b * inner(grad(u), grad(v)) * dx(1)
       + inner(v, lmbdainv(s_mid, mu, no_0, nw_0) * Kinv * u) * dx(0)
   

   @tfkonopka

   L3 = (
       phi * r * (s - s0) * dx(0)
       + r * (s - s0) * dx(1)
       - dt * inner(grad(r), F(s_mid, mu, no_0, nw_0) * u) * dx(0)
       - dt * inner(grad(r), F(s_mid, mu, no_1, nw_1) * u) * dx(1)
   

   implies that no_0 and nw_0 apply to the exterior region, while no_1 and nw_1 apply to the interior region [0.2,0.6]\times [0.3, 0.7]. However, in your usage of the code I supplied, i.e.

   @tfkonopka

   mf_cells = MeshFunction("size_t", mesh, mesh.topology().dim(), 1)
   

   @tfkonopka

   no = {1: 1.0, 2: 2.0}
   
   VVV = FunctionSpace(mesh, "DG", 0)
   
   for i in [1]:
       sub = Subdomains(i)
       sub.mark(mf_cells, 2)
       subBnd = SubdomainsBoundary(i)
       subBnd.mark(mf_facets, 1)
       dS = Measure("dS", mesh, subdomain_data=mf_facets)
       noo = PiecewiseConstant(no, mf_cells)
   

   you have used the opposite parameters. Here, the inner domain is marked 2 and the surroundings are marked 1. The line no = {1: 1.0, 2: 2.0} specifies that domains marked 1 should have no = 1.0, and the domains marked 2 should have no = 2.0. This is the opposite of how it is implemented in L1 and L2.

3. The definition of mf_cells is unnecessary, because it contains the same information as domains. Likewise, the definition of mf_facets is unnecessary because you don’t use any restrictions on dS (e.g. dS(1)).

4. Several terms in the weak form can be simplified. For example,

   @tfkonopka

       - div(v) * p * dx(1)
       - div(v) * p * dx(0)
   

   can be written - div(v) * p * dx,

   @tfkonopka

   L2 = q * div(u) * dx(1) + q * div(u) * dx(0)
   

   can be written q * div(u) * dx,

   @tfkonopka

           un("+") * F(s_mid, mu, noo_proj, noo_proj)("+")
           - un("-") * F(s_mid, mu, noo_proj, noo_proj)("-")
   

   can be written jump(un * F(s_mid, mu, noo_proj, noo_proj)), and

   @tfkonopka

       - dt * inner(grad(r), F(s_mid, mu, no_0, nw_0) * u) * dx(0)
       - dt * inner(grad(r), F(s_mid, mu, no_1, nw_1) * u) * dx(1)
   

   can be combined into a single statement using the PiecewiseConstants:

       - dt * inner(grad(r), F(s_mid, mu, noo_proj, noo_proj) * u) * dx
   

   This doesn’t affect the solution but it makes the code a little cleaner.

5. This code

   @tfkonopka

   un = 0.5 * (inner(u0, n) + sqrt(inner(u0, n) * inner(u0, n))) 
   un_h = 0.5 * (inner(u0, n) - sqrt(inner(u0, n) * inner(u0, n)))
   

   appears to be a complicated way to write

   un = 0.5 * (inner(u0, n) + abs(inner(u0, n))) 
   un_h = 0.5 * (inner(u0, n) - abs(inner(u0, n)))
   

Implementing all these changes and slightly refining the mesh (30 x 30) gives me much better convergence:

from fenics import *
import numpy as np
# from mshr import *
import ufl


# Total mobility
def lmbdainv(s, mu_rel, no, nw):
    return 1.0 / ((1.0 / mu_rel) * s ** nw + (1.0 - s) ** no)


# Fractional flow function
def F(s, mu_rel, no, nw):
    return s ** nw / (s ** nw + mu_rel * (1.0 - s) ** no)


class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[1], (0.3, 0.7)) and between(x[0], (0.2, 0.6))


def tensor_jump(v, n):
    return jump(ufl.outer(v, n))


class PiecewiseConstant(UserExpression):
    def __init__(self, values, markers, **kwargs):
        self._values = values
        self._markers = markers
        super().__init__(**kwargs)

    def eval_cell(self, values, x, cell):
        values[0] = self._values[self._markers[cell.index]]

    def value_shape(self):
        return tuple()

Kinv = Constant(1 / (1e-6))
# Time step
dt = Constant(0.5)
mu = 1
mu_b = 1
phi = Constant(0.4)
pin = 6000
pout = 1000
sbar = 1

t = 0.0
T = 300 * float(dt)

mesh = UnitSquareMesh(30, 30, "crossed")
n = FacetNormal(mesh)

order = 1

BDM = FiniteElement("BDM", mesh.ufl_cell(), order)
DG = FiniteElement("DG", mesh.ufl_cell(), order - 1)
Element = MixedElement([BDM, DG, DG])
mixed_space = FunctionSpace(mesh, Element)

# Function spaces and functions
V = TestFunction(mixed_space)
dU = TrialFunction(mixed_space)
U = Function(mixed_space)
U0 = Function(mixed_space)

v, q, r = split(V)
u, p, s = split(U)
u0, p0, s0 = split(U0)

s_mid = 0.5 * (s0 + s)

# ============= DEFINITION OF SPATIALLY-VARYING PARAMETERS =====================
marker_inner = 1
marker_outer = 0

no_outer = 1
nw_outer = 1
no_inner = 2
nw_inner = 2

obstacle = Obstacle()
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(marker_outer)
obstacle.mark(domains, marker_inner)

no = {marker_inner: no_inner, marker_outer: no_outer}
nw = {marker_inner: nw_inner, marker_outer: nw_outer}

VVV = FunctionSpace(mesh, "DG", 0)

noo = PiecewiseConstant(no, domains)
noo_proj = project(noo, VVV)
nww = PiecewiseConstant(nw, domains)
nww_proj = project(nww, VVV)
# =========== END DEFINITION OF SPATIALLY-VARYING PARAMETERS ===================

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
left = AutoSubDomain(lambda x: near(x[0], 0.0))
right = AutoSubDomain(lambda x: near(x[0], 1.0))
bottom = AutoSubDomain(lambda x: near(x[1], 0.0))
top = AutoSubDomain(lambda x: near(x[1], 1.0))

# Define boundary markers
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dx = Measure("dx", domain=mesh, subdomain_data=domains)

bc1 = DirichletBC(mixed_space.sub(0), Constant((6.0e-3, 0.0)), boundaries, 1)
bc2 = DirichletBC(mixed_space.sub(0), Constant((0.0, 0.0)), boundaries, 2)
# # bc3 = DirichletBC(VQ.sub(0), Constant((0.0, 0.0)), boundaries, 3)
bc4 = DirichletBC(mixed_space.sub(0), Constant((0.0, 0.0)), boundaries, 4)

bcs = [bc1, bc2, bc4]  # velocity BC

alpha = 35
h = CellDiameter(mesh)
h2 = ufl.Min(h("+"), h("-"))

stab = (
    mu * (alpha / h2) * inner(tensor_jump(u, n), tensor_jump(v, n)) * dS
    - mu * inner(avg(grad(u)), tensor_jump(v, n)) * dS
    - mu * inner(avg(grad(v)), tensor_jump(u, n)) * dS
)

L1 = (
    mu_b * inner(grad(u), grad(v)) * dx(marker_inner)
    + inner(v, lmbdainv(s_mid, mu, no_outer, nw_outer) * Kinv * u) * dx(marker_outer)
    - div(v) * p * dx
    # + pin * dot(v, n) * ds(1)
    + pout * dot(v, n) * ds(3)
    + stab
)

L2 = q * div(u) * dx
un = 0.5 * (inner(u0, n) + abs(inner(u0, n)))
un_h = 0.5 * (inner(u0, n) - abs(inner(u0, n)))

stabilisation = (
    dt("+")
    * inner(
        jump(r),
        jump(un * F(s_mid, mu, noo_proj, nww_proj))
    )
    * dS
    + dt * r * un_h * sbar * ds(1)
)

# -----------------------------------------------------------------------

L3 = (
    phi * r * (s - s0) * dx(0)
    + r * (s - s0) * dx(1)
    - dt * inner(grad(r), F(s_mid, mu, noo_proj, nww_proj) * u) * dx
    + dt * r * F(s_mid, mu, no_outer, nw_outer) * un * ds
    + stabilisation
)

# Total L
L = L1 + L2 + L3

J = derivative(L, U, dU)

problem = NonlinearVariationalProblem(L, U, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)

s_file = File("teste/saturation.pvd")

i = 1
while t < T:
    print("Iteration", i)
    t += float(dt)
    U0.assign(U)
    solver.solve()
    u, p, s = U.split()
    s_file << s
    i += 1

[tfkonopka]

Dear @conpierce8

During this week I realized that the physics is correct. So, your solution for marking the domains work perfectly. Thank you

Best regardings
Thiago Konopka