Magnetic fluid equation of porous media

Hey guys,
I am doing a program design of a porous medium magnetic fluid model. Now there is a problem that needs help. The following is my program code。

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

#setting parameters
M = 1.0
Lambda = 1.0
K = 1.0   # K^{-1} in paper
niu_w = 1.0
niu_f = 2.0
miu = 1.0
tau = 1.0
beta = 1.0
B = 1.0        # SAV constant for ensuring sqrt{}  meaningful
k_w = 0
k_f = 1.0
EPS = 0.5

# Artifical force
f_phi = Expression('cos(pi*x[0])*cos(pi*x[1])*cos(t) + 2*M*pow(pi, 2)*cos(pi*x[0])*cos(pi*x[1])*sin(t) - pi*pow(cos(pi*x[0]), 2)*sin(pi*x[1])*pow(sin(t), 2)*cos(pi*(x[1] + 0.5)) - pi*cos(pi*x[1])*pow(sin(pi*x[0]), 2)*pow(sin(t), 2)*sin(pi*(x[1] + 0.5))', degree=2, M=M, t=0)

f_w = Expression('cos(pi*x[0])*cos(pi*x[1])*sin(t) - 2*Lambda*pow(pi, 2)*cos(pi*x[0])*cos(pi*x[1])*sin(t) - (Lambda*cos(pi*x[0])*cos(pi*x[1])*sin(t)*(pow(cos(pi*x[0]),2)*pow(cos(pi*x[1]),2)*pow(sin(t),2)- 1))/pow(EPS,2)', degree=2, Lambda=Lambda, EPS=EPS, t=0)

f_u = Expression(('2*sin(t)*(2*x[1]-1)-miu*sin(t+x[0])*sin(t)+sin(pi*x[0])*cos(t)*sin(pi*(x[1]+0.5))-pi*cos(pi*x[0])*pow(cos(pi*x[1]),2)*sin(pi*x[0])*pow(sin(t),2)+K*sin(pi*x[0])*sin(t)*sin(pi*(x[1]+0.5))*(1/(exp(-(cos(pi*x[0])*cos(pi*x[1])*sin(t))/EPS)+1)+1)',
'2*sin(t)*(2*x[0]-1)-miu*sin(t+x[1])*sin(t)+cos(pi*x[0])*cos(t)*cos(pi*(x[1]+0.5))-pi*pow(cos(pi*x[0]),2)*cos(pi*x[1])*sin(pi*x[1])*pow(sin(t),2)+K*cos(pi*x[0])*sin(t)*cos(pi*(x[1]+0.5))*(1/(exp(-(cos(pi*x[0])*cos(pi*x[1])*sin(t))/EPS)+1)+1)'), degree=2, K=K, EPS=EPS, miu=miu, t=0)

f_m = Expression(('cos(t + x[1]) + (sin(t + x[1]) - (x[1]*sin(t))/(exp(-(cos(pi*x[0])*cos(pi*x[1])*sin(t))/EPS) + 1))/tau - beta*sin(t + x[0])*(x[1]*sin(t + x[0])*sin(t) - sin(t + x[1])*sin(t)*(x[0] - 0.5)) + cos(pi*x[0])*cos(t + x[1])*sin(t)*cos(pi*(x[1] + 0.5))',
'cos(t + x[0]) + (sin(t + x[0]) - (sin(t)*(x[0] - 0.5))/(exp(-(cos(pi*x[0])*cos(pi*x[1])*sin(t))/EPS) + 1))/tau + beta*sin(t + x[1])*(x[1]*sin(t + x[0])*sin(t) - sin(t + x[1])*sin(t)*(x[0] - 0.5)) + sin(pi*x[0])*cos(t + x[0])*sin(t)*sin(pi*(x[1] + 0.5))'), degree=2, EPS=EPS, beta=beta, tau=tau, t=0)

h_a = Expression(('sin(t+x[1])+x[1]*sin(t)', 'sin(t+x[0])+(x[0]-0.5)*sin(t)'), degree=2, t=0)

h_b = Expression(('cos(t+x[1])+x[1]*cos(t)', 'cos(t+x[0])+(x[0]-0.5)*cos(t)'), degree=2, t=0)

