Need help with a model


184
views
0
3 months ago by

Hi!

I'm a MSc-student at University of Bergen (Norway) and I'm writing my thesis with use of FEniCS. I'm using it to study PDEs for fluid flow in poreous media.

My mesh is 2D, with (0, 1) x (0, 0.1) with 100 cells in the x-direction and 10 cells in the y-direction. As this is a time-dependent problem, I'm first using an initial condition to solve Darcy's equation with an initial condition for pressure and viscosity (dependent of concentration) and a fixed permeability tensor ((2,0),(0,1)). Darcy's equation:
\[
\begin{equation}
u = - K \nabla p
\end{equation}
\]
After the flux is calculated from this, the flux is used in convection–diffusion equation:
\[
\begin{equation}
\frac{\delta c}{\delta t} - \nabla \cdot \left( D \nabla c \right) + \nabla \cdot \left( u c \right) = f_c
\end{equation}
\]

Here \( p \) is pressure, \( u \) is flux (vector), \( K \) permability tensor and \( c \) consentration (number between 0 and 1, defined in each cell). For these calculations, CG-3 elements is used for the pressure, DRT-2 elements for flux and P-2 elements for consentration.

After these calculations is done, everything looks okey and the flow seams to be correct compared to my constructed solutions. But for our use, I want to project this 2D solution onto a 1D plot, where the y-axis is an averaged value for each x-point. Is there any examples on this or does anyone know how to do this most efficiently?

Community: FEniCS Project

2 Answers


1
3 months ago by
pf4d  
You can iterate through your 2D solution and average across a dimension.  Just call the function like this :

p = f(x,y)​


Otherwise, you can perform the integration across the domain numerically by solving a variational problem :
\[
  v(x,y) = \int_a^x u(s,y) ds = U(x,y) - U(a,y)
\]
where \(U\) is the anti-derivative of \(v\):
\[
  \frac{\partial U}{\partial x} = u(x,y).
\]
The integral at the start \(x = a\) will be zero, and so \(v(x,y) = U(x,y)\) and the BVP is therefore
\[
  \frac{\partial v}{\partial x} = u(x,y), \hspace{5mm} v(a,y) = 0.
\]
with associated variational problem consisting of finding \(v\) such that
\[
  \int_{\Omega} \frac{\partial v}{\partial x} \phi \mathrm{d}\Omega = \int_{\Omega} u \phi \mathrm{d}\Omega
\]
for all test functions \(\phi\).

An example in one dimension with \(u = \cos(x)\), analytical solution \(v = \sin(x)\) :

from fenics import *

mesh = IntervalMesh(1000,0,2*pi)
Q    = FunctionSpace(mesh, 'CG', 1)
x    = SpatialCoordinate(mesh)[0]
t    = mesh.coordinates()[:,0]

u    = interpolate(Expression('cos(x[0])', degree=2), Q)
s    = interpolate(Expression('sin(x[0])', degree=2), Q)
v    = TrialFunction(Q)
phi  = TestFunction(Q)

def left(x, on_boundary):
  return on_boundary and abs(x[0]) < 1e-14

# integral is zero on the left
bcs = DirichletBC(Q, 0.0, left)

a      = v.dx(0) * phi * dx
L      = u * phi * dx
v      = Function(Q)
solve(a == L, v, bcs)

uf = u.vector().array()[::-1]
vf = v.vector().array()[::-1]
sf = s.vector().array()[::-1]

r  = vf - sf

#===============================================================================
# plot :

from pylab import *
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset

mpl.rcParams['font.family']          = 'serif'
mpl.rcParams['text.usetex']          = True
mpl.rcParams['text.latex.preamble']  = ['\usepackage{fouriernc}']

fig = figure(figsize=(5,3.5))
ax  = fig.add_subplot(111)

ax.plot(t, uf, 'k',       lw=2.0, label=r"$u$")
ax.plot(t, vf, 'r--',     lw=2.0, label=r"$v$")
ax.plot(t, r,  '#880cbc', lw=2.0, label=r"$\epsilon$")

ax.set_xlim(0,2*pi)
leg = ax.legend(loc='upper center', fontsize='medium')
leg.get_frame().set_alpha(0.0)
ax.set_xlabel(r'$x$')
ax.grid()

tight_layout()
savefig("integration.pdf")


This is easily translatable to 2- or 3-dimensional space.

You want the average, so you'll have to divide the integral afterwards by the width of the dimension you integrated -- which may vary with the other dimension -- depending on your mesh.

edit: fixed the script above and description of example problem solved.

0
4 weeks ago by
Why are you using cos(x[0]), sin(x[0])? And if it were a 3D problem, would a TensorFunctionSpace be needed instead of a FunctionSpace?
There is a link "add comment" below my answer, this is the appropriate place to ask questions regarding my answer.

I use trig functions so that I can verify the numerical integration.

In 3D, why do you think a tensor function space should be used, and for what?  I do not follow.
written 4 weeks ago by pf4d  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »