Offline Phase: SGREEDY algorithm for sensor placement

This notebook implements the following algorithm for sensor placement:

• Stabilised GREEDY algorithm (SGREEDY)

import numpy as np
import os
from IPython.display import clear_output
import pickle

import gmsh
from mpi4py import MPI
from dolfinx.io.gmshio import model_to_mesh
from dolfinx.fem import FunctionSpace

from pyforce.tools.write_read import ImportH5, StoreFunctionsList as store
from pyforce.tools.functions_list import FunctionsList
from pyforce.offline.sensors import SGREEDY

path_off ='./Offline_results/'

The geometry is imported from “ANL11A2_octave.geo”, generated with GMSH. Then, the mesh is created with the gmsh module.

gdim = 2

model_rank = 0
mesh_comm = MPI.COMM_WORLD

# Initialize the gmsh module
gmsh.initialize()

# Load the .geo file
gmsh.merge('ANL11A2_octave.geo')
gmsh.model.geo.synchronize()

# Set algorithm (adaptive = 1, Frontal-Delaunay = 6)
gmsh.option.setNumber("Mesh.Algorithm", 6)

gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")

# Domain
domain, ct, ft = model_to_mesh(gmsh.model, comm = mesh_comm, rank = model_rank, gdim = gdim )
gmsh.finalize()

fuel1_marker    = 1
fuel2_marker    = 2
fuel_rod_marker = 3
refl_marker     = 4

void_marker     = 10
sym_marker      = 20

tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)

clear_output()

Importing Snapshots

The snapshots are loaded and stored into suitable data structures.

# Defining the functional space
V = FunctionSpace(domain, ("Lagrange", 1))

# Define the variables to load
var_names = [
             'phi_1',
             'phi_2'
             ]

tex_var_names = [
                 r'\phi_1',
                 r'\phi_2'
                 ]

# Snapshot path
path_FOM = './Snapshots/'

################ Importing Snapshots ########################
test_snaps = list()

for field_i in range(len(var_names)):
    test_snaps.append(FunctionsList(V))

    tmp_FOM_list, _ = ImportH5(V, path_FOM+'test_snap_'+var_names[field_i], var_names[field_i])

    for mu in range(len(tmp_FOM_list)):
        test_snaps[field_i].append(tmp_FOM_list(mu))

    del tmp_FOM_list

test_params = list()
for field_i in range(len(var_names)):
    with open(path_FOM+'./test.param', 'rb') as f:
        test_params.append(pickle.load(f))

Let us import the POD modes

bf = dict()

for field_i in range(len(var_names)):

    bf[var_names[field_i]] = ImportH5(V, path_off+'/BasisFunctions/basisPOD_'+var_names[field_i], 'POD_'+var_names[field_i])[0]

SGREEDY for POD basis

A greedy procedure is set up to maximise the inf-sup constant \(\beta_{N,M}\), for 3 different values if the point spread

s = [0.1, 1., 2.5]
sgreedy_pod = dict()

for field_i in range(len(var_names)):
    sgreedy_pod[var_names[field_i]] = [SGREEDY(domain, bf[var_names[field_i]], V, var_names[field_i], s_jj) for s_jj in s]

Since the snapshots belong to \(H^1\subset L^2\), the Riesz representation of a functional is sought in this space, endowed with the inner product \begin{equation*} \left(u,v\right)_{H^1}=\int_\Omega \nabla u\cdot \nabla v\,d\Omega + \int_\Omega u\cdot v\,d\Omega \end{equation*} The results will be compared with the approximation properties of the case with the Riesz representation in \(L^2\), whose inner product is \begin{equation*} \left(u,v\right)_{L^2} = \int_\Omega u\cdot v\,d\Omega \end{equation*}

Nmax = 10
Mmax = 30

sam_every = 2

is_H1 = [False, True]
fun_space_label = ['L2', 'H1']

for field_i in range(len(var_names)):
    for jj in range(len(s)):
        for kk in range(len(is_H1)):

            print('SGREEDY for '+var_names[field_i]+' with s={:.2f}'.format(s[jj])+' and Riesz representation in '+fun_space_label[kk])

            sgreedy_pod[var_names[field_i]][jj].generate(Nmax, Mmax, tol = 0.5, sampleEvery = sam_every, verbose=False, is_H1 = is_H1[kk])
            store(domain, sgreedy_pod[var_names[field_i]][jj].basis_sens,
                  'SGREEDYPOD_' +var_names[field_i]+'_s_{:.2e}'.format(s[jj]),
                  path_off+'/BasisSensors/sensorsSGREEDYPOD_' + var_names[field_i]+'_s_{:.2e}_'.format(s[jj])+fun_space_label[kk])

Generating sensors (sampled every 2 cells): 3435.000 / 3435.00 - 0.002 s/it
m = 30, n = 10 | beta_n,m = 0.084779
Generating sensors (sampled every 2 cells): 3435.000 / 3435.00 - 0.003 s/it
m = 30, n = 10 | beta_n,m = 0.146526
Generating sensors (sampled every 2 cells): 3435.000 / 3435.00 - 0.002 s/it
m = 30, n = 10 | beta_n,m = 0.180025
Generating sensors (sampled every 2 cells): 3435.000 / 3435.00 - 0.003 s/it
m = 30, n = 10 | beta_n,m = 0.223101
Generating sensors (sampled every 2 cells): 3435.000 / 3435.00 - 0.002 s/it
m = 30, n = 10 | beta_n,m = 0.393625
Generating sensors (sampled every 2 cells): 3435.000 / 3435.00 - 0.003 s/it
m = 30, n = 10 | beta_n,m = 0.398683
Generating sensors (sampled every 2 cells): 3435.000 / 3435.00 - 0.002 s/it
m = 30, n = 10 | beta_n,m = 0.085134
Generating sensors (sampled every 2 cells): 3435.000 / 3435.00 - 0.003 s/it
m = 30, n = 10 | beta_n,m = 0.148458
Generating sensors (sampled every 2 cells): 3435.000 / 3435.00 - 0.002 s/it
m = 30, n = 10 | beta_n,m = 0.173819
Generating sensors (sampled every 2 cells): 3435.000 / 3435.00 - 0.003 s/it
m = 30, n = 10 | beta_n,m = 0.214982
Generating sensors (sampled every 2 cells): 3435.000 / 3435.00 - 0.002 s/it
m = 30, n = 10 | beta_n,m = 0.357794
Generating sensors (sampled every 2 cells): 3435.000 / 3435.00 - 0.003 s/it
m = 30, n = 10 | beta_n,m = 0.371680