# The problem of equation coefficients being functions of the previous time step

**URL:** <https://fenicsproject.discourse.group/t/the-problem-of-equation-coefficients-being-functions-of-the-previous-time-step/15048>\
**Category:** dolfinx\
**Tags:** variational-form\
**Created:** [June 22, 2024, 2:50am UTC](https://fenicsproject.discourse.group/t/the-problem-of-equation-coefficients-being-functions-of-the-previous-time-step/15048 "2024-06-22T02:50:00Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![wanghao](https://avatars.discourse-cdn.com/v4/letter/w/96bed5/32.png) [@wanghao](https://fenicsproject.discourse.group/u/wanghao)\
**Post date:** [June 22, 2024, 2:50am UTC](https://fenicsproject.discourse.group/t/the-problem-of-equation-coefficients-being-functions-of-the-previous-time-step/15048/1 "2024-06-22T02:50:00Z")

</div>

Hi everyone!  
The coefficient **c** in a simple equation.  
For example:

```auto
(un,phin)=ufl.TrialFunctions(W)     
(v,fai)=ufl.TestFunctions(W)
c=10
FF1 = inner(delta*un,v)*dx-inner(delta*un_1,v)*dx
FF1 +=(-c*inner(un,v)*dx)
FF1 +=inner(delta*phin,fai)*dx-inner(delta*phin_1,fai)*dx
FF1 +=inner(phin,fai)*dx
aa1=ufl.lhs(FF1)
LL1=ufl.rhs(FF1)
a1=fem.form(aa1)
L1=fem.form(LL1)

```

**c** is constant .Where **un** is unknown. **un\_1** is known.  
If c is a function of un\_1 and un\_2, how should this code be rewritten？  
**c=c(un,un\_1)**

The complete code is as follows：

```auto
import numpy as np
import ufl
import ufl.coefficient
import dolfinx
from dolfinx import cpp as _cpp
from dolfinx import fem,io 
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
                         extract_function_spaces, form, 
                         locate_dofs_geometrical, locate_dofs_topological) 
from dolfinx.io import XDMFFile 
from dolfinx.mesh import (CellType, GhostMode, create_rectangle, 
                          locate_entities_boundary,create_unit_square,refine,locate_entities,meshtags) 
from dolfinx import mesh 
from ufl import div, dx, grad, inner,dot,outer,ds,dS,avg,jump,rot,cross 
from mpi4py import MPI 
from petsc4py import PETSc as _PETSc 
from petsc4py import PETSc 
from dolfinx.io import XDMFFile 
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,  
                               create_vector, set_bc) 
import matplotlib.pyplot as plt 
msh = mesh.create_interval(MPI.COMM_WORLD, 100,[0,1]) 
 
u_element = ufl.FiniteElement("CG", 'interval',1) 
phi_element = ufl.FiniteElement("CG",'interval', 1) 

TH = ufl.MixedElement([u_element, phi_element])
W = fem.FunctionSpace(msh, TH)

W0, map0 = W.sub(0).collapse()
W1, map1 = W.sub(1).collapse()

(un,phin)=ufl.TrialFunctions(W) #
(v,fai)=ufl.TestFunctions(W)    

#---------------------------------------------------------------------------------------------
un_1=fem.Function(W0)
phin_1=fem.Function(W1)

un_2=fem.Function(W0)
phin_2=fem.Function(W1)

def initrial(x):
    values = np.zeros((1,x.shape[1]))
    values[0, :] = 2
    return values

def initrial_phi(x):
    values = np.zeros((1,x.shape[1]))
    values[0, :] =1
    return values

un_1.interpolate(initrial)
phin_1.interpolate(initrial_phi)
un_2.interpolate(initrial)
phin_2.interpolate(initrial_phi)

def bound(x):
    return np.isclose(x[0],1)|np.isclose(x[0],0)

bcs=[]
                               
#-------------------------------------------------------------------------------------------------
from petsc4py import PETSc
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, 
                               create_vector, set_bc)
delta=0.001

c=10
FF1 = inner(delta*un,v)*dx-inner(delta*un_1,v)*dx
FF1 +=(-c*inner(un,v)*dx)
FF1 +=inner(delta*phin,fai)*dx-inner(delta*phin_1,fai)*dx
FF1 +=inner(phin,fai)*dx
aa1=ufl.lhs(FF1)
LL1=ufl.rhs(FF1)
a1=fem.form(aa1)
L1=fem.form(LL1)

# Convert a to the matrix form
A1 = assemble_matrix(a1, bcs)
A1.assemble()
b1 = create_vector(L1)
set_bc(b1,bcs)

solver1 = PETSc.KSP().create(msh.comm)
solver1.setOperators(A1)
solver1.setType("preonly")
solver1.getPC().setType("lu")
solver1.getPC().setFactorSolverType("mumps")
pc2 = solver1.getPC()
pc2.setHYPREType("boomeramg")#"boomeramg"
opts = PETSc.Options() 
opts["ksp_error_if_not_converged"] = 1
opts["mat_mumps_icntl_14"] = 80  
opts["mat_mumps_icntl_24"] = 1  
opts["mat_mumps_icntl_25"] = 0  
solver1.setFromOptions()

t=0
xdmfu=XDMFFile(msh.comm, "EXAMPLE_code/velocityTT.xdmf", "w")
xdmfphi=XDMFFile(msh.comm, "EXAMPLE_code/c.xdmf", "w")
xdmfu.write_mesh(msh)
xdmfphi.write_mesh(msh)

timestep=200
E=np.zeros(timestep)
fenli2=np.zeros(timestep)
wn=Function(W)
for n in range(timestep):
    t += 0.01
    A1.zeroEntries()
    fem.petsc.assemble_matrix(A1, a1, bcs=bcs)  
    A1.assemble()
    with b1.localForm() as b_loc:
        b_loc.set(0)
    assemble_vector(b1, L1)  
    apply_lifting(b1, [a1], [bcs])                 
    b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b1, bcs)
    solver1.solve(b1, wn.vector)
    u_h=wn.sub(0).collapse()
    phi_h=wn.sub(1).collapse()

    un_1.x.array[:]=u_h.x.array
    phin_1.x.array[:]=phi_h.x.array

    comm = u_h.function_space.mesh.comm
    error = form((0.5*phi_h **2+0.5*u_h** 2) * ufl.dx)
    
    E[n] = np.abs(comm.allreduce(fem.assemble_scalar(error), MPI.SUM))
    xdmfu.write_function(u_h,t)
    xdmfphi.write_function(phi_h,t)

    print('t is',n)
    print(E[n])
xdmfu.close()
xdmfphi.close()

```

Thank you to anyone who has commented！！

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 22, 2024, 4:05am UTC](https://fenicsproject.discourse.group/t/the-problem-of-equation-coefficients-being-functions-of-the-previous-time-step/15048/2 "2024-06-22T04:05:32Z")

</div>

> [@wanghao](#):
>
> c=c(un,un\_1)

What is the explicit formula for c?

Consider  
[https://jsdokken.com/dolfinx-tutorial/chapter2/nonlinpoisson\_code.html](https://jsdokken.com/dolfinx-tutorial/chapter2/nonlinpoisson_code.html)  
and the definition of q(u)

---

<div class="post-metadata">

**Author:** ![wanghao](https://avatars.discourse-cdn.com/v4/letter/w/96bed5/32.png) [@wanghao](https://fenicsproject.discourse.group/u/wanghao)\
**Post date:** [June 22, 2024, 5:24am UTC](https://fenicsproject.discourse.group/t/the-problem-of-equation-coefficients-being-functions-of-the-previous-time-step/15048/3 "2024-06-22T05:24:52Z")

</div>

Thanks dokken!  
q(u) Solved some of the problems.  
If I define：

```auto
def c(un_1):
    return c=np.sqrt(un_1**2)

```

AttributeError: ‘Sum’ object has no attribute ‘sqrt’

TypeError: loop of ufunc does not support argument 0 of type Sum which has no callable sqrt method.

**other case:**  
**c** is a scalar

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/5/5df08ad5ee51c6e6ba4c0102485f15e005b2be79.jpeg)  
How can I solve the above situation?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 22, 2024, 5:42am UTC](https://fenicsproject.discourse.group/t/the-problem-of-equation-coefficients-being-functions-of-the-previous-time-step/15048/4 "2024-06-22T05:42:08Z")

</div>

> [@wanghao](#):
>
> ```auto
> def c(un_1):
> return c=np.sqrt(un_1**2)
> 
> ```

This should be

```auto
def c(un_1):
    return ufl.sqrt(un_1**2)

```