# Exact solution
phi_e = Expression('sin(t)*cos(pi*x[0])*cos(pi*x[1])', degree=2, t=0)       #exact solution of theta
w_e = Expression('sin(t)*cos(pi*x[0])*cos(pi*x[1])', degree=2, t=0)       #exact solution of omega
u_e = Expression(('sin(t)*sin(pi*x[0])*sin(pi*(x[1]+0.5))', 'sin(t)*cos(pi*x[0])*cos(pi*(x[1]+0.5))'), degree=2, t=0)   # exact solution of u
p_e = Expression('sin(t)*(2*x[0]-1)*(2*x[1]-1)', degree=2, t=0)               #exact solution of p
m_e = Expression(('sin(t+x[1])', 'sin(t+x[0])'), degree=2, t=0)               #exact solution of m
varphi_e = Expression('(x[0]-0.5)*x[1]*sin(t)', degree=2, t=0)                #exact solution of varphi
h_e = Expression(('x[1]*sin(t)', '(x[0]-0.5)*sin(t)'), degree=2, t=0)         #exact solution of h

# Denfine dt and mesh
T = 0.5
m = 32       #time step
dt = 1.0/m
k = Constant(dt)    # \delta t
n = 32             # space mesh
mesh = UnitSquareMesh(n, n)
domain = 1.0 * 1.0       #domain

# Denfine the finite element and the mixed functionspace
P1 = FiniteElement("CG", mesh.ufl_cell(), 1)        # P1 element
TH1 = MixedElement([P1, P1])             # Mixed element P1*P1
W1 = FunctionSpace(mesh, TH1)            # Mixed functionspace P1*P1

# Denfin the functionspace
P1Space = FunctionSpace(mesh, "CG", 1)            #P1  Space
P1_dgSpace = VectorFunctionSpace(mesh, "DG", 1)   #P1-dg  Space

P1VectorSpace = VectorFunctionSpace(mesh, "CG", 1)   #P1-vector  Space
P2Space = FunctionSpace(mesh, "CG", 2)            #P2 Space

# solve the average
int = Function(P2Space)
Unit = Constant(1.0)
average = (1.0/domain)*dot(int, Unit)*dx              #solve (1/|Omega|)*int_{Omega} ... dx

#Define  function
def F(a):
return (1.0/4)*(pow(a, 2)-1)*(pow(a, 2)-1)      # define the F(phi)
def f(a):
return a*(pow(a, 2)-1)   #define the f(fai)
def vcross(x, y):
return x[0] * y[1] - x[1] * y[0]  # 向量叉乘
# def h(a):
def niu(a):
return niu_w + (niu_f-niu_w)*(1.0/(1.0+exp(-(a/EPS))))     #define the \nu(phi)
def kappa(a):
return k_w+(k_f-k_w)*(1.0/(1.0+exp(-(a/EPS))))         #define the \kappa(phi)
def intergal(a):
return dot(F(a), 1)*dx                  #define      int_{Omega} F(phi)*dx

# Define the initial value
phi0 = interpolate(phi_e, P1Space)         # phi0
w0 = interpolate(w_e, P1Space)
# w0_trial = TrialFunction(P1Space)
# y0_test = TestFunction(P1Space)
# w0 = Function(P1Space)                     # w0
# aw0 = dot(w0_trial, y0_test)*dx
# solve(aw0 == Lw0, w0)                      # solve new w0

u0 = interpolate(u_e, P1_dgSpace)        # u0
p0 = Constant(0)                         # p0
m0 = interpolate(m_e, P1VectorSpace)     # m0

ha0 = interpolate(h_a, P1VectorSpace)
varphi0_trial = TrialFunction(P2Space)
psi0_test = TestFunction(P2Space)
varphi0 = Function(P2Space)
solve(avarphi0 == Lvarphi0, varphi0)

int.assign(varphi0)
value0 = assemble(average)
varphi0.vector()[:] = varphi0.vector() - value0

r0 = sqrt( (1.0/(EPS*EPS))*assemble(intergal(phi0)) + B )   # r0
print('r0= ', r0)

# Define the old value
phi_old = Function(P1Space)        # phi^{n}
w_old = Function(P1Space)          # w^{n}
u_old = Function(P1_dgSpace)       # u^{n}
p_old = Function(P2Space)          # p^{n}
m_old = Function(P1VectorSpace)    # m^{n}
varphi_old = Function(P2Space)     # varphi^{n}
r_old = 1.0

