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

139

views

0

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)
a = m1*p*v*dx(500) + m2*dot(grad(p),grad(v))*dx(500)\
+ ap*f1*p*v('+')*dS(777) + ap*f1*p*v('-')*dS(777)\
+ ap*f2*dot(grad(p),grad(v('+')))*dS(777)\
+ ap*f2*dot(grad(p),grad(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

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

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

a = m1*p*v*dx(500) + m2*dot(grad(p),grad(v))*dx(500)

, but when this part is added

+ ap*f1*p*v('+')*dS(777) + ap*f1*p*v('-')*dS(777)\ + ap*f2*dot(grad(p),grad(v('+')))*dS(777)\ + ap*f2*dot(grad(p),grad(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:

Thank you very much for your comments I will keep investigating your recommendations.

Meen

Thank you very much for your comments I will keep investigating your recommendations.

Meen

written
3 months ago by
Teeratorn Kadeethum

Please login to add an answer/comment or follow this question.

`p_n`