Linear elasticity workflow example with meshing of complicated components - LLNL/scaleupROM GitHub Wiki
TL;DR
- A full workflow example for realistic structural analysis of a frame structure with more complicated component geometries.
- We run GMSH to generate the meshes of the components.
- We plot the correlation between relative error, speedup and the number of bases for the ROM.
- The full code is provided at the end of the page.
NB: The reader is assumed to previously have read the linear elasticity workflow example for a simpler FOM
Full order model
The full order model is a two dimensional frame structure subjected to a wind load Pw = 0.5 kN/m and a gravity load Pg. The material of the structure is steel, with a Young's modulus E = 211 GPa. The frame structure has a width w = 0.01m, and a density of 7870 kg/m3, and the gravitational acceleration is assumed to be g = 9.81m/s2, leading to a gravity load Pg = 0.01 * 7870 * 9.81 = 0.772 kN/m2.
Above figure: Full Order Model. Applied wind load in red. Bottom boundary is clamped. Blue text denotes the gravity load. NB: the gravity load is applied at all components
The FOM components in this case consists of a 200mm x 200mm joint component, connecting beam and column components. The beam component has been designed to efficiently carry a uniformly distributed load. It has a length of 6000mm, a smallest heigh of 200mm and a largest height of 400mm. The column is designed for improved buckling resistance. It has a height of 4000mm, a smallest width of 200mm and a largest width of 400mm.
Above figure: FOM components. Dimensions in black text (unit: m), boundary attributes in blue text.
Meshing of components
The complicated component geometry makes it important to properly define the meshes. This can be done quite straightforwardly with the GMSH software. In this section, we provide a walkthrough for how to create the mesh for the beam component using an GMSH script (to be saved in .geo file format). The meshing of the column, and joint components are done in a similar way. The beam geometry is described in the below figure. Note that there are four different boundary attributes that will be assigned to the mesh. Also note that one of the edges is a spline curve.
Above figure: Beam component geometry. Points in local coordinate system in red, boundary attributes in blue text.
The following code defines the unmeshed geometry of the beam:
SetFactory("OpenCASCADE");
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {3, -0.2, 0, 1.0};
//+
Point(3) = {6, 0, 0, 1.0};
//+
Point(4) = {6, 0.2, 0, 1.0};
//+
Point(5) = {0, 0.2, 0, 1.0};
//+
Spline(1) = {1, 2, 3};
//+
Line(2) = {3, 4};
//+
Line(3) = {4, 5};
//+
Line(4) = {5, 1};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
The initial command,SetFactory("OpenCASCADE"), specifies what geometry kernel to use. Another possibility is to use the GMSH built-in kernel with SetFactory("Built-in"). The following commands are quite self explanatory, first points are added and then line / curve segments are drawn between the points. We then create a closed curve from the cruve segments with Curve Loop. Finally, the closed curve is used to define a planar surface with the Plane Surface command.
Physical Curve(1) = {1};
Physical Curve(2) = {2};
Physical Curve(3) = {3};
Physical Curve(4) = {4};
Physical Surface(1) = {1};
The command Physical Curve is used to assign boundary attributes to the curves. More specifically, Physical Curve(i) = {S}, where i is an integer and S is a list of integers corresponding to the curve indicies. In this case curve 1 gets boundary attribute 1 etc... The command Physical Surface works in a similar way for surfaces.
Mesh.Algorithm = 1;
Mesh.AnisoMax = 0;
Mesh.MshFileVersion = 2.2;
Mesh 2; // This is important
RefineMesh;
RefineMesh;
Save "optbeam.msh";
Exit;
The first four lines does not need to be changed to create a different mesh file. The curious reader can consult the GMSH docs for more information. Especially the line Mesh.MshFileVersion = 2.2; needs to be kept as compatibility with MFEM is broken otherwise. The command RefineMesh; performs a uniform mesh refinement step.
The program to generate the beam mesh, optbeam.geo is presented in its entirety below:
SetFactory("OpenCASCADE");
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {3, -0.2, 0, 1.0};
//+
Point(3) = {6, 0, 0, 1.0};
//+
Point(4) = {6, 0.2, 0, 1.0};
//+
Point(5) = {0, 0.2, 0, 1.0};
//+
Spline(1) = {1, 2, 3};
//+
Line(2) = {3, 4};
//+
Line(3) = {4, 5};
//+
Line(4) = {5, 1};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
Physical Curve(1) = {1};
Physical Curve(2) = {2};
Physical Curve(3) = {3};
Physical Curve(4) = {4};
Physical Surface(1) = {1};
Mesh.Algorithm = 1;
Mesh.AnisoMax = 0;
Mesh.MshFileVersion = 2.2;
Mesh 2; // This is important
RefineMesh;
RefineMesh;
Save "optbeam.msh";
Exit;
Finally, we run the following BASH script to create the meshes with GMSH, then convert them to an MFEM compatible format and copy them over to the correct folders.
cd ~/build/utils
gmsh ~/utils/gmsh/optjoint.geo -2
gmsh ~/utils/gmsh/optbeam.geo -2
gmsh ~/utils/gmsh/optcol.geo -2
./gmsh2mfem -m ~/utils/gmsh/optjoint.msh -o 1 -fnc
./gmsh2mfem -m ~/utils/gmsh/optbeam.msh -o 1 -fnc
./gmsh2mfem -m ~/utils/gmsh/optcol.msh -o 1 -fnc
rm ~/examples/linelast/meshes/optjoint.mesh
rm ~/examples/linelast/meshes/optbeam.mesh
rm ~/examples/linelast/meshes/optcol.mesh
rm ~/build/examples/linelast/meshes/optjoint.mesh
rm ~/build/examples/linelast/meshes/optbeam.mesh
rm ~/build/examples/linelast/meshes/optcol.mesh
mv ~/utils/gmsh/optjoint.msh.mfem ~/examples/linelast/meshes/optjoint.mesh
mv ~/utils/gmsh/optbeam.msh.mfem ~/examples/linelast/meshes/optbeam.mesh
mv ~/utils/gmsh/optcol.msh.mfem ~/examples/linelast/meshes/optcol.mesh
cp ~/examples/linelast/meshes/optjoint.mesh ~/build/examples/linelast/meshes
cp ~/examples/linelast/meshes/optbeam.mesh ~/build/examples/linelast/meshes
cp ~/examples/linelast/meshes/optcol.mesh ~/build/examples/linelast/meshes
Snapshot generation
For the snapshot generation sequence, we use a small FOM and the same overall strategy with four different training scenarios. In each scenario, one side of the model has prescribed Dirichlet BCs and the remaining three sides of the model are subjected to parameterized external loadings as Neumann BCs:
$$ F_i(x,y) = \begin{pmatrix} f_i^x + \alpha_i^x * \sin\left(2\pi \cdot \lambda_i^x \cdot \left(\frac{x}{6.0} + \gamma_i^x \right)\right) \text{if} :p_i^x \geq 0.5 \\ f_i^y + \alpha_i^y * \sin\left(2\pi \cdot \lambda_i^y \cdot \left(\frac{y}{6.0} + \gamma_i^y \right)\right) \text{if} :p_i^y \geq 0.5 \end{pmatrix}, i \in {1, 2, 3, 4}$$
The FOM of this example differs from the previous one in that there is a body force corresponding to the gravitational load present. This motivates us to include a parameterized body force to the training sets, this force is evenly distributed over the small FOM training model for each snapshot generation case and denoted F4 in the below picture.
Above figure: Parameterized Problem Permutations. Prescribed displacements in red. Applied external forces in blue. Body force in purple. BC setup ID in black.
The following ranges where used to sample our parameters:
$$
\begin{align*}
-50 \cdot 10^3 &\leq f_i^x, f_i^y \leq 50 \cdot 10^3, i \in {1, 2, 3} \
-5 \cdot 10^3 &\leq f_i^x, f_i^y \leq 5 \cdot 10^3 , i \in {4} \
-50 \cdot 10^3 &\leq \alpha_i^x, \alpha_i^y \leq 50 \cdot 10^3, i \in {1, 2, 3} \
-5 \cdot 10^3 &\leq \alpha_i^x, \alpha_i^y \leq 5 \cdot 10^3, i \in {4} \
0.1 &\leq \lambda_i^x, \lambda_i^y \leq 1, i \in {1, \ldots, 4} \
0 &\leq \gamma_i^x, \gamma_i^y \leq 1, i \in {1, \ldots, 4} \
0 &\leq p_i^x, p_i^y \leq 1, i \in {1, \ldots, 4} \
\end{align*}
$$
# Subcomponent sampling
cd ~/scaleupROM
# Create configs for componentwise training
python3 utils/python/linelast_cwtrain.py cwtrain_optROM
# Copy to examples and build folders
for i in {0..3}
do
rm examples/linelast/config/linelast.comp_train$i.h5
cp linelast.comp_train$i.h5 examples/linelast/config
rm build/examples/linelast/config/linelast.comp_train$i.h5
cp linelast.comp_train$i.h5 build/examples/linelast/config
rm linelast.comp_train$i.h5
done
cd build/examples/linelast
# Sample generation
n_samp=20000
for i in {0..3}
do
../../bin/main -i linelast.opt_comp_train.yml -f main/mode=sample_generation:mesh/component-wise/global_config="config/linelast.comp_train$i.h5":sample_generation/file_path/prefix="linelast_cwtrain$i":sample_generation/random_sample_generator/number_of_samples=$n_samp:single_run/linelast_cwtrain/lambda=$l:single_run/linelast_cwtrain/mu=$m
done
# Rename samples to fit format
for i in {0..3}
do
echo "Rename $i"
mv linelast_cwtrain$i\_sample_optjoint_snapshot.000000 linelast_cwtrain_sample_optjoint_snapshot$i\.000000
mv linelast_cwtrain$i\_sample_optbeam_snapshot.000000 linelast_cwtrain_sample_optbeam_snapshot$i\.000000
mv linelast_cwtrain$i\_sample_optcol_snapshot.000000 linelast_cwtrain_sample_optcol_snapshot$i\.000000
done
Train ROM
The ROM is trained in virtually the same manner as the previous linear elasticity example:
../../bin/main -i linelast.opt_comp_train.yml -f main/mode=train_rom:basis/tags/0/snapshot_format/maximum=3:basis/tags/1/snapshot_format/maximum=3:basis/tags/2/snapshot_format/maximum=3
The spectrum of the ROM basis vectors can be visualized with the following Python code:
import os
os.chdir('/Users/larsson4/repos/scaleupROM/build/examples/linelast')
def get_svs(filename):
with open(filename, 'r') as file:
return [float(line.strip()) for line in file]
svs_j = get_svs('linelast_cwtrain_optjoint_sv.txt')
svs_b = get_svs('linelast_cwtrain_optbeam_sv.txt')
svs_c = get_svs('linelast_cwtrain_optcol_sv.txt')
import matplotlib.pyplot as plt
plt.plot(range(len(svs_j)), svs_j, label='Joint')
plt.plot(range(len(svs_b)), svs_b, label='Beam')
plt.plot(range(len(svs_c)), svs_c, label='Column')
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.title('Advanced component SV spectrum')
Above figure: Singular value spectrum for advanced component ROMs
Large scale prediction
We run the basis scaling in a similar way to the previous example for linear elasticity. A notable change here is that we also scale the relative size of the bases for the beam and column component ROMs proportional to their relative number of DOFs, compared to the joint component. I.e. our n is the number of bases for the joint component, with the largest possible number being 384, which is the number of DOFs in the FOM joint component. The FOM beam component has approximately 3x as many DOFs as the joint component, so nb, the number of bases for the beam component, gets scaled by three compared to n. Similarily, since the FOM column component has approximately 2x as many DOFs as the joint component, the number of bases for the column component, nc, gets scaled by 2x compared to n.
# Run basis scaling
mkdir basis_scaling
for n in 1 2 3 4 6 7 9 12 15 20 25 32 41 53 68 87 111 142 183 234 300 384
do
nb=$(($n*3))
nc=$(($n*2))
echo "Prediction $n"
# Building ROM
../../bin/main -i linelast.opt_comp_train.yml -f main/mode=build_rom:main/use_rom=true:basis/tags/0/number_of_basis=$n:basis/tags/1/number_of_basis=$nb:basis/tags/2/number_of_basis=$nc:mesh/component-wise/global_config="config/linelast.optfom.h5":single_run/linelast_cwtrain/lambda=$l:single_run/linelast_cwtrain/mu=$m
# Single run with larger FOM configuration
../../bin/main -i linelast.opt_comp_train.yml -f main/use_rom=true:mesh/component-wise/global_config="config/linelast.optfom.h5":parameterized_problem/name="linelast_frame_wind":single_run/linelast_frame_wind/qwind_f=500.0:single_run/linelast_frame_wind/g=9.81:basis/tags/0/number_of_basis=$n:basis/tags/1/number_of_basis=$nb:basis/tags/2/number_of_basis=$nc -o "basis_scaling/comparison$n.h5"
done
Plot results
We run the same plotting command as in the previous linear elasticity example:
# Plot results
python3 ../../../utils/python/linelast_cwtrain.py bscale basis_scaling scaling.png
Full program
cd ~/build/utils
gmsh ~/utils/gmsh/optjoint.geo -2
gmsh ~/utils/gmsh/optbeam.geo -2
gmsh ~/utils/gmsh/optcol.geo -2
./gmsh2mfem -m ~/utils/gmsh/optjoint.msh -o 1 -fnc
./gmsh2mfem -m ~/utils/gmsh/optbeam.msh -o 1 -fnc
./gmsh2mfem -m ~/utils/gmsh/optcol.msh -o 1 -fnc
rm ~/examples/linelast/meshes/optjoint.mesh
rm ~/examples/linelast/meshes/optbeam.mesh
rm ~/examples/linelast/meshes/optcol.mesh
rm ~/build/examples/linelast/meshes/optjoint.mesh
rm ~/build/examples/linelast/meshes/optbeam.mesh
rm ~/build/examples/linelast/meshes/optcol.mesh
mv ~/utils/gmsh/optjoint.msh.mfem ~/examples/linelast/meshes/optjoint.mesh
mv ~/utils/gmsh/optbeam.msh.mfem ~/examples/linelast/meshes/optbeam.mesh
mv ~/utils/gmsh/optcol.msh.mfem ~/examples/linelast/meshes/optcol.mesh
cp ~/examples/linelast/meshes/optjoint.mesh ~/build/examples/linelast/meshes
cp ~/examples/linelast/meshes/optbeam.mesh ~/build/examples/linelast/meshes
cp ~/examples/linelast/meshes/optcol.mesh ~/build/examples/linelast/meshes
# Subcomponent sampling
cd ~
# Create configs for componentwise training
python3 utils/python/linelast_cwtrain.py cwtrain_optROM
# Copy to examples and build folders
for i in {0..3}
do
rm examples/linelast/config/linelast.comp_train$i.h5
cp linelast.comp_train$i.h5 examples/linelast/config
rm build/examples/linelast/config/linelast.comp_train$i.h5
cp linelast.comp_train$i.h5 build/examples/linelast/config
rm linelast.comp_train$i.h5
done
cd build/examples/linelast
# Sample generation
n_samp=20000
for i in {0..3}
do
../../bin/main -i linelast.opt_comp_train.yml -f main/mode=sample_generation:mesh/component-wise/global_config="config/linelast.comp_train$i.h5":sample_generation/file_path/prefix="linelast_cwtrain$i":sample_generation/random_sample_generator/number_of_samples=$n_samp:single_run/linelast_cwtrain/lambda=$l:single_run/linelast_cwtrain/mu=$m
done
# Rename samples to fit format
for i in {0..3}
do
echo "Rename $i"
mv linelast_cwtrain$i\_sample_optjoint_snapshot.000000 linelast_cwtrain_sample_optjoint_snapshot$i\.000000
mv linelast_cwtrain$i\_sample_optbeam_snapshot.000000 linelast_cwtrain_sample_optbeam_snapshot$i\.000000
mv linelast_cwtrain$i\_sample_optcol_snapshot.000000 linelast_cwtrain_sample_optcol_snapshot$i\.000000
done
../../bin/main -i linelast.opt_comp_train.yml -f main/mode=train_rom:basis/tags/0/snapshot_format/maximum=3:basis/tags/1/snapshot_format/maximum=3:basis/tags/2/snapshot_format/maximum=3
# Run basis scaling
mkdir basis_scaling
for n in 1 2 3 4 6 7 9 12 15 20 25 32 41 53 68 87 111 142 183 234 300 384
do
nb=$(($n*3))
nc=$(($n*2))
echo "Prediction $n"
# Building ROM
../../bin/main -i linelast.opt_comp_train.yml -f main/mode=build_rom:main/use_rom=true:basis/tags/0/number_of_basis=$n:basis/tags/1/number_of_basis=$nb:basis/tags/2/number_of_basis=$nc:mesh/component-wise/global_config="config/linelast.optfom.h5":single_run/linelast_cwtrain/lambda=$l:single_run/linelast_cwtrain/mu=$m
# Single run with larger FOM configuration
../../bin/main -i linelast.opt_comp_train.yml -f main/use_rom=true:mesh/component-wise/global_config="config/linelast.optfom.h5":parameterized_problem/name="linelast_frame_wind":single_run/linelast_frame_wind/qwind_f=500.0:single_run/linelast_frame_wind/g=9.81:basis/tags/0/number_of_basis=$n:basis/tags/1/number_of_basis=$nb:basis/tags/2/number_of_basis=$nc -o "basis_scaling/comparison$n.h5"
done
# Plot results
python3 ../../../utils/python/linelast_cwtrain.py bscale basis_scaling scaling.png