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

dolfinx
1

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

(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：

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！！

2

What is the explicit formula for c?

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

3

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

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
image1597×499 67 KB

How can I solve the above situation?

4

This should be

def c(un_1):
    return ufl.sqrt(un_1**2)
1 Like