https://github.com/hpc-maths/GenEO
Tip revision: 8cfb272518a47a1ea20715e6c0e5182632e7b591 authored by Loic Gouarin on 21 May 2024, 17:47:31 UTC
Create AUTHORS.txt
Create AUTHORS.txt
Tip revision: 8cfb272
post_treatment.py
import pyvista as pv
import numpy as np
import h5py
import os
import re
import json
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from io import BytesIO
import pandas as pd
pv.set_plot_theme("document")
def get_cell_value(y, E1, E2, stripe_nb,Ly):
stripe_nb = 2
if stripe_nb == 0:
return E2
elif stripe_nb == 1:
if 1./7 <= y - np.floor(y) <= 2./7:
return E1
else:
return E2
elif stripe_nb == 2:
if (1./7 <= y - np.floor(y) <= 2./7) or (3./7 <= y - np.floor(y) <= 4./7):
return E1
else:
return E2
elif stripe_nb == 3:
if (1./7 <= y - np.floor(y) <= 2./7) or (3./7 <= y - np.floor(y) <= 4./7) or (5./7 <= y - np.floor(y) <= 6./7):
return E1
else:
return E2
def plot_solution(path, data, E1, E2, stripe_nb, Ly):
scale = 100000
p = pv.Plotter(off_screen=True)
filename = os.path.join(path, 'solution_2d.vts')
if os.path.exists(filename):
grid = pv.read(os.path.join(path, 'solution_2d.vts'))
E = np.zeros((grid.n_cells))
for ic in range(grid.n_cells):
cell = grid.cell_points(ic)
center_y = .5*(cell[0, 1] + cell[2, 1])
E[ic] = get_cell_value(center_y, E1, E2, stripe_nb,Ly)
grid.cell_arrays.update({'E': E})
disp = grid.point_arrays['Unnamed']
print(grid.x.shape, disp.shape)
grid.x.T.flat[:] += scale*disp[:, 0]
grid.y.T.flat[:] += scale*disp[:, 1]
p.add_mesh(grid, scalars='E', show_edges=True, show_scalar_bar=False, n_colors=2, edge_color='#464646', cmap=['#dbd9d3', '#a2acbd'])
print(os.path.join(path, 'solution.png'))
p.show(cpos='xy', screenshot=os.path.join(path, 'solution.png'))
def plot_coarse_vec(path, data, E1, E2, stripe_nb,Ly):
scale = 0.5
for k, v in data.items():
print(k, v)
ncol = int(np.ceil(np.sqrt(len(v))))
p = pv.Plotter(shape=(ncol, ncol), off_screen=True)
with open(os.path.join(path, f'properties_{k}.txt')) as json_file:
prop = json.load(json_file)
iplot = 0
for i in range(ncol):
for j in range(ncol):
p_individual = pv.Plotter(off_screen=True)
if iplot < len(v):
h5 = h5py.File(os.path.join(path, v[iplot]))
coords = h5['coordinates']
coarse_vec = h5['coarse_vec']
nx = 0
dim = len(coords[:].shape)
if dim == 1:
coords = coords[:].copy()
coords.shape = (coords.size//2, 2)
xx = coords[:, 0]
while xx[nx] < xx[nx + 1]:
nx += 1
nx += 1
ny = coords.size//2//nx
x = np.zeros((ny, nx))
x.flat = coords[:, 0] #+ scale*coarse_vec[::2]
y = np.zeros((ny, nx))
y.flat = coords[:, 1] #+ scale*coarse_vec[1::2]
scalar = np.zeros((ny, nx))
grid = pv.StructuredGrid(x, y, np.zeros((ny, nx)))
E = np.zeros((grid.n_cells))
for ic in range(grid.n_cells):
cell = grid.cell_points(ic)
center_y = .5*(cell[0, 1] + cell[1, 1])
E[ic] = get_cell_value(center_y, E1, E2, stripe_nb,Ly)
grid.cell_arrays.update({'E': E})
if dim == 1:
coarse_vec = coarse_vec[:].copy()
coarse_vec.shape = (coarse_vec.size//2, 2)
grid.x[:] += scale*np.reshape(coarse_vec[:, 0], (ny, nx, 1))
grid.y[:] += scale*np.reshape(coarse_vec[:, 1], (ny, nx, 1))
p.subplot(i, j)
p.add_text(f"{prop['eigs'][iplot]}", position='upper_edge', font_size=4)
p.add_mesh(grid, show_edges=True, show_scalar_bar=False, n_colors=2, edge_color='#464646', cmap=['#dbd9d3', '#a2acbd'])
p_individual.add_mesh(grid, show_edges=True, show_scalar_bar=False, n_colors=2, edge_color='#464646', cmap=['#dbd9d3', '#a2acbd'])
p_individual.show(cpos='xy', window_size=[800, 800], screenshot=os.path.join(path, f'coarse_vec_{k}_{iplot}.png'))
iplot += 1
p.show(cpos='xy', screenshot=os.path.join(path, f'coarse_vec_{k}.png'))
def plot_eigenvalues(path, eigs):
fig, ax = plt.subplots()
for i, eig in enumerate(eigs):
ax.scatter(i*np.ones_like(eig), eig)
ax.set_xlabel('Subdomain number')
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
ax.set_ylabel('Generalized eigenvalues (log scale)')
ax.set_yscale('log')
fig.savefig(os.path.join(path, 'eigenvalues.png'), dpi=300)
def plot_condition_number(path, condition):
fig, ax = plt.subplots()
for k, v in condition.items():
data = np.asarray(v)
index = np.argsort(data[:, 0])
data = data[index]
ax.plot(data[:, 1], data[:, 2], label=k, marker='.')
ax.set_xlabel('Coarse space size')
ax.set_ylabel('Condition number (log scale)')
ax.set_yscale('log')
ax.legend()
fig.savefig(os.path.join(path, 'condition_number.png'), dpi=300)
regexp = re.compile('coarse_vec_(\d+)_(\d+).h5')
regexp_case_9 = re.compile('case_9_*')
regexp_case_8 = re.compile('case_8_*')
regexp_case_7 = re.compile('case_7_*')
regexp_case_6 = re.compile('case_6_*')
regexp_case_5 = re.compile('case_5_*')
regexp_case_4 = re.compile('case_4_*')
regexp_case_2 = re.compile('case_2_*')
regexp_case_3 = re.compile('case_3_*')
path = 'output.d'
fig, ax = plt.subplots()
condition = {
'classical GenEO': [],
'AWG': [],
# 'pcnew_neg': []
}
dfs = []
for (dirpath, dirnames, filenames) in os.walk(path):
if 'results.json' in filenames:
data = {}
for f in filenames:
res = regexp.search(f)
if res:
i, rank = res.groups()
d = data.get(rank, [])
d.append(f)
data[rank] = d
with open(os.path.join(dirpath, 'results.json')) as json_file:
results = json.load(json_file)
res_case_4 = regexp_case_4.search(os.path.split(dirpath)[-1])
if res_case_4:
name = open(os.path.join(dirpath, 'name.txt')).read()
data_case4 = {'name': name,
'kappa': results['kappa'][0],
'labdamin': results['lambdamin'][0],
'lambdamax': results['lambdamax'][0],
'V0dim': int(results['V0dim']),
# 'vneg': results['vneg'],
'sum_gathered_nneg': int(results['sum_gathered_nneg']),
}
dfs.append(data_case4)
# plot_coarse_vec(dirpath, data, results['E1'], results['E2'], results['stripe_nb'],results['Ly'])
# plot_solution(dirpath, data, results['E1'], results['E2'], results['stripe_nb'],results['Ly'])
# plot_eigenvalues(dirpath, results['GenEOV0_gathered_Lambdasharp'])
#####CASE 2
# res_case_2 = regexp_case_2.search(os.path.split(dirpath)[-1])
# if res_case_2:
# #name = open(os.path.join(dirpath, 'name.txt')).read()
# data = {'nu': results['nu1'],
# 'kappa': results['kappa'][0],
# 'iter': len(results['precresiduals']),
# 'lambdamin': results['lambdamin'][0],
# 'lambdamax': results['lambdamax'][0],
# 'V0dim': int(results['V0dim']),
# 'sum_gathered_nneg': int(results['sum_gathered_nneg']),
# }
# dfs.append(data)
# #plot_coarse_vec(dirpath, data, results['E1'], results['E2'], results['stripe_nb'],results['Ly'])
# plot_eigenvalues(dirpath, results['GenEOV0_gathered_Lambdasharp'])
#####END CASE 2
#####CASE 3
#res_case_3 = regexp_case_3.search(os.path.split(dirpath)[-1])
#if res_case_3:
# name = open(os.path.join(dirpath, 'name.txt')).read()
# data = {' ': name,
# 'E_1': results['E1'],
# 'E_2': results['E2'],
# 'kappa': results['kappa'][0],
# 'iter': len(results['precresiduals']),
# 'lambdamin': results['lambdamin'][0],
# 'lambdamax': results['lambdamax'][0],
# 'V0dim': int(results['V0dim']),
# 'sum_gathered_nneg': int(results['sum_gathered_nneg']),
# }
# dfs.append(data)
# #plot_coarse_vec(dirpath, data, results['E1'], results['E2'], results['stripe_nb'],results['Ly'])
# plot_eigenvalues(dirpath, results['GenEOV0_gathered_Lambdasharp'])
#####END CASE 3
###### I deleted CASE 4 by accident, it is in Loic's commit from June 10th 2021
#######CASE 1
# if os.path.basename(dirpath).split('_')[1] == '1':
# if os.path.basename(dirpath).split('_')[2] == 'pcbnn':
# condition['classical GenEO'].append((results['taueigmax'], results['V0dim'], results['kappa'][0]))
# else:
# condition['AWG'].append((results['taueigmax'], results['V0dim'], results['kappa'][0]))
######## END CASE 1
#######CASE 9
# if os.path.basename(dirpath).split('_')[1] == '9':
# if os.path.basename(dirpath).split('_')[2] == 'pcbnn':
# condition['classical GenEO'].append((results['taueigmax'], results['V0dim'], results['kappa'][0]))
# else:
# condition['AWG'].append((results['taueigmax'], results['V0dim'], results['kappa'][0]))
######## END CASE 9
#####CASE 5 and 8
# res_case_5 = regexp_case_5.search(os.path.split(dirpath)[-1])
# res_case_8 = regexp_case_8.search(os.path.split(dirpath)[-1])
# if res_case_5 or res_case_8:
# name = open(os.path.join(dirpath, 'name.txt')).read()
# data = {' ': name,
# 'nu': results['nu1'],
# 'rtol': results['Aposrtol'],
# 'kappa': results['kappa'][0],
# 'iter': len(results['precresiduals']),
# 'lambdamin': results['lambdamin'][0],
# 'lambdamax': results['lambdamax'][0],
# 'V0dim': int(results['V0dim']),
# 'sum_gathered_nneg': int(results['sum_gathered_nneg']),
# }
# dfs.append(data)
# #plot_coarse_vec(dirpath, data, results['E1'], results['E2'], results['stripe_nb'],results['Ly'])
# plot_eigenvalues(dirpath, results['GenEOV0_gathered_Lambdasharp'])
#####END CASE 5 and 8
# ax.plot(results['precresiduals'])
######CASE 6
res_case_6 = regexp_case_6.search(os.path.split(dirpath)[-1])
if res_case_6:
name = open(os.path.join(dirpath, 'name.txt')).read()
data = {' ': name,
'nbstripes': results['stripe_nb'],
'kappa': results['kappa'][0],
'iter': len(results['precresiduals']),
'lambdamin': results['lambdamin'][0],
'lambdamax': results['lambdamax'][0],
'V0dim': int(results['V0dim']),
'sum_gathered_nneg': int(results['sum_gathered_nneg']),
}
dfs.append(data)
#plot_coarse_vec(dirpath, data, results['E1'], results['E2'], results['stripe_nb'],results['Ly'])
# plot_eigenvalues(dirpath, results['GenEOV0_gathered_Lambdasharp'])
#######END CASE 6
#####CASE 7
res_case_7 = regexp_case_7.search(os.path.split(dirpath)[-1])
if res_case_7:
name = open(os.path.join(dirpath, 'name.txt')).read()
data = {' ': name,
'N': len(results['minV0_gathered_dim']),
'kappa': results['kappa'][0],
'iter': len(results['precresiduals']),
'lambdamin': results['lambdamin'][0],
'lambdamax': results['lambdamax'][0],
'V0dim': int(results['V0dim']),
'sum_gathered_nneg': int(results['sum_gathered_nneg']),
}
dfs.append(data)
#plot_coarse_vec(dirpath, data, results['E1'], results['E2'], results['stripe_nb'],results['Ly'])
plot_eigenvalues(dirpath, results['GenEOV0_gathered_Lambdasharp'])
#####END CASE 7
# ax.set_xlabel('iteration number')
# ax.set_ylabel('residual')
# ax.set_yscale('log')
# fig.savefig(os.path.join(path, 'residuals.png'), dpi=300)
#plot_condition_number(path, condition) #CASE 1 and 9
df = pd.DataFrame(dfs)
print(df.to_latex())