phi_a = Function(P1Space)          # phi_a^{n+1}
phi_b = Function(P1Space)          # phi_B^{n+1}
w_a = Function(P1Space)            # w_a^{n+1}
w_b = Function(P1Space)            # w_b^{n+1}
utilde_a = Function(P1_dgSpace)      # save utilde_a
utilde_b = Function(P1_dgSpace)        # save utilde_b
varphi_a = Function(P2Space)           #save varphi_a
varphi_b = Function(P2Space)          #save varphi_b
m_a = Function(P1VectorSpace)           #save m_a
m_b = Function(P1VectorSpace)           #save m_b
phi = Function(P1Space)            # save phi^{n+1}
w = Function(P1Space)            # save w^{n+1}
utilde = Function(P1_dgSpace)      #save utilde^{n+1}
varphi = Function(P2Space)         #save varphi^{n+1}
m = Function(P1VectorSpace)        #save m^{n+1}
p = Function(P2Space)             #save p^{n+1}
u = Function(P1_dgSpace)          #save u^{n+1}

S_n = 1.0
r = 1.0

# Define trial and test function
(phia_trial, wa_trial) = TrialFunctions(W1)
(X_test, Y_test) = TestFunctions(W1)
G1 = Function(W1)             # save phi_a^{n+1}  and w_a^{n+1}
Lphi_a = (1.0/k)*dot(phi_old, X_test)*dx + dot(f_phi, X_test)*dx + dot(f_w, Y_test)*dx

(phib_trial, wb_trial) = TrialFunctions(W1)
G2 = Function(W1)             # save phi_b^{n+1}  and w_b^{n+1}
Lphi_b = (1.0/S_n)*dot(u_old*phi_old, grad(X_test))*dx + (1.0/S_n)*(Lambda/(EPS*EPS))*dot(f(phi_old), Y_test)*dx

Phn_trial = TrialFunction(P1_dgSpace)
n_test = TestFunction(P1_dgSpace)
Ph_n = Function(P1_dgSpace)
aPh = dot(Phn_trial, n_test)*dx

utildea_trial = TrialFunction(P1_dgSpace)
v_test = TestFunction(P1_dgSpace)
autilde_a = (1.0/k)*dot(utildea_trial, v_test)*dx + K*dot(niu(phi_old)*utildea_trial, v_test)*dx
Lutilde_a = (1.0/k)*dot(u_old, v_test)*dx - dot(grad(p_old), v_test)*dx + dot(f_u, v_test)*dx

utildeb_trial = TrialFunction(P1_dgSpace)
autilde_b = (1.0/k)*dot(utildeb_trial, v_test)*dx + K*dot(niu(phi_old)*utildeb_trial, v_test)*dx

varphia_trial = TrialFunction(P2Space)
psi_test = TestFunction(P2Space)

varphib_trial = TrialFunction(P2Space)

ma_trial = TrialFunction(P1VectorSpace)
z_test = TestFunction(P1VectorSpace)
am_a = (1.0/k)*dot(ma_trial, z_test)*dx + (1.0/tau)*dot(ma_trial, z_test)*dx
Lm_a = (1.0/k)*dot(m_old, z_test)*dx + (1.0/tau)*dot(kappa(phi_old)*grad(varphi_a), z_test)*dx + dot(f_m, z_test)*dx

mb_trial = TrialFunction(P1VectorSpace)
am_b = (1.0/k)*dot(mb_trial, z_test)*dx + (1.0/tau)*dot(mb_trial, z_test)*dx

eta_a = (1.0/(2*S_n*pow(EPS, 2)))*dot(f(phi_old), phi_a)*dx - (k/(2*Lambda*S_n))*dot(u_old*phi_old, grad(w_a))*dx + (k/(2*Lambda*S_n))*dot(phi_old*grad(w_old), utilde_a)*dx \
Eta_a = assemble(eta_a)

eta_b = (1.0/(2*S_n*pow(EPS, 2)))*dot(f(phi_old), phi_b)*dx - (k/(2*Lambda*S_n))*dot(u_old*phi_old, grad(w_b))*dx + (k/(2*Lambda*S_n))*dot(phi_old*grad(w_old), utilde_b)*dx \
Eta_b = assemble(eta_b)

