### Mixed dimensional fluid flow in porous media - Discontinuous type Argument must be restricted

139
views
0
3 months ago by
I model a flow in porous media that has one fracture as attached:

// Characteristic size
Lm = 0.1;

// Size of the block (m)
Lx = 0.199;
Ly = 0.1;

// Fracture aperture (0.025cm)
fa = 2.5e-5;

// Density of mesh near Points
lcf = fa;

lcm = Lm/5;

// Geometrical entities = 1d and 2D regions indices
// These are the regions in CSMP
// #2D_REGIONS
MATRIX = 500;
// #2D_REGIONS

// #1D_REGIONS
LEFT = 1;
RIGHT = 2;
BOTTOM = 3;
TOP = 4;
OPEN_FRACTURES = 777;
CLOSED_FRACTURES = 778;
// #1D_REGIONS

// Points
Point(1) = {0,  0, 0, lcm};
Point(2) = {Lx, 0, 0, lcm};
Point(3) = {Lx,  Ly, 0, lcm};
Point(4) = {0,  Ly, 0, lcm};
Point(5) = {0.1,  0.02, 0, lcm};
Point(6) = {0.1,  0.08, 0, lcm};

// Define the bounded domains
Line(501) = {1, 2};
Line(502) = {2, 3};
Line(503) = {3, 4};
Line(504) = {4, 1};
Line(505) = {5, 6};

Line Loop(600) = {501, 502, 503, 504};
Plane Surface(600) = {600};

Line{505} In Surface {600};	// To include the faults in the triangulation

// Gradual mesh refinement starting from BOTTOM_FRAC
Field[1] = Attractor;
Field[1].NNodesByEdge = 10;
Field[1].EdgesList = {505};
//Field[1].EdgesList = {507, 505, 506, 504, 503, 508, 501, 502};

Field[2] = Threshold;
Field[2].IField = 1;
Field[2].LcMin = 50 * fa;	// lcf
Field[2].LcMax = lcm;
Field[2].DistMin = fa;
Field[2].DistMax = Lm;

Background Field = 2;

Physical Line(LEFT) = {504};
Physical Line(RIGHT) = {502};
Physical Line(TOP) = {503};
Physical Line(BOTTOM) = {501};

Physical Line(OPEN_FRACTURES) = {505};

Physical Surface(MATRIX) = {600};
​

Then the python script is attached as follows:

from dolfin import *
import time
import os

#parameters
T = 100.0
num_steps = 200
dt = T/num_steps
ct = 2.e-9 #Pa-1
rho = 900 #kg/m3
phi = 0.3 #fraction

m1 = phi*ct*rho/dt

k_m = 1.e-17  #m2
#k_f = 1.e-18  #m2
mu = 2.e-3 #Pa.s

m2 = rho*k_m/mu

cf = 1.e-9 #Pa-1
phif = 1.0 #fraction
f1 = phif*cf*rho/dt

ap = 1.e-4 #m
k_f = (ap**2)/12 #m2
f2 = rho*k_f/mu

os.system('dolfin-convert test_1fr.msh test.xml')

mesh = Mesh("test.xml")

subdomains = MeshFunction("size_t",mesh,"test_physical_region.xml")
boundaries = MeshFunction("size_t",mesh,"test_facet_region.xml")

dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dS = Measure('dS', domain=mesh, subdomain_data=boundaries)

g_B = Constant(0.0)
g_T = Constant(0.0)

ul = Constant(19.e6)
ur = Constant(5.e6)

f = Constant(0.0)

V = FunctionSpace(mesh, "Lagrange", 1)

bcs = [DirichletBC(V, ul, boundaries, 1),
DirichletBC(V, ur, boundaries, 2)]

p_0 = Constant(5.e6)
p_n = project(p_0,V)

p =  TrialFunction(V)
v = TestFunction(V)

+ ap*f1*p*v('+')*dS(777) + ap*f1*p*v('-')*dS(777)\

l = (m1*p_n + f)*v*dx(500) + g_T*v*ds(4) + g_B*v*ds(3)\
+ (ap*f1*p_n + f)*v('+')*dS(777) + (ap*f1*p_n + f)*v('-')*dS(777)

vtkfile = File('test1/solution.pvd')
#vtkfile << (u_n,t)

p = Function(V)
t = 0

for n in range(num_steps):

t += dt

solve(a == l, p, bcs)

vtkfile << (p,t)

p_n.assign(p)

The fracture is modelled as a line element. I guess the problem is when I assemble 'a'.

The problem occurs as 'Discontinuous type Argument must be restricted'. Do you have any recommendations to fix this issue?

Thank you!

Meen

Community: FEniCS Project
you need to choose a side (+ or -) for p_n
written 3 months ago by Nate
Dear Nate:

Thank you very much for your help!

I chose the side for p_n and also used avg(v); however, the problem persisted. If I have only this part of the code, the script works fine

, but when this part is added

+ ap*f1*p*v('+')*dS(777) + ap*f1*p*v('-')*dS(777)\

, the problem occurs.

Do you have any further recommendations?

Thank you very much!!

Meen

written 3 months ago by Teeratorn Kadeethum
1
You need to choose a side (+ or -) for all functions defined on your domain. This includes p. Examine the DG demos.
written 3 months ago by Nate
I think I fix it; thank you for your help!!
written 3 months ago by Teeratorn Kadeethum
Just a comment:
v lives in V = FunctionSpace(mesh, "Lagrange", 1).
That means v is a continuous function and hence v('-') does not make that much sense.

If you want a discontinous function you should take V = FunctionSpace(mesh, "DG", 1).
written 3 months ago by Lukas O.
Dear Lucus: