Hi all,
I am solving a diffusion equation in a linear DG0 domain split into two: in each of the subdomains the diffusion coefficent assumes a different value. I want to implement a jump on an internal boundary which is modulated by a partitioning coefficient P and scaled by a parameter k. Thus, the flux at the internal bondary is given by -k(u('-') -P*u('+')). I am starting implementing it without considering the parameter k: (u('-') -P*u('+')). My weak form and the boundary condition are as follow:
F = (u_n-ui)*vi*dx - D_*dt/d_avg*dot(jump(ui,n),jump(vi,n))*dS
F += -P*jump(u_n)*jump(vi)*dS(120)
Yet, the jump across the boundary grows at every time step and never meets the partitioning coefficient P. Is there something wrong in my understanding of the jump operator? Any help/suggestion would be very helpful.
Here is my complete code, inclusive of mesh:
from re import X
import numpy as np
import pandas as pd
import sys
import pyvista
import matplotlib.pyplot as plt
import gmsh
import dolfinx
from dolfinx import fem, plot, geometry, mesh
from dolfinx.io import XDMFFile, gmshio
from ufl import (TestFunction, TrialFunction,jump, grad, inner, Measure, avg, CellDiameter,
FacetNormal, dx, dS, dot, lhs, rhs, FiniteElement, MixedElement)
from mpi4py import MPI
from petsc4py import PETSc
# Initializing parameters
t = 0 # Start time
T = 10 # End time
num_steps = 100 # Number of time steps
dt = (T-t)/num_steps # Time step size
#Geometrical parameters
Lb=0 # Left boundary
R= 1 # Radius of droplet
Rb = 10 # Right boundary
#Physical parameters
I_conc = 0.0 # Initial concentration inside
O_conc = 0.1 # Initial concentration outside
P = 5 # Partition coefficient
Di = 0.02 # Diffusion coefficient inside
Do = 0.2 #Diffusion coefficient outside
k = 1 # Interfacial resistance
### MESH GENERATION
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
model_rank = 0
comm = MPI.COMM_WORLD
gdim=1
gmsh.option.setNumber("Mesh.Algorithm", 7)
if comm.rank == model_rank:
gmsh.model.add("Linear_1D")
gmsh.model.setCurrent("Linear_1D")
Left = gmsh.model.occ.addPoint(Lb,0,0,tag=1)
Interface = gmsh.model.occ.addPoint(R,0,0,tag=2)
Right = gmsh.model.occ.addPoint(Rb,0,0,tag=3)
Line_I=gmsh.model.occ.addLine(Left,Interface,tag=4)
Line_O=gmsh.model.occ.addLine(Interface,Right,tag=5)
gmsh.model.occ.synchronize()
# add a physical group Inside
gmsh.model.addPhysicalGroup(gdim, [Line_I], tag=6)
gmsh.model.setPhysicalName(1, 6, "Inside")
# add a physical group Inside
gmsh.model.addPhysicalGroup(gdim, [Line_O], tag=7)
gmsh.model.setPhysicalName(1, 7, "Outside")
# add a physical group to interface
gmsh.model.addPhysicalGroup(gdim, [Interface], tag=8)
gmsh.model.setPhysicalName(1, 8, "Interface")
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",0.01)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",0.01)
gmsh.model.mesh.generate(gdim)
gmsh.write("Linear_1D_DG.msh")
domain, cell_markers, facet_markers= gmshio.model_to_mesh(gmsh.model,
comm, model_rank,
partitioner=dolfinx.mesh.create_cell_partitioner(dolfinx.mesh.GhostMode.shared_facet))
domain.name = "Linear_1D_DG"
cell_markers.name = f"{domain.name}_cells"
facet_markers.name = f"{domain.name}_facets"
gmsh.finalize()
# Creating xdmf to write domain and solution
xdmf = XDMFFile(domain.comm, "Linear_1D_DG.xdmf", "w")
xdmf.write_mesh(domain)
# Creating space function
V = fem.FunctionSpace(domain, ("DG", 0))
# Defining subdomains with mesh tags
cells_0 = cell_markers.find(6)
cells_1 = cell_markers.find(7)
D_ = fem.Function(V)
D_.x.array[cells_0] = np.full_like(cells_0, Di, dtype=PETSc.ScalarType)
D_.x.array[cells_1] = np.full_like(cells_1, Do, dtype=PETSc.ScalarType)
# Creating Trial and Test functions
ui = TrialFunction(V)
vi = TestFunction(V)
# mesh related functions
n = FacetNormal(domain)
d = CellDiameter(domain)
d_avg = (d('+') + d('-'))/2
# FRAP initial conditions
u_n = fem.Function(V)
u_n.x.array[cells_0] = np.full_like(cells_0, I_conc, dtype=PETSc.ScalarType)
u_n.x.array[cells_1] = np.full_like(cells_1, O_conc, dtype=PETSc.ScalarType)
# Define solution variable
uh = fem.Function(V)
uh.x.array[cells_0] = np.full_like(cells_0, I_conc, dtype=PETSc.ScalarType)
uh.x.array[cells_1] = np.full_like(cells_1, O_conc, dtype=PETSc.ScalarType)
### DEFINIG VARIATIONAL FORMULATION
F = (u_n-ui)*vi*dx - D_*dt/d_avg*dot(jump(ui,n),jump(vi,n))*dS
## BOUNDARY CONDITIONS
boundary_R =[(120, lambda x:np.isclose(x[0],R))]
fdim= domain.topology.dim - 1
facet_indices, facet_markers = [], []
for (marker, locator) in boundary_R:
facets = dolfinx.mesh.locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(domain, fdim,
facet_indices[sorted_facets], facet_markers[sorted_facets])
domain.topology.create_connectivity(domain.topology.dim-1,
domain.topology.dim)
xdmf.write_meshtags(facet_tag)
dS = Measure("dS", domain=domain, subdomain_data=facet_tag)
F += -P*jump(u_n)*jump(vi)*dS(120)
bc=[]
###
# Write initial solution to file
xdmf.write_function(uh, t)
### Finishing setting up variational formulation
a = fem.form(lhs(F))
L = fem.form(rhs(F))
A = fem.petsc.create_matrix(a)
b = fem.petsc.create_vector(L)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.BCGS)
solver.getPC().setType(PETSc.PC.Type.JACOBI)
###
for n in range(num_steps):
# update time
t += dt
# solve linear system
A.zeroEntries()
fem.petsc.assemble_matrix(A, a, bcs=bc)
A.assemble()
with b.localForm() as loc:
loc.set(0)
fem.petsc.assemble_vector(b, L)
fem.petsc.apply_lifting(b, [a], [bc])
b.ghostUpdate(PETSc.InsertMode.ADD_VALUES, PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, bc)
solver.solve(b, uh.vector)
uh.x.scatter_forward()
# update solution
u_n.x.array[:] = uh.x.array
# save solution
xdmf.write_function(uh, t)
xdmf.close()
Thank you in advance!
Marta