Vector heat equation need help

variational formulation
1 
#vectorheat
from __future__ import print_function
from dolfin import *
import numpy as np
from math import *
import os, sys
import matplotlib.pyplot as plt
import sympy as sp
from sympy import exp, sin, pi, cos

T = 1          # final time
num_steps = 10000    # number of time steps
dt = T / num_steps # time step size

# define symbol expressions of u,p,f using sympyw23 
x,y,t = sp.symbols("x[0],x[1],t")
phsi=cos(t)*sin(4*pi*(x+y))*((4*y*(1-y))**4)
u1=sp.diff(phsi,y,1)
u2=-sp.diff(phsi,x,1)
f1=sp.diff(u1,t,1)-sp.diff(u1,x,2)-sp.diff(u1,y,2)
f2=sp.diff(u2,t,1)-sp.diff(u2,x,2)-sp.diff(u2,y,2)
# transform the expressions to ccode style
u1 = sp.printing.ccode(u1)
u2 = sp.printing.ccode(u2)
f1 = sp.printing.ccode(f1)
f2 = sp.printing.ccode(f2)
phsip = sp.printing.ccode(phsi)

class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        y[0] = x[0] - 1
        y[1] = x[1]
# Create mesh and define function spaces
mesh = UnitSquareMesh(32, 32)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Q = FunctionSpace(mesh, P1, constrained_domain=PeriodicBoundary())
Mh= FunctionSpace(mesh, "RT", 1, constrained_domain=PeriodicBoundary()) 
# Define trial and test functions  
u = TrialFunction(Mh)
v = TestFunction(Mh)
sigm =TrialFunction(Q)
sigmT =TestFunction(Q) 
sigm_=Function(Q)
u_=  Function(Mh)
f = Expression((f1,f2), degree=5, t=0)
u_e = Expression((u1,u2), degree=5,t=0)
u_n = interpolate(u_e, Mh)
k   = Constant(dt)

# Define variational problem for step 0
a0=inner(sigm,sigmT)*dx
L0=inner(u_n,curl(sigmT))*dx
# Define variational problem for step 1
F1=inner((u-u_n)/k,v)*dx+inner(curl(sigm_),v)*dx+inner(div(u),div(v))*dx-inner(f,v)*dx
a1=lhs(F1)
L1=rhs(F1)
# Time-stepping
t = 0
for n in range(num_steps):
     t += dt
     solve(a0 == L0, sigm_)
     f.t=t
     solve(a1 == L1, u_)
     u_e.t = t
     error=errornorm(u_e, u_, "L2")
     print('t = %.2f: error = %.3g' % (t, error))
     u_n.assign(u_)

I want to solve the vector heat equation, the above is the code and the eqution is from the paper HIGH-ORDER METHODS FOR A PRESSURE POISSON EQUATION REFORMULATION OF THE NAVIER-STOKES EQUATIONS WITH ELECTRIC BOUNDARY CONDITIONS https://arxiv.org/abs/2002.09801. In this paper, the author is use Fenics, so I implement it and want to get the same results. But I meet some trouble. The equation and varitional form as

1
11303×321 53.9 KB
2 .
After run the code , take dt=0.001, T=0.01, and I get
3

So I really need someone help me.

2

As your error seems to go to infinity, there is probably and error in your variational form, which causes instability. I would analyze the solution at the first step, and then try to reduce your example. As you have the initial condition, the solution after one solve should not change a lot. This is how you should start addressing your issue.

3

I do not have idea for this issue, I just follow the author idea from the paper. This problem make my work come difficulty, so I really want someone help me. Thanks a lot!

4

Have you even looked at what I suggested?
If you look at the error after 1 time-step, you observe that the solution is off the actual solution with about an magnitude of 0.1

Then, if you look at the error at later time steps, you observe instabilities in the top and bottom of your domain.

I’ve slightly modified your problem, using dolfin’s own differentation instead of sympy:

from dolfin import *
import numpy as np

dt = 1e-4
class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        y[0] = x[0] - 1
        y[1] = x[1]
# Create mesh and define function spaces
mesh = UnitSquareMesh(128, 128)

x = SpatialCoordinate(mesh)
t_c = Constant(0)
psi = cos(t_c)*sin(4*pi*(x[0]+x[1]))*(4*x[1]*(1-x[1]))**4
grad_psi = grad(psi)
u_e = as_vector((grad_psi[1], -grad_psi[0]))
f = diff(u_e, t_c)-div(grad(u_e))
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Q = FunctionSpace(mesh, P1, constrained_domain=PeriodicBoundary())
Mh= FunctionSpace(mesh, "RT", 1, constrained_domain=PeriodicBoundary())
# Define trial and test functions
u = TrialFunction(Mh)
v = TestFunction(Mh)
sigm =TrialFunction(Q)
sigmT =TestFunction(Q)
sigm_=Function(Q)
u_=  Function(Mh)
#f = Expression((f1,f2), degree=5, t=0)
#u_e = Expression((u1,u2), degree=5,t=0)
u_n = project(u_e, Mh)
k   = Constant(dt)

# Define variational problem for step 0
a0=inner(sigm,sigmT)*dx
L0=inner(u_n,curl(sigmT))*dx
# Define variational problem for step 1
F1=inner((u-u_n)/k,v)*dx+inner(curl(sigm_),v)*dx+inner(div(u),div(v))*dx-inner(f,v)*dx
a1=lhs(F1)
L1=rhs(F1)
# Time-stepping
t = 0
error_f = File("error.pvd")
errorfunc = Function(Mh)
while t < 10*dt:
    t += dt
    t_c.assign(t)
    solve(a0 == L0, sigm_)
    solve(a1 == L1, u_)
    error=np.sqrt(assemble(inner(u_-u_e, u_-u_e)*dx))
    print('t = %.2e: error = %.3g' % (t, error))
    u_n.assign(u_)
    errorfunc.assign(project(u_n-u_e, Mh))
    error_f << errorfunc
5

yeah, I just follow the authors idea and implement the code. At first I can not believe this thing even though I take small time steps, the error blow up very fast, at first I just think the code is wrong from the paper, but I can not find anything wrong. Now I will check the variational form, and analyze it stable or not.

6

one thing I think it seems the EBC bc condition, I do not apply in the finite element space. But I have no idea for this.