Vector heat equation

variational formulation
1

When solve this problem, I have trouble, here is the problem . Numerical Results for the Vector Heat Equation with EBC. Before presenting the results for the PPE
reformulation, we show the convergence results for the vector heat equation (VHE) with EBC as a benchmark.
This provides both: a baseline for the PPE reformulation convergence study, and a verification of the code.
Let the problem domain be Ω = [0,1) × [0, 1] with periodic b.c. applied in the x-direction and EBC in the
y-direction. Hence, the problem reads as
ut = ∆u + f in (0,1) ×(0,1)
n × u = 0;div u = 0 on [0,1) ×{0,1}
u(0; y) = u(1; y) for 0<y<1

which the problem is defined in the paper HIGH-ORDER METHODS FOR A PRESSURE POISSON EQUATIONREFORMULATION OF THE NAVIER-STOKES EQUATIONS WITH ELECTRIC BOUNDARY CONDITIONS section3.4.1

I write the code as the follow
#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(4pi*(x+y))((4y*(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)

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(1):

 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_)

but I got a complete wrong answer. I need help for it.

2

Please format your code using ``` to make it readable.

3

sorry I don’t get you mean, is it the follow?

‘’#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(4pi*(x+y))((4y*(1-y))**4)

u1=sp.diff(phsi,y,1)

#u1=pisin(2pi*y)sin(pix)**2

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 = RectangleMesh(Point(0, 0), Point(1, 1), 16, 16)

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(1):

 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_)''
4

As explained in this post: Solution of Mixed formulation

5

yeah, I follow the paper and write the code, but I dont get the right answer. Maybe I define a wrong boundary condition. Or I define a wrong space of the RT element. But I have no idea about this. So I want some help me to solve this issue.

6

My point was that how to use the Markdown syntax is explained in this post, such that your code becomes readable.

7

I try it but there is something need to change ``` and can run

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(1):
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_)