2D case - Figure 12: Space variability training#

Two-parameter 2D problem with heterogeneous stiffness - Investigation of the convergence of the reduced-order model and of the evolution of the size of the reduced-order basis.

Libraries import#

import sys  
import torch
import torch.nn as nn
from neurom.HiDeNN_PDE import MeshNN, NeuROM, MeshNN_2D, MeshNN_1D
import neurom.src.Pre_processing as pre
from neurom.src.PDE_Library import Strain, Stress,VonMises_plain_strain
from neurom.src.Training import Training_NeuROM_multi_level
import neurom.Post.Plots as Pplot
import time
import os
import torch._dynamo as dynamo
from importlib import reload  
import tomllib
import numpy as np
import argparse

torch.manual_seed(0)

* WARNING: could not load tikzplotlib

<torch._C.Generator at 0x108df9670>

Load the config file#

    Configuration_file = 'Configurations/config_2D_ROM_SV.toml'

    with open(Configuration_file, mode="rb") as file:
        config = tomllib.load(file)

Definition of the space domain and mechanical proprieties of the structure#

The initial Material parameters, the geometry, mesh and the boundary conditions are set.

# Overwrites the multi-level training
config["training"]["multiscl_max_refinment"] = 3
# MaxElemSize2D = config["interpolation"]["MaxElemSize2D"] = 0.5

# Material parameters definition

Mat = pre.Material(             flag_lame = False,                                  # If True should input lmbda and mu instead of E and nu
                                coef1     = config["material"]["E"],                # Young Modulus
                                coef2     = config["material"]["nu"]                # Poisson's ratio
                )


# Create mesh object
MaxElemSize = pre.ElementSize(
                                dimension     = config["interpolation"]["dimension"],
                                L             = config["geometry"]["L"],
                                order         = config["interpolation"]["order"],
                                np            = config["interpolation"]["np"],
                                MaxElemSize2D = config["interpolation"]["MaxElemSize2D"]
                            )
Excluded = []
Mesh_object = pre.Mesh( 
                                config["geometry"]["Name"],                         # Create the mesh object
                                MaxElemSize, 
                                config["interpolation"]["order"], 
                                config["interpolation"]["dimension"]
                        )

Mesh_object.AddBorders(config["Borders"]["Borders"])
Mesh_object.AddBCs(                                                                 # Include Boundary physical domains infos (BCs+volume)
                                config["geometry"]["Volume_element"],
                                Excluded,
                                config["DirichletDictionryList"]
                    )                   

Mesh_object.MeshGeo()                                                               # Mesh the .geo file if .msh does not exist
Mesh_object.ReadMesh()       
Mesh_object.ExportMeshVtk()

Parametric study definition#

The hypercube describing the parametric domain used for the tensor decomposition is set-up here

ParameterHypercube = torch.tensor([ [   config["parameters"]["para_1_min"],
                                        config["parameters"]["para_1_max"],
                                        config["parameters"]["N_para_1"]],
                                    [   config["parameters"]["para_2_min"],
                                        config["parameters"]["para_2_max"],
                                        config["parameters"]["N_para_2"]]])

Initialisation of the surrogate model#

ROM_model = NeuROM(                                                                 # Build the surrogate (reduced-order) model
                    Mesh_object, 
                    ParameterHypercube, 
                    config,
                    config["solver"]["n_modes_ini"],
                    config["solver"]["n_modes_max"]
                    )

Training the model#

ROM_model.Freeze_Mesh()                                                             # Set space mesh coordinates as untrainable
ROM_model.Freeze_MeshPara()                                                         # Set parameters mesh coordinates as untrainable

ROM_model.TrainingParameters(   
                                loss_decrease_c = config["training"]["loss_decrease_c"], 
                                Max_epochs = config["training"]["n_epochs"], 
                                learning_rate = config["training"]["learning_rate"]
                            )

ROM_model.train()                                                                   # Put the model in training mode
ROM_model = Training_NeuROM_multi_level(ROM_model,config, Mat)         
# ROM_model,Mesh_object = Training_NeuROM_multi_level(ROM_model,config, Mat)         

ROM_model = ROM_model[0]

Plotting area#

Reproducing figure 8, a-b-c

tikz = False

import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
import matplotlib
plt.rcParams['text.usetex'] = False

Modes_flag = ROM_model.training_recap["Mode_vect"]
error = ROM_model.training_recap["Loss_vect"]
decay = ROM_model.training_recap["Loss_decrease_vect"]
threshold = config["training"]["loss_decrease_c"]

name = 'FigXa'


# plot Fig 8a

fig = plt.figure()
ax = fig.add_subplot(111)