eta_c = (1.0/(2*S_n*pow(EPS, 2)))*dot(f(phi_old), phi_old)*dx
Eta_c = assemble(eta_c)

phi_trial = TrialFunction(P1Space)
x_test = TestFunction(P1Space)
aphi = dot(phi_trial, x_test)*dx
Lphi = dot(phi_a + r*phi_b, x_test)*dx

w_trial = TrialFunction(P1Space)
y_test = TestFunction(P1Space)
aw = dot(w_trial, y_test)*dx
Lw = dot(w_a + r*w_b, y_test)*dx

utilde_trial = TrialFunction(P1_dgSpace)
autilde = dot(utilde_trial, v_test)*dx
Lutilde = dot(utilde_a + r*utilde_b, v_test)*dx

varphi_trial = TrialFunction(P2Space)
avarphi = dot(varphi_trial, psi_test)*dx
Lvarphi = dot(varphi_a + r*varphi_b, psi_test)*dx

m_trial = TrialFunction(P1VectorSpace)
am = dot(m_trial, z_test)*dx
Lm = dot(m_a + r*m_b, z_test)*dx

p_trial = TrialFunction(P2Space)
q_test = TestFunction(P2Space)

u_trial = TrialFunction(P1_dgSpace)
au = dot(u_trial, v_test)*dx

phi_old.assign(phi0)
w_old.assign(w0)
u_old.assign(u0)
p_old.assign(p0)
varphi_old.assign(varphi0)
m_old.assign(m0)
r_old = r0

# Work in cycle
t = dt
while t < T + 1.0e-8:
print('t= ', t)

f_phi.t = t
f_w.t = t
f_u.t = t
f_m.t = t
h_a.t = t
h_b.t = t

S_n = sqrt( (1.0/(EPS*EPS))*assemble(intergal(phi_old)) + B)  # define S_n

solve(aphi_a == Lphi_a, G1)             # solve phi_a^{n+1} and w_a^{n+1}
(phia_up, wa_up) = G1.split(True)
phi_a.assign(phia_up)
w_a.assign(wa_up)

solve(aphi_b == Lphi_b, G2)             # solve phi_b^{n+1} and w_b^{n+1}
(phib_up, wb_up) = G2.split(True)
phi_b.assign(phib_up)
w_b.assign(wb_up)

solve(aPh == LPh, Ph_n)                     # solve the L^{2} projection P(h)

solve(autilde_a == Lutilde_a, utilde_a)   # solve utilde_a^{n+1}

solve(autilde_b == Lutilde_b, utilde_b)   # solve utilde_b^{n+1}

solve(avarphi_a == Lvarphi_a, varphi_a)   # solve varphi_a^{n+1}
int.assign(varphi_a)
value1 = assemble(average)
varphi_a.vector()[:] = varphi_a.vector() - value1     # # varphi_a在[0,1]上的积分为0

solve(avarphi_b == Lvarphi_b, varphi_b)   # solve varphi_b^{n+1}
int.assign(varphi_b)
value2 = assemble(average)
varphi_b.vector()[:] = varphi_b.vector() - value2   # varphi_b在[0,1]上的积分为0

solve(am_a == Lm_a, m_a)        # solve m_a^{n+1}

solve(am_b == Lm_b, m_b)        # solve m_b^{n+1}

Eta_a = assemble(eta_a)
Eta_b = assemble(eta_b)
Eta_c = assemble(eta_c)

r = (Eta_a - Eta_c + r_old) / (1.0 - Eta_b)     #solve r^{n+1}

solve(aphi == Lphi, phi)         # combine phi

solve(aw == Lw, w)               # combine w

solve(autilde == Lutilde, utilde)   # combine utilde

solve(avarphi == Lvarphi, varphi)   # combine varphi
int.assign(varphi)
value3 = assemble(average)
varphi.vector()[:] = varphi.vector() - value3    # varphi^{n+1}在[0,1]上的积分为0

solve(am == Lm, m)               # combine m

solve(ap == Lp, p)               # solve p^{n+1}
int.assign(p)
value4 = assemble(average)
p.vector()[:] = p.vector() - value4   ## p^{n+1}在[0,1]上的积分为0

solve(au == Lu, u)           # solve u^{n+1}

print('r=', r)
print('S_n=', S_n)

phi_old.assign(phi)
w_old.assign(w)
u_old.assign(u)
p_old.assign(p)
varphi_old.assign(varphi)
m_old.assign(m)
r_old = r
t = t + dt

# Compute the error
phi_e.t = T
u_e.t = T
p_e.t = T
m_e.t = T
varphi_e.t = T
h_e.t = T

error_phi_L2 = errornorm(phi_e, phi, 'L2')
error_phi_H1 = errornorm(phi_e, phi, 'H1')
error_uL2 = errornorm(u_e, u, 'L2')
error_pL2 = errornorm(p_e, p, 'L2')
error_pH1 = errornorm(p_e, p, 'H1')
error_mL2 = errornorm(m_e, m, 'L2')
error_hL2 = (h_e - h)**2*dx
E_h = sqrt(abs(assemble(error_hL2)))
phie = interpolate(phi_e, P1Space)
error_r = abs(r - sqrt((1.0/(EPS*EPS))*assemble(intergal(phie)) + B))

print('error_phi_L2 =', error_phi_L2)
print('error_phi_H1 =', error_phi_H1)
print('error_uL2 =', error_uL2)
print('error_pL2 =', error_pL2)
print('error_pH1 = ', error_pH1)
print('error_mL2 =', error_mL2)
print('error_hL2 = ', E_h)
print('error_r = ', error_r)
print('T = ',  T, ',', 'M =', M, ',', 'Lambda =', Lambda, ',', 'K =', K, ',', 'niu_w = ', niu_w, ',', 'niu_f = ', niu_f, ',', 'miu = ', miu, ',', 'tau = ', tau, ',', 'beta = ', beta, ',', 'B = ', B, ',', 'k_f = ', k_f, ',', 'EPS = ', EPS)

# Polt
plt.figure()
plot(phi, title="phi")

phi_e_inp = interpolate(phi_e, P1Space)
plt.figure()
plot(phi_e_inp, title="phi_e")

plt.figure()
plot(p, title="p")

pe_inp = interpolate(p_e, P2Space)
plt.figure()
plot(pe_inp, title="pe")

plt.figure()
plot(u, title="u")

ue_inp = interpolate(u_e, P1_dgSpace)
plt.figure()
plot(ue_inp, title="ue")

plt.figure()
plot(m, title="m")

me_inp = interpolate(m_e, P1VectorSpace)
plt.figure()
plot(me_inp, title="me")

plt.figure()
plot(varphi, title="varphi")

varphie_inp = interpolate(varphi_e, P2Space)
plt.figure()
plot(varphie_inp, title="varphie")

plt.show()


In theory, the value of r and Sn should be similar, but after I run the program, the value of r and Sn is very different, and r is gradually decreasing. How to modify here ?
The following is the result when the running grid is divided into 1.0 / 128

r0=  1.4142135623730951

t=  0.0078125
r= 1.4142135623730951
S_n= 1.4142135623730951

t=  0.015625
r= 1.414187910672441
S_n= 1.4142058154844845

t=  0.0234375
r= 1.4141286188315734
S_n= 1.414179140332573

t=  0.03125
r= 1.414030554744989
S_n= 1.414131289019339

t=  0.0390625
r= 1.4138895806313179
S_n= 1.4140617500037005

t=  0.046875
r= 1.4137018768898324
S_n= 1.413970510188084

t=  0.0546875
r= 1.4134637320441845
S_n= 1.4138576606630635

t=  0.0625
r= 1.4131714905061583
S_n= 1.4137232960791615

t=  0.0703125
r= 1.4128215480407609
S_n= 1.4135674998836976

t=  0.078125
r= 1.4124103595575561
S_n= 1.4133903490259934

t=  0.0859375
r= 1.41193444907977
S_n= 1.413191920139061

t=  0.09375
r= 1.4113904195071811
S_n= 1.4129722937287472

t=  0.1015625
r= 1.410774961890038
S_n= 1.4127315565220093

t=  0.109375
r= 1.4100848643572623
S_n= 1.4124698026585278

