Jump() in Expression over interface in DG Poisson


581
views
0
7 months ago by

I have a question regarding how to implement a jump condition in FEniCS over an internal boundary. I want to use jump(u) within the expression for g (simplified in the code below). I want to use jump(u) on the interface marked '1' within an expression containing an exponential--in a manner similar to the expression for g below. Is this possible? If not, would it be better to extract the nodal values along interface '1', do my computation on the array? The error I get with g assigned as it is in the code below is "TypeError: expected default arguments for member variables to be scalars or GenericFunctions."

from fenics import *
import matplotlib.pyplot as plt
import numpy as np

ofilename = 'model_axisymm2D.xml'
mesh = Mesh(ofilename)
subdomains = MeshFunction("size_t", mesh, "%s_physical_region.xml"%ofilename.split('.')[0])
boundaries = MeshFunction("size_t", mesh, "%s_facet_region.xml"%ofilename.split('.')[0])

V = FunctionSpace(mesh, 'DG', 2)
v = TestFunction(V)
u = TrialFunction(V)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dS = Measure('dS', domain=mesh, subdomain_data=boundaries)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)

r = Expression('x[0]', degree=1, domain=mesh)
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2
alpha = Constant(10.0)
gamma = Constant(alpha)

VOLT = 40.0

g = Expression('alpha*exp(U*U)', alpha=1.0, U = jump(u), degree=2)

Fu_interior = inner(grad(v), grad(u))*r*dx

Fu_edges = -inner(jump(v,n), avg(grad(u)))*r*dS(0) \
          - inner(avg(grad(v)), jump(u,n))*r*dS(0) \
          + alpha/h_avg*inner(jump(u), jump(v))*r*dS(0)
          
Fu_edges += - inner(avg(grad(v)), jump(u,n)-g*n('+'))*r*dS(1) \
          + alpha/h_avg*inner(jump(u)-g, jump(v))*r*dS(1)

Fu_bc = -inner(grad(v), (u-VOLT)*n)*r*ds(4) \
        - inner(v*n, grad(u))*r*ds(4) \
        + (gamma/h)*v*(u-VOLT)*r*ds(4)

Fu_bc += -inner(grad(v), (u+VOLT)*n)*r*ds(6) \
         - inner(v*n, grad(u))*r*ds(6) \
         + (gamma/h)*v*(u+VOLT)*r*ds(6)
         
Fu = Fu_interior + Fu_edges + Fu_bc
         
a = lhs(Fu)
l = rhs(Fu)

u = Function(V)

solve(a==l, u)

p = plot(u)
plt.colorbar(p)
Community: FEniCS Project

1 Answer


2
7 months ago by
Update: I think I found the answer to my own problem here. Something to the tune of
g = exp(1e-9*jump(u0)*jump(u0))​
 works--I just have to use the solution from a previous timestep u0 to get my jump condition to work, which is not in the MWE I posted above.
Dear Dan,
I was just curious about the physics of your problem.
What is it?
Thanks in advance.
written 7 months ago by Murilo Moreira  
1
I was playing around with a Boltzmann distribution, if it helps.
written 7 months ago by Dan Sweeney  
I see!
Thank you very much!
written 7 months ago by Murilo Moreira  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »