Why using different Function spaces I get same solution with different scaling?

1

Hello.
I am implementing a code for the 2D wave equation. I have a Mixed formulation:my function space is the product of space for the pressure (u1) and a vector space for the 2D velocity (u2).
The code seem to work, however if I try different function space combination, such as P1-RT1, P2-RT2 or P1_-BDM1 I get same results but scaled differently. For example, if I use P1-RT1 after 0.4 seconds my simulation gives me a minimum pressure of -0.256, while if I change Lmix to P2-RT2 I get the same pressure distribution but the minimum is -0.64.
Does anyone can explain me why? Shoud not get the same solution?
Here is my code:

from fenics import *
    import dolfin
    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D 
    from matplotlib import cm
    from matplotlib.ticker import LinearLocator, FormatStrFormatter
    from matplotlib import rc

# Constants 
dt = 0.0005 # time step size
FT = 0.5    # final time
t = 0       # initialize time

# Mesh
N = 20
M = 20
mesh = RectangleMesh(Point(0,0), Point(1,1), N, M, "left")
n = FacetNormal(mesh)
# plot(mesh)

# Finite element spaces
P1 = FiniteElement("P-", triangle, 1, 0)
P2 = FiniteElement("P-", triangle, 2, 0)
P3 = FiniteElement("P-", triangle, 3, 0)

RT1 = FiniteElement("P-", triangle, 1, 1)
RT2 = FiniteElement("P-", triangle, 2, 1)
RT3 = FiniteElement("P-", triangle, 3, 1)

P1_ = FiniteElement("P", triangle, 1, 0)
P2_ = FiniteElement("P", triangle, 2, 0)
P3_ = FiniteElement("P", triangle, 3, 0)

BDM1 = FiniteElement("P", triangle, 1, 1)
BDM2 = FiniteElement("P", triangle, 2, 1)
BDM3 = FiniteElement("P", triangle, 3, 1)

# Define mixed Functional Space
Lmix = FunctionSpace(mesh, P1*RT1)

# Trial functions
u1, u2 = TrialFunctions(Lmix)

# Test functions
v1, v2 = TestFunctions(Lmix)

# Functions - previous step - boundary inputs
u_n = Function(Lmix) # previous step (respect to n+1)
u_b = Function(Lmix) # boundary inputs
u_s = Function(Lmix) # solution (to be updated at each time step in the iteration loop)
# Compute the solution at each time step

while(t < 0.4):
  
    #1 Update "gaussian" boundary conditions in velocity on Omega left
    u_D = Expression(("0", "pow(sin(pi*t), 2)*near(x[0],0)*(t<1)", "0"), degree = 2,t = t)  
    u_b = project(u_D, Lmix)
    
    # Variational problem for harmonic bc
    F = v1*u1*dx - v1*u_n.split()[0]*dx + dot(v2, u2)*dx - dot(v2, u_n.split()[1])*dx - dt*dot(grad(v1), u2)*dx + \
    	dt*dot(v2, grad(u_n.split()[0]))*dx - dt*v1*dot(u_b.split()[1], n)*ds
    
    a, L = lhs(F), rhs(F)
    solve(a == L, u_s)
                                             
    # Assign solution to previous step vector
    u_n.assign(u_s)
    print(t)
    
    # Plot u_1 (pressure)
    p = plot(u_s.split()[0])
    plt.colorbar(p)
    plt.show()
    
    # Update current time
    t += dt

verticesD = {}

for cell in cells(mesh):
    local_to_global_map0 = Lmix.sub(0).dofmap().cell_dofs(cell.index())
    index = 0
    for vertex in vertices(cell):
        verticesD[tuple(vertex.point().array()[0:2])] = local_to_global_map0[index]
        index += 1


xx = np.linspace(0,1,M+1)
yy = np.linspace(0,1,N+1)

XX, YY = np.meshgrid(xx,yy)
vector = u_s.split()[0].vector()[:]
Z = np.zeros((N+1,M+1))
for i in range (0, N+1):
    for j in range(0,M+1):
        Z[j][i] =  vector[verticesD.get((i*1/N, j*1/M), "none")]
        
fig = plt.figure()
ax = fig.gca(projection='3d')
        
c = ax.plot_surface(XX, YY, Z, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)

ax.view_init(60, 250)

fig.colorbar(c, shrink=0.5, aspect=5)
plt.show()

Moreover, in my code I implemented an external harmonic boundary input sin^2(pi*t) that act on the left boundary (x=0). Is there an intelligent way to implement such an input only on point, the origin for example? Because if I put in the code something like

  u_D = Expression(("0", "pow(sin(pi*t), 2)*(near(x[0],0) && near(x[1],0))*(t<1)", "0"), degree = 2,t = t)  
  u_b = project(u_D, Lmix)

it does not work. I tried also PointSource function in the loop but I do not get proper results.

a, L = lhs(F), rhs(F)
A, b = assemble_system(a, L)

# create a DirichletBC
delta = PointSource(Lmix, Point(0, 0), pow(sin(pi*t), 2))
delta.apply(b)

# Define solver (iterative)
solver = LUSolver(A,method='mumps')
solver.parameters['symmetric'] = True
solver.parameters['report'] = True
solver.parameters['verbose'] = True

solver.solve(u_s.vector(),b)