How to use VectorFunctionSpace correctly?

variational formulation
1

I am trying to solve a time-dependent PDE (parabolic PDE) involving vector-valued functions on a 2D mesh.

I tried solving the problem using the VectorFunctionSpace functionality but it seems to be giving incorrect results. To test, I wrote two versions of the code for the 1-dimensional vector case.

In the first code, I used a “FunctionSpace” for a scalar function space on a 2D mesh.
In the second code, I used a “VectorFunctionSpace” of vector-dimension 1, on the 2D mesh.

The first code, (scalar function space) case gives the correct expected result. However, the second code which was expected to give the same result seems to be giving a different incorrect result. Could you please help me understand what might be going wrong? I would like to use the VectorFunctionSpace code to deal with vector-valued functions, but since it doesn’t even work for the 1d vector case, it seems to be giving incorrect results for the higher dimensional vector case as well.

The two codes written in python are as follows:

#Code 1 (Scalar function case, works correctly)
from dolfin import *
from mshr import *
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import ScalarFormatter, FormatStrFormatter

dx = dx(metadata={'quadrature_degree': 3})
Tf = 1e-3 # Final time
Nsteps = 5e4
dt = Tf/Nsteps 
t = 0. # initial time

Domain = Rectangle(Point(-20.,0.),Point(20.,10.))
Mesh = generate_mesh(Domain,60)
V = FunctionSpace(Mesh, 'P',2)
u = Function(V)
v = TestFunction(V)

phis = [ Expression('0.5*x[0]*x[0]',degree=2,domain=Mesh) ]
ic = Expression(str(1./40.),degree=2,domain=Mesh)
u0 = project(ic,V)
#PDE
J = Dx(u,0) + u*Dx(phis[0],0) 
ai_phi = (u*v+ dt*J*Dx(v,0))*dx 
Li_phi = u0*v*dx
F = ai_phi - Li_phi

#Time Loop for PDE solving 

Jac = derivative(F,u)
problem = NonlinearVariationalProblem(F, u, [],Jac)
solver  = NonlinearVariationalSolver(problem)

fig3 = plt.figure(3)
ax3 = fig3.add_subplot(111, projection='3d')
ax3.view_init(elev=45., azim=-65.)
plt.ion()
c2 = None
plotcount=0
for n in range(10):
        # Solve the problem
        solver.solve()
        ###### Update intial condition for next time step ######
        u0.assign(u)
        t=t+dt

        plotcount = plotcount + 1
        #region ############## plot the solution ####################
        if plotcount >=5:
                plt.figure(3)
                fig3.clear()
                ax3 = fig3.add_subplot(111, projection='3d')
                ax3.view_init(elev=45., azim=-65.)
                ax3.zaxis.set_major_formatter(ScalarFormatter(useOffset=False,useMathText=False))
                
                h2=plot(u,title=r'$\theta$')
                c2 = plt.colorbar(h2)
                c2.solids.set_edgecolor("face")
                
                
                plt.show(block=False)
                plotcount = 0
                print(n)
                print(t)

plt.savefig('u.png', bbox_inches='tight')
plt.show(block=True)

Output for code 1:
u

#Code 2 (Vector Function case, incorrect result)
from dolfin import *
from mshr import *
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import ScalarFormatter, FormatStrFormatter

dx = dx(metadata={'quadrature_degree': 3})
Tf = 1e-3 # Final time
Nsteps = 5e4
dt = Tf/Nsteps 
t = 0. # initial time

Domain = Rectangle(Point(-20.,0.),Point(20.,10.))
Mesh = generate_mesh(Domain,60)
Ndim = 1  # Number of dimensions
V = VectorFunctionSpace(Mesh, 'P',2,Ndim)  # Main difference from code 1
u = Function(V)
v = TestFunction(V)

phis = [ Expression('0.5*x[0]*x[0]',degree=2,domain=Mesh) ]
ic = Expression(tuple([str(1./40.)]*Ndim),degree=2,domain=Mesh)
u0 = project(ic,V)
#PDE
J = Dx(u[0],0) + u[0]*Dx(phis[0],0) 
ai_phi = (u[0]*v[0]+ dt*J*Dx(v[0],0))*dx 
Li_phi = u0[0]*v[0]*dx
F = ai_phi - Li_phi

#Time Loop for PDE solving 

Jac = derivative(F,u)
problem = NonlinearVariationalProblem(F, u, [],Jac)
solver  = NonlinearVariationalSolver(problem)

fig3 = plt.figure(3)
ax3 = fig3.add_subplot(111, projection='3d')
ax3.view_init(elev=45., azim=-65.)
plt.ion()
c2 = None
plotcount=0
for n in range(10):
        # Solve the problem
        solver.solve()
        ###### Update intial condition for next time step ######
        u0.assign(u)
        t=t+dt

        plotcount = plotcount + 1
        #region ############## plot the solution ####################
        if plotcount >=5:
                plt.figure(3)
                fig3.clear()
                ax3 = fig3.add_subplot(111, projection='3d')
                ax3.view_init(elev=45., azim=-65.)
                ax3.zaxis.set_major_formatter(ScalarFormatter(useOffset=False,useMathText=False))
                
                h2=plot(u[0],title=r'$\theta$')
                
                c2 = plt.colorbar(h2)
                c2.solids.set_edgecolor("face")
                
                
                plt.show(block=False)
                plotcount = 0
                print(n)
                print(t)

plt.savefig('u1.png', bbox_inches='tight')
plt.show(block=True)

Output for code 2 (incorrect result):
u1

Could you please help me understand why the use of VectorFunctionSpace leads to an incorrect result.

Thanks in advance.

1 Like
2

Apparently both the codes were correct. The discrepancy seen was due to approximation errors due to the coarse grids (both results with the coarse grid were incorrect and inconsistent). Making both grids finer and running the solver seems to give the same results. Best,
Sanket

2 Likes