Can a UFL statement be evaluated at each iteration?


249
views
0
7 months ago by

I'm trying to implement a time-dependent DG scheme that looks similar to here. The following bit of code is causing a bit of trouble for me, I think, because it appears that with each iteration, the expression for g gets longer and longer. Is it possible to condense the expression for g0 during the "g = g0" step where this expression would be evaluated rather than adding all of the UFL expressions together with each iteration? u0 is the solution from a previous iteration, if that helps. Thanks!

g0 = 1.0

for t in time:
    g = exp(1e-12*jump(u0)) + g0
    ...
    g0 = g
Community: FEniCS Project
I have another related question, and maybe my approach to the problem is incorrect, but what I'm trying to do is enforce a flux jump condition on an internal interface that depends on the difference in u across the interface with penalty terms similar to
$-\int_{\Gamma}\left\{\nabla v\right\}\cdot\left(\left[\left[u\right]\right]-g\right)dS+\frac{\alpha}{h}\int_{\Gamma}\left[\left[v\right]\right]\left(\left[\left[u\right]\right]-g\right)dS$Γ{v}·([[u]]g)dS+αh Γ[[v]]([[u]]g)dS .
I'm getting unreasonably large values for [[u]] and I'm wondering if it has something to do with [[u]] not being evaluated how I am thinking it is. Am I able to extract u('+') - u('-') on an interface $\Gamma$Γ that I can then use in another function for an update within each time step? I only need this term on the  $\Gamma$Γ and the unreasonably large values I get when I project it during the current update step in the code above make me think I'm approaching this problem incorrectly. Any help or suggestions are welcome and appreciated!
written 7 months ago by Dan Sweeney  

1 Answer


2
7 months ago by
Try replacing g0=g by
g0.assign(project(g, u0.ufl_function_space()))​

---
I was questioning the sense of the statement g=g0 ;)

After your comment, I realized I made an error in the final line and I fixed it in the code in the original post. I also tried your suggestion, but unfortunately I got the error: 
'Sum' object has no attribute 'assign'​
written 7 months ago by Dan Sweeney  
of course.. my mistake.
g = Function(u0.ufl_function_space())
g0 = Function(u0.ufl_function_space())
g0.interpolate(Constant(1.0))

time = range(5)
for t in time:
    g.assign(g0+project(exp(1e-12*jump(u0)), u0.ufl_function_space()))
#    print g.vector().array()[:]
    g0.assign(g)​
written 7 months ago by klunkean  
Thanks for your quick response! That looks like what I'm trying to do, but I get the following error
ufl.log.UFLException: Discontinuous type Coefficient must be restricted.​
Is this something that has to do with the jump() not being defined except on the edges of the elements?
written 7 months ago by Dan Sweeney  
1
A jump does not have to be restricted in principal. Your error probably gets thrown somewhere else...
This mwe runs for me:
from fenics import *
mesh = UnitSquareMesh(3,3)
ele = FiniteElement('DG', mesh.ufl_cell(), 0)

V = FunctionSpace(mesh, ele)

u0 = Function(V)

g = Function(u0.ufl_function_space())
g0 = Function(u0.ufl_function_space())
g0.interpolate(Constant(1.0))

time = range(5)
for t in time:
    g.assign(g0+project(exp(1e-12*jump(u0)), u0.ufl_function_space()))
    print g.vector().array()[:]
    g0.assign(g)​
written 7 months ago by klunkean  
Thanks for the help! I realized that I hadn't properly restricted g on the interface. My working example code is below.
from fenics import *
import numpy as np

ofilename = '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', 1)
u = TrialFunction(V)
v = TestFunction(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=2, domain=mesh)
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2
alpha = Constant(10.0)
gamma = Constant(alpha)
p0 = Constant(1.0)
u0 = interpolate(Constant(1.0), V)
g = Function(u0.ufl_function_space())
g0 = Function(u0.ufl_function_space())
g0.interpolate(Constant(1.0))

time = np.linspace(0.0, 1.0, 4)
for t in time:
    g.assign(g0+project(exp(1e-12*jump(u0)), u0.ufl_function_space()))
    print(g.vector().array()[:])
    g0.assign(g)
    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)-jump(g*n))*r*dS(1) \
              + alpha/h_avg*inner(jump(u)-jump(g), jump(v))*r*dS(1)
              
    Fu_bc = -inner(grad(v), (u-p0)*n)*r*ds(4) \
            - inner(v*n, grad(u))*r*ds(4) \
            + (gamma/h)*v*(u-p0)*r*ds(4)
            
    Fu_bc += -inner(grad(v), (u+p0)*n)*r*ds(6) \
             - inner(v*n, grad(u))*r*ds(6) \
             + (gamma/h)*v*(u+p0)*r*ds(6)   
             
    Fu = Fu_interior + Fu_edges + Fu_bc       
    a = lhs(Fu)
    l = rhs(Fu)
    u0 = Function(V)
    solve(a==l, u0)
plot(u0)​
written 7 months ago by Dan Sweeney  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »