Question with enriching basis function near boundary


129
views
0
10 weeks ago by
Dear friends,

I am doing some simulation with fluid. The resolution of boundary layer is important for solution. In order to improve the accuracy, I want to enrich the basis function near boundary with some special exponential function. I wonder whether this is possible to implement with FEniCS. Hope you reply soon. Thank you for your help all the time.


Yours sincerely,
Qiming Zhu
Community: FEniCS Project

1 Answer


1
10 weeks ago by
There was some work a while back on partition of unity methods for adding enrichment functions into FEniCS,

https://repository.tudelft.nl/islandora/object/uuid:49e4dcc2-d22c-429e-a70d-73c30d72ba4a

but it looks pretty involved, with modifications of DOLFIN and FFC that I don't think were ever merged back into the main project.  (Maybe someone with more knowledge can correct me here.)

On the other hand, if there are only a few "special" functions for the whole problem, you could try to jury-rig something with "Real" type elements, e.g.,

from dolfin import *

# You need really good quadrature to fully resolve the special exponential
# function.
dx = dx(metadata={"quadrature_degree":20})

mesh = UnitIntervalMesh(5)
VE = FiniteElement("Lagrange",mesh.ufl_cell(),1)
RE = FiniteElement("Real",mesh.ufl_cell(),0)
V = FunctionSpace(mesh,MixedElement([VE,RE]))
V_unenriched = FunctionSpace(mesh,VE)

u_mixed = TrialFunction(V)
v_mixed = TestFunction(V)

# advection--diffusion parameters
kappa = Constant(0.01)
u_adv = Constant((1.0,))

# Peclet number
Pe = sqrt(inner(u_adv,u_adv))/kappa

# Define a special function to get an accurate boundary layer
x = SpatialCoordinate(mesh)[0]
specialFunction = exp(Pe*x) - (1.0-x) - x*exp(Pe)

# Helper functions to extract parts from elements of the mixed function space
def linearPart(u_mixed):
    u_linear,_ = split(u_mixed)
    return u_linear
def specialPart(u_mixed):
    _,u_specialCoeff = split(u_mixed)
    return u_specialCoeff*specialFunction
def mixedToFunc(u_mixed):
    return linearPart(u_mixed) + specialPart(u_mixed)

# What to plug into the variational problem
u = mixedToFunc(u_mixed)
v = mixedToFunc(v_mixed)

# Define variational forms for advection--diffusion (with no stabilization)
def a(u,v):
    return kappa*inner(grad(u),grad(v))*dx + inner(grad(u),u_adv)*v*dx
def L(v):
    return Constant(0.0)*v*dx

lhs = a(u,v)
rhs = L(v)

def getBCs(V):
    bc_left = DirichletBC(V, 0.0,"on_boundary && near(x[0],0.0)")
    bc_right = DirichletBC(V, 1.0,"on_boundary && near(x[0],1.0)")
    return [bc_left, bc_right]

# Solve for enriched approximate solution
u_mixed = Function(V)
solve(lhs==rhs,u_mixed,getBCs(V.sub(0)))

# Solve using the linear space only, for comparison:
u_unenriched = TrialFunction(V_unenriched)
v_unenriched = TestFunction(V_unenriched)
lhs = a(u_unenriched,v_unenriched)
rhs = L(v_unenriched)
u_unenriched = Function(V_unenriched)
solve(lhs==rhs,u_unenriched,getBCs(V_unenriched))

# Refined mesh for plotting
refinedMesh = UnitIntervalMesh(1000)
V_refined = FunctionSpace(refinedMesh,"Lagrange",1)

# Plot results, projected to the refined mesh:
import matplotlib.pyplot as plt
plot(project(mixedToFunc(u_mixed),V_refined))
plot(project(u_unenriched,V_refined))
#plot(project(specialPart(u_mixed),V_refined))
#plot(project(linearPart(u_mixed),V_refined))
plt.autoscale()
plt.show()
​
Thank you for your detailed reply very much. I appreciate your help all the time. I would try to implement my problem with your brilliant idea and "Real" type element. Thank you for your help.
written 10 weeks ago by zqm1992  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »