# UserExpression update in next time step gone wrong

**URL:** <https://fenicsproject.discourse.group/t/userexpression-update-in-next-time-step-gone-wrong/4880>\
**Category:** variational formulation\
**Created:** [January 11, 2021, 2:39pm UTC](https://fenicsproject.discourse.group/t/userexpression-update-in-next-time-step-gone-wrong/4880 "2021-01-11T14:39:23Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![heitorvc](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/heitorvc/32/754_2.png) [@heitorvc](https://fenicsproject.discourse.group/u/heitorvc)\
**Post date:** [January 11, 2021, 2:39pm UTC](https://fenicsproject.discourse.group/t/userexpression-update-in-next-time-step-gone-wrong/4880/1 "2021-01-11T14:39:23Z")

</div>

Hello,

I’m working on a free surface problem using the ALE method, and so far I’ve been able to update the surface curvature and normal vector by redefining my variational problem at each time step. I apologize as couldn’t figure out how to develop a MWE for this application, so I kindly ask to refer [here](https://github.com/heitorvitorc/masters_thesis/blob/main/passive_surface_fenics_v1.py) for the complete code. The code is still being developed, so it is still computationally inneficient (and a little fuzzy). I’m using the fenics tutorial for the navier-stokes equation using the IPCS approach, found [here](https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft08_navier_stokes_cylinder.py). I have the following stress boundary condition to represent the free passive surface:  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/a/a3f86147ad6d2981f565d27c1e1e1efe82569005.png)

Currently I’m redefining the variational problem at each time step using the python function below to apply the boundary condition:

```
# Interface stress 
def IST(sigma, kappa, n):
    return sigma*kappa*n

```

The variational formulation becomes

```
F1 = rho*dot((u - u_n) / k, v)*dx \
    + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
    + inner(TT(U, p_n, mu, I), DD(v))*dx \
    - dot(mu*nabla_grad(U)*n, v)*ds \
    - dot(f, v)*dx \
    + dot(p_n*n, v)*ds
b_int = inner(IST(sigma, kappa, n), v)*ds(surface_marker)
F1 = F1 - b_int

```

By using this method I generate physically coherent results, such as the velocity field and pressure field of a perturbed interface as described below. The results present a pressure difference due to the surface curvature, expressing the capillarity effects that are expected for this transport problem. The capillary effects will thus induce a fluid displacement in both on the surface and the bulk of the domain.

Velocity field:

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/d2baed36dc75e4f4c2eefe06b48eb77381c366f9.png)  
Pressure field:  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/4/48afe4c1fc18eabd66cc753caaad6b1bc38c1296.png)

Now i’m trying to implement a UserExpression to be updated at each time step, so that I don’t need to redefine the variational problem for every step of the temporal loop. I’m following this [post](https://fenicsproject.discourse.group/t/how-to-update-an-userexpression-for-the-next-time-step/4369/3), in special this [solution](https://fenicsproject.discourse.group/t/how-to-update-an-userexpression-for-the-next-time-step/4369/4) to implement a UserExpression class. I am also studying [this post](https://fenicsproject.discourse.group/t/update-expression-for-dirichletbc/2412) on how to update time-dependent boundary conditions, as it is also related to my question. I came up to the following class for my problem:

```
class InterfaceTension(UserExpression):
    def __init__ (self, sigma, kappa, n, **kwargs):
        super(). __init__ (kwargs)
        self.sigma = sigma
        self.kappa = kappa
        self.n = n

    def eval(self, values, x):
        values = self.sigma * self.kappa * self.n 

    def update(self, t, u):
        self.kappa = kappa
        self.n = n

    def value_shape(self):
        return (2,)

```

Then I instantiate the class by  
`IST = InterfaceTension(sigma, kappa, n)`

The variational form then becomes

```
F1 = rho*dot((u - u_n) / k, v)*dx \
    + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
    + inner(TT(U, p_n, mu, I), DD(v))*dx \
    - dot(mu*nabla_grad(U)*n, v)*ds \
    - dot(f, v)*dx \
    + dot(p_n*n, v)*ds

b_int = inner(IST, v)*ds(surface_marker)
F1 = F1 - b_int

```

After solving the problem for the first time-step, I need to update the new curvature and normal vector associated with the surface (inside the time loop):

```
# inside the temporal loop
kappa, n = curvature(mesh, ds, surface_marker) # see github to access this function
IST.kappa = kappa
IST.n = n

```

This approach returns a non-physical result, different from the results obtained by using the first approach (redefining the variational problem at each time step):

Velocity field from the second approach:

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/5/501c7093a41a2cf13b14e3d803866265f6a68198.png)  
Pressure field from the second approach:  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/f8f99b252f3ac0385d1c4f1ba0e4b3a90b754361.png)

Does anyone know what can be happening to create this difference between the results obtained using the two approaches? Honestly I have no idea why this is happening. I am very grateful for any insights.

Thanks in advance.

---

<div class="post-metadata">

**Author:** ![heitorvc](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/heitorvc/32/754_2.png) [@heitorvc](https://fenicsproject.discourse.group/u/heitorvc)\
**Post date:** [January 16, 2021, 1:20pm UTC](https://fenicsproject.discourse.group/t/userexpression-update-in-next-time-step-gone-wrong/4880/2 "2021-01-16T13:20:33Z")

</div>

Hello. I switched to a mixed element formulation and got coherent results.Thanks anyway.