### Can a UFL statement be evaluated at each iteration?

249

views

0

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

### 1 Answer

2

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:

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.

$-\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!