t=  0.1171875
r= 1.4093170208441723
S_n= 1.4121871342497427

t=  0.125
r= 1.40846843969748
S_n= 1.4118836616154369

t=  0.1328125
r= 1.4075362521804875
S_n= 1.4115595033618282

t=  0.140625
r= 1.406517720867978
S_n= 1.4112147863828683

t=  0.1484375
r= 1.4054102479020651
S_n= 1.4108496458237658

t=  0.15625
r= 1.4042113830708816
S_n= 1.4104642250248425

t=  0.1640625
r= 1.4029188316673795
S_n= 1.4100586754540525

t=  0.171875
r= 1.4015304620834195
S_n= 1.4096331566318487

t=  0.1796875
r= 1.4000443130935962
S_n= 1.4091878360500254

t=  0.1875
r= 1.3984586007833693
S_n= 1.4087228890853771

t=  0.1953125
r= 1.396771725076741
S_n= 1.4082384989084076

t=  0.203125
r= 1.3949822758198587
S_n= 1.4077348563873093

t=  0.2109375
r= 1.3930890383784407
S_n= 1.407212159987287

t=  0.21875
r= 1.391090998708873
S_n= 1.4066706156653337

t=  0.2265625
r= 1.3889873478650303
S_n= 1.4061104367604265

t=  0.234375
r= 1.3867774859055253
S_n= 1.4055318438793059

t=  0.2421875
r= 1.384461025169007
S_n= 1.4049350647777654

t=  0.25
r= 1.382037792888353
S_n= 1.4043203342376596

t=  0.2578125
r= 1.379507833118138
S_n= 1.4036878939395057

t=  0.265625
r= 1.3768714079535125
S_n= 1.4030379923309333

t=  0.2734375
r= 1.3741289980226528
S_n= 1.4023708844909075

t=  0.28125
r= 1.3712813022391266
S_n= 1.401686831989864

t=  0.2890625
r= 1.3683292368049065
S_n= 1.400986102745791

t=  0.296875
r= 1.3652739334592163
S_n= 1.4002689708764173

t=  0.3046875
r= 1.362116736973009
S_n= 1.3995357165474651

t=  0.3125
r= 1.3588592018934875
S_n= 1.3987866258171575

t=  0.3203125
r= 1.355503088547701
S_n= 1.398021990477051

t=  0.328125
r= 1.3520503583188697
S_n= 1.3972421078892323

t=  0.3359375
r= 1.3485031682135893
S_n= 1.3964472808200579

t=  0.34375
r= 1.3448638647424989
S_n= 1.395637817270472

t=  0.3515625
r= 1.34113497714122
S_n= 1.3948140303030547

t=  0.359375
r= 1.337319209962443
S_n= 1.3939762378658813

t=  0.3671875
r= 1.3334194350738735
S_n= 1.3931247626133236

t=  0.375
r= 1.329438683100274
S_n= 1.3922599317239042

t=  0.3828125
r= 1.325380134351142
S_n= 1.391382076715329

t=  0.390625
r= 1.3212471092785079
S_n= 1.3904915332568084

t=  0.3984375
r= 1.317043058511857
S_n= 1.3895886409788183

t=  0.40625
r= 1.3127715525195578
S_n= 1.3886737432804

t=  0.4140625
r= 1.308436270947835
S_n= 1.3877471871341847

t=  0.421875
r= 1.3040409916899243
S_n= 1.3868093228892306

t=  0.4296875
r= 1.299589579739034
S_n= 1.3858605040718193

t=  0.4375
r= 1.2950859758793574
S_n= 1.384901087184382

t=  0.4453125
r= 1.290534185269625
S_n= 1.3839314315027063

t=  0.453125
r= 1.2859382659735514
S_n= 1.3829518988715321

t=  0.4609375
r= 1.2813023174908689
S_n= 1.3819628534987236

t=  0.46875
r= 1.2766304693418318
S_n= 1.3809646617481495

t=  0.4765625
r= 1.2719268697566168
S_n= 1.3799576919314775
t=  0.484375
r= 1.2671956745195159
S_n= 1.378942314098959

t=  0.4921875
r= 1.2624410360157536
S_n= 1.377918899829465

t=  0.5
r= 1.2576670925265407
S_n= 1.3768878220198562