## First curve
ax.invert_yaxis()
g1 = ax.semilogy(-torch.tensor(error), color='#d95319ff')
ax.set_ylabel(r'$ - J\left(u\left(x\right)\right)$',color='#d95319ff')
ax.tick_params(axis='y', colors='#d95319ff', which='both')
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
ax.set_xlabel(r'Epochs')

## Second curve
ax2 = ax.twinx()
g2 = ax2.plot(Modes_flag, color='#247ab5ff')
ax2.set_ylabel(r'$m$',color='#247ab5ff')
ax2.tick_params(axis='y', colors='#247ab5ff')
ax2.yaxis.set_major_locator(MaxNLocator(integer=True))

if tikz:
    import tikzplotlib
    tikzplotlib.save('Results/'+name+'_zoom.tex')
plt.savefig('Results/'+name+'_zoom.pdf', transparent=True, bbox_inches = "tight")

plt.show() 
plt.clf() 





name = 'FigXb'

# pre processing for padding and zooming on epochs after rough training during the first epochs
Zoom_depth = np.min(np.where(np.array(Modes_flag) == np.array(Modes_flag)[0]+1))
Zoom_start_index = int(np.floor(0.9*Zoom_depth))
second_stages_epochs = len(error) - len(Modes_flag)
Modes_flag.extend([Modes_flag[-1]]*second_stages_epochs)
x_indexes = np.arange(len(Modes_flag[Zoom_start_index:]))+Zoom_start_index

# plot Fig 8b

fig = plt.figure()
ax = fig.add_subplot(111)

## First curve
ax.invert_yaxis()
g1 = ax.semilogy(x_indexes,-torch.tensor(error[Zoom_start_index:]), color='#d95319ff')
ax.set_ylabel(r'$ - J\left(u\left(x\right)\right)$',color='#d95319ff')
ax.tick_params(axis='y', colors='#d95319ff', which='both')
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
ax.set_xlabel(r'Epochs')

## Second curve
ax2 = ax.twinx()
g2 = ax2.plot(x_indexes,Modes_flag[Zoom_start_index:], color='#247ab5ff')
ax2.set_ylabel(r'$m$',color='#247ab5ff')
ax2.tick_params(axis='y', colors='#247ab5ff')
ax2.yaxis.set_major_locator(MaxNLocator(integer=True))

if tikz:
    import tikzplotlib
    tikzplotlib.save('Results/'+name+'_zoom.tex')
plt.savefig('Results/'+name+'_zoom.pdf', transparent=True, bbox_inches = "tight")

plt.show() 
plt.clf() 


name = 'FigXc'

# plot Fig 8c

ax = plt.gca()
ax.semilogy(decay,color='#d95319ff')
ax.tick_params(axis='y', colors='#d95319ff')
ax.set_ylabel(r'd log($J\left(u\left(x\right)\right)$)',color='#d95319ff')
plt.axhline(threshold,color = 'k')
ax.xaxis.set_major_locator(MaxNLocator(integer=True))


# plt.ylim((0.01,20))
ax2 = plt.gca().twinx()
ax2.plot(Modes_flag,color='#247ab5ff')
ax2.set_ylabel(r'$m$',color='#247ab5ff')
ax2.tick_params(axis='y', colors='#247ab5ff')
ax2.yaxis.set_major_locator(MaxNLocator(integer=True))

if tikz:
    import tikzplotlib
    tikzplotlib.save('Results/'+name+'.tex')
plt.savefig('Results/'+name+'.pdf', transparent=True, bbox_inches = "tight")
plt.show() 

plt.clf()

ROM_model = ROM_model[0]

Export plot data#

The data are saved in a csv file so that they can be plotted in the article using pgfplot.

import pandas as pd
epochs = list(range(len(Modes_flag)))
N = len(decay)
Name = "Fig_SV"

df_full = pd.DataFrame(np.stack((epochs[:N],decay[:N],Modes_flag[:N],(-torch.tensor(error)[:N]).tolist()) ,axis=1), columns=['epochs', 'Decay','Modes',"Loss"])
df_truncated = pd.DataFrame(np.stack((Modes_flag[Zoom_start_index:],(-torch.tensor(error)[Zoom_start_index:]).tolist(), epochs[Zoom_start_index:]) ,axis=1), columns=["Modes_truncated","Loss_truncated","epochs_truncated"])
df_combined = pd.concat([df_full, df_truncated], axis=1)
df_combined = df_combined.astype('float64')
df_combined.to_csv('Results/'+Name+'.csv')

Save trained model#

Save the state dictionary of the trained model

torch.save(ROM_model.state_dict(), 'Pretrained_models/2D_ROM_SV')

2D case - Figure 12: Space variability training

Contents