Issue with the get_local. max() function

General
1

Hey guys. I have an issue when using get_local.max() function

So first for my code, I am using a 1D interval

mesh = IntervalMesh(1000, 0, 1)

And after solving the problem. I want to find the max value and plot it, so I use

print(theta_list[999].vector().get_local().max())
plt.ylim(0, 3.14)
plot(theta_list[999])
plt.show()

To my surprise. The local max output is 53.64738743162838 but the plot clearly shows that the max on the interval (0, 1) is below 2. the plot looks like this

image
image635×560 21 KB

Any help and suggestions are appreciated!

2

Some more details. I add local_range and max_index. It seems they are not consistent with the graph.

The codes are

max_index = np.argmax(theta_list[999].vector()[:])
x= ME.tabulate_dof_coordinates()
print(x[max_index])
print(theta_list[999].vector().get_local().max())
print(theta_list[999].vector().local_range())
plt.ylim(0, 3.14)
plot(theta_list[999])
plt.show()

outputs are

[0.499]
53.64669017236389
(0, 2002)

but the plot looks like this

image
image635×560 21.5 KB

In case you need the whole codes are

from dolfin import *
import matplotlib.pyplot as plt
from tqdm import tqdm
import random
import numpy as np

# Define 1D mesh and space
set_log_level(40)
mesh = IntervalMesh(1000, 0, 1)
Element1 = FiniteElement("CG", mesh.ufl_cell(), 1)
Element2 = FiniteElement("CG", mesh.ufl_cell(), 1)
W_elem = MixedElement([Element1, Element2])
ME = FunctionSpace(mesh, W_elem)

# Define functions
u_old = Function(ME)
u_new = Function(ME)
theta_old, vy_old = split(u_old)
theta_new, vy_new = split(u_new)
vy_mid = 0.5 * vy_old + 0.5 * vy_new
theta_mid = 0.5 * theta_old + 0.5 * theta_new
tf = TestFunction(ME)
q, v = split(tf)
kapa_ga = 1
lamb = 1.5
f1 = lamb * cos(2 * theta_mid)
alpha = -450
f2 = alpha * cos(2 * theta_new)
dt = 0.01


# Boundary condition
def boundary(x, on_boundary):
    return on_boundary


bc = DirichletBC(ME, [Constant(1.57), Constant(0)], boundary)


# Initial Condition
class InitialConditions(UserExpression):
    def eval(self, values, x):
        values[0] = 1.57+0.01*random.gauss(0, 1)
        values[1] = 0

    def value_shape(self):
        return (2,)


# Set initial conditions to old function
u_init = InitialConditions()
u_old.interpolate(u_init)

# Weak statement of the equations
F0 = (theta_new - theta_old) / dt * q * dx + kapa_ga*theta_mid.dx(0) * q.dx(0) * dx - 0.5 * (f1 + 1) * vy_mid.dx(0) * q * dx
F1 = (-vy_new.dx(0) * v.dx(0)) * dx + f2 * theta_new.dx(0) * v * dx
F = F0 + F1

t = 0.0
T = 1000 * dt
theta_list = []
vy_list = []

for t in tqdm(np.arange(0, 10.01, 0.01)):
    solve(F == 0, u_new, bc)
    u_old.vector()[:] = u_new.vector()
    u_copy = Function(ME)
    u_copy.vector()[:] = u_new.vector()
    theta_list.append(u_copy.split()[0])
    vy_list.append(u_copy.split()[1])

max_index = np.argmax(theta_list[999].vector()[:])
x= ME.tabulate_dof_coordinates()
print(x[max_index])
print(theta_list[999].vector().get_local().max())
print(theta_list[999].vector().local_range())
plt.ylim(0, 3.14)
plot(theta_list[999])
plt.show()