# Weak form: Implementing source term - Strain history

**URL:** <https://fenicsproject.discourse.group/t/weak-form-implementing-source-term-strain-history/7204>\
**Category:** variational formulation\
**Created:** [December 2, 2021, 9:56am UTC](https://fenicsproject.discourse.group/t/weak-form-implementing-source-term-strain-history/7204 "2021-12-02T09:56:51Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![JoelTomasek](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/joeltomasek/32/2302_2.png) [@JoelTomasek](https://fenicsproject.discourse.group/u/JoelTomasek)\
**Post date:** [December 2, 2021, 9:56am UTC](https://fenicsproject.discourse.group/t/weak-form-implementing-source-term-strain-history/7204/1 "2021-12-02T09:56:51Z")

</div>

Dear all,

I would like to implement the following equation in a Navier Stokes/Convection-diffusion solver in FEniCS:  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/0/07f3b48607dca139b9ff813d187ceaca92be010b.png)

The goal of this equation is to get an insight in the “strain-history” of particles inside the velocity field.To that end, the scalar stress is defined as the Frobenius norm of the rate of strain tensor, which is added as a source term in the convection diffusion equation.

I tried to implement this, but I am not certain if my approach and specifically the variational formulation is correct. Could you tell me if at first glance I made any mistakes, or whether this is the correct approach? I added the code below:

```auto

V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
S = FunctionSpace(mesh,'P',1)

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
s = TrialFunction(S)
z = TestFunction(S)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
s_n = Function(S)
s_ = Function(S)

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

def strainNorm(u):
    return Constant(sqrt(2.)*4e-3)*sqrt(inner(epsilon(u),epsilon(u)))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))
	
# Define expressions used in variational forms
U = 0.5*(u_n + u)
n = FacetNormal(mesh)
f = Constant((0, 0))
fs = Constant(0.)
k = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
   + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Define variational problem for step 4

F4 = rho*((s - s_n) / k)*z*dx \
   + rho*dot(u_n, grad(s))*z*dx \
   - fs*z*dx
a4 = lhs(F4)
L4 = rhs(F4)

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')
	
    # Step 4: Scalar solve step
    A4 = assemble(a4)
    [bc.apply(A4) for bc in bcs]
    b4 = assemble(L4)
    [bc.apply(b4) for bc in bcs]
    solve(A4,s_.vector(),b4)

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)
    s_n.assign(s_)

```

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [December 2, 2021, 3:41pm UTC](https://fenicsproject.discourse.group/t/weak-form-implementing-source-term-strain-history/7204/2 "2021-12-02T15:41:01Z")

</div>

If `F4` is meant to correspond to the advection equation in the post, you would need to replace `fs` with \eta instead of `Constant(0)`. Treating the velocity explicitly, as in the advection term of `F4`, this would be

```auto
fs = strainNorm(u_n)

```

---

<div class="post-metadata">

**Author:** ![JoelTomasek](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/joeltomasek/32/2302_2.png) [@JoelTomasek](https://fenicsproject.discourse.group/u/JoelTomasek)\
**Post date:** [December 3, 2021, 9:53am UTC](https://fenicsproject.discourse.group/t/weak-form-implementing-source-term-strain-history/7204/3 "2021-12-03T09:53:58Z")

</div>

Thank you for your reply and checking the implementation. Indeed I meant to write fs = strainNorm(u\_n), but I forgot to make the edit.

---

<div class="post-metadata">

**Author:** ![JoelTomasek](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/joeltomasek/32/2302_2.png) [@JoelTomasek](https://fenicsproject.discourse.group/u/JoelTomasek)\
**Post date:** [December 3, 2021, 10:28am UTC](https://fenicsproject.discourse.group/t/weak-form-implementing-source-term-strain-history/7204/4 "2021-12-03T10:28:01Z")

</div>

Upon running this implementation in a benchmark case, I notice that I get much coarser results for the scalar solution than the results from this paper [Potential fluid mechanic pathways of platelet activation | SpringerLink](https://link.springer.com/article/10.1007/s10237-012-0417-4)  
With the implementation above, is it true that I only calculate strain for one point in a cell?

If I want to get a higher resolution for the strain computation and the scalar solution, is it sufficient to change S to FunctionSpace(mesh,‘P’, 2)?

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [December 3, 2021, 3:47pm UTC](https://fenicsproject.discourse.group/t/weak-form-implementing-source-term-strain-history/7204/5 "2021-12-03T15:47:46Z")

</div>

> I get much coarser results for the scalar solution than the results from this paper

It looks like the linked paper is using Lagrangian particles instead of an Eulerian discretization of advection. This is more complicated to implement, but introduces less diffusion. (See, e.g., [LEoPart](https://bitbucket.org/jakob_maljaars/leopart/src/master/) for a FEniCS-based implementation of Lagrangian particles.)

> If I want to get a higher resolution for the strain computation and the scalar solution, is it sufficient to change S to FunctionSpace(mesh,‘P’, 2)?

I don’t think this is likely to improve results without first working on other aspects of the discretization. The first-order time discretization may be introducing significant error and, for pure advection, it’s best to include some form of stabilization in the spatial discretization.

---

<div class="post-metadata">

**Author:** ![JoelTomasek](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/joeltomasek/32/2302_2.png) [@JoelTomasek](https://fenicsproject.discourse.group/u/JoelTomasek)\
**Post date:** [December 4, 2021, 12:42pm UTC](https://fenicsproject.discourse.group/t/weak-form-implementing-source-term-strain-history/7204/6 "2021-12-04T12:42:56Z")

</div>

Thank you for your suggestions. I have decided to take a look at the implementation of LEoPart. For my full application, I am using Oasis to solve incompressible NS in optimized fashion ([GitHub - mikaem/Oasis](https://github.com/mikaem/Oasis)). Is it possible to implement LEoPart in this formulation, or does it require that I solve NS with a discontinuous Galerkin approach?

In the meantime I will attempt to introduce stabilization in the spatial discretization and improve time discretization for the Eularian approach.  
Below I attached my new version of the code, with the addition of SUPG stabilization terms.  
I never implemented SUPG before, I am currently performing some tests to see how this implementation performs. Would you say that the code below provides a correct implementation?

```auto
# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
S = FunctionSpace(mesh,'P',1)

# Define inlet flow profile
plug_velocity=0.02e-6/(math.pi*0.001375**2)
inflow_profile = ('2*plug_velocity*velocity_scaling*((0.001375-x[1])*(0.001375+x[1])) / pow(0.001375, 2)', '0')
vely=Expression(inflow_profile, plug_velocity=plug_velocity,velocity_scaling=velocity_scaling, degree=2)

# Define boundary conditions
bcu_inflow = DirichletBC(V, vely, inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcs_inflow = DirichletBC(S, Constant(0.),inflow)
bcu = [bcu_inflow, bcu_walls]
bcp = [bcp_outflow]
bcs = [bcs_inflow]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
s = TrialFunction(S)
z = TestFunction(S)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
s_n = Function(S)
s_ = Function(S)

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))
def strainNorm(u):
    return Constant(sqrt(2.)*4e-3)*sqrt(inner(epsilon(u),epsilon(u)))
# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define expressions used in variational forms
U = 0.5*(u_n + u)
s_mid = 0.5*(s_n + s)
n = FacetNormal(mesh)
f = Constant((0, 0))
fs = strainNorm(u_n)
k = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)
h = CellDiameter(mesh)

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
   + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Define variational problem for step 4
F4 = ((s - s_n) / k)*z*dx \
   + dot(u_, grad(s_mid))*z*dx \
   - fs*z*dx
   
#Add SUPG stabilization terms
vnorm = sqrt(dot(u_, u_))
# Residual
r = (s - s_n)/k + (dot(u_, grad(s_mid)) - fs)
F4 += (h/(2.0*vnorm))*dot(u_, grad(z))*r*dx

a4 = lhs(F4)
L4 = rhs(F4)

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

s=s_n

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt
    if n%int(0.05*1/dt)==0:
        print(t)
    flowsec=0.04e-6
    baseflow = 0.06e-6
    vely.plug_velocity = (baseflow+flowsec*math.sin(4*math.pi*t-0.5*math.pi))/(math.pi*0.001375**2)

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')
	
    # Step 4: Scalar solve step
    A4 = assemble(a4)
    [bc.apply(A4) for bc in bcs]
    b4 = assemble(L4)
    [bc.apply(b4) for bc in bcs]
    solve(A4,s_.vector(),b4)

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)
    s_n.assign(s_)
    

```