Eparish1/opinf add neural networks 2

.github/workflows/ci.yaml

@@ -33,14 +33,29 @@ jobs:
         with:
           version: ${{ matrix.version }}
 
+      - name: Set up Python
+        uses: actions/setup-python@v2
+        with:
+          python-version: '3.10'  # Specify your Python version
+
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip 
+          pip install torch
+          pip install numpy
+
       - name: Cache Julia artifacts
         uses: julia-actions/cache@v2
 
       - name: Build project
         uses: julia-actions/julia-buildpkg@v1
+        env:
+          PYTHON: python
 
       - name: Run tests
         uses: julia-actions/julia-runtest@v1
+        with:
+          test_args: "run-nn-opinf-tests"
 
       - name: Process coverage
         uses: julia-actions/julia-processcoverage@v1

Project.toml

@@ -20,6 +20,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
+PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 
 [compat]
 DelimitedFiles = "1"

examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-1-for-test.yaml

@@ -0,0 +1,40 @@
+type: single
+input mesh file: cuboid-1.g
+output mesh file: cuboid-1.e
+model:
+  type: solid mechanics
+  material:
+    blocks:
+      fine: hyperelastic
+    hyperelastic:
+      model: neohookean
+      elastic modulus: 1.0e+09
+      Poisson's ratio: 0.25
+      density: 1000.0
+time integrator:
+  type: Newmark
+  β: 0.25
+  γ: 0.5
+boundary conditions:
+  Dirichlet:
+    - node set: nsx-
+      component: x
+      function: "0.0"
+    - node set: nsy-
+      component: y
+      function: "0.0"
+    - node set: nsz-
+      component: z
+      function: "0.0"
+  Schwarz overlap:
+    - side set: ssz+
+      source: cuboid-2-for-test.yaml
+      source block: coarse
+      source side set: ssz-
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  minimum iterations: 1
+  maximum iterations: 16
+  relative tolerance: 1.0e-10
+  absolute tolerance: 1.0e-06

examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-1.yaml

@@ -0,0 +1,40 @@
+type: single
+input mesh file: cuboid-1.g
+output mesh file: cuboid-1.e
+model:
+  type: solid mechanics
+  material:
+    blocks:
+      fine: hyperelastic
+    hyperelastic:
+      model: neohookean
+      elastic modulus: 1.0e+09
+      Poisson's ratio: 0.25
+      density: 1000.0
+time integrator:
+  type: Newmark
+  β: 0.25
+  γ: 0.5
+boundary conditions:
+  Dirichlet:
+    - node set: nsx-
+      component: x
+      function: "0.0"
+    - node set: nsy-
+      component: y
+      function: "0.0"
+    - node set: nsz-
+      component: z
+      function: "0.0"
+  Schwarz overlap:
+    - side set: ssz+
+      source: cuboid-2.yaml
+      source block: coarse
+      source side set: ssz-
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  minimum iterations: 1
+  maximum iterations: 16
+  relative tolerance: 1.0e-10
+  absolute tolerance: 1.0e-06

examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-2-for-test.yaml

@@ -0,0 +1,41 @@
+type: single
+input mesh file: cuboid-2.g
+output mesh file: cuboid-2.e
+model:
+  type: neural network opinf rom
+  model-directory: neural-network-model/
+  ensemble-size: 0
+  material:
+    blocks:
+      coarse: hyperelastic
+    hyperelastic:
+      model: neohookean
+      elastic modulus: 1.0e+09
+      Poisson's ratio: 0.25
+      density: 1000.0
+time integrator:
+  type: Newmark
+  β: 0.25
+  γ: 0.5
+boundary conditions:
+  OpInf Dirichlet:
+    - node set: nsz+
+      component: z
+      function: "1.0 * t"
+      model-directory: neural-network-model/ 
+      ensemble-size: 0
+
+  OpInf Schwarz overlap:
+    - side set: ssz- 
+      source: cuboid-1-for-test.yaml
+      source block: fine 
+      source side set: ssz+
+      model-directory: neural-network-model/ 
+      ensemble-size: 0
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  minimum iterations: 1
+  maximum iterations: 1
+  relative tolerance: 1.0e-10
+  absolute tolerance: 1.0e-06

examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-2.yaml

@@ -0,0 +1,41 @@
+type: single
+input mesh file: cuboid-2.g
+output mesh file: cuboid-2.e
+model:
+  type: neural network opinf rom
+  model-directory: neural-network-model/
+  ensemble-size: 5
+  material:
+    blocks:
+      coarse: hyperelastic
+    hyperelastic:
+      model: neohookean
+      elastic modulus: 1.0e+09
+      Poisson's ratio: 0.25
+      density: 1000.0
+time integrator:
+  type: Newmark
+  β: 0.25
+  γ: 0.5
+boundary conditions:
+  OpInf Dirichlet:
+    - node set: nsz+
+      component: z
+      function: "1.0 * t"
+      model-directory: neural-network-model/ 
+      ensemble-size: 5
+
+  OpInf Schwarz overlap:
+    - side set: ssz- 
+      source: cuboid-1.yaml
+      source block: fine 
+      source side set: ssz+
+      model-directory: neural-network-model/ 
+      ensemble-size: 5
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  minimum iterations: 1
+  maximum iterations: 16
+  relative tolerance: 1.0e-10
+  absolute tolerance: 1.0e-06

examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboids-for-test.yaml

@@ -0,0 +1,11 @@
+type: multi
+domains: ["cuboid-1-for-test.yaml", "cuboid-2-for-test.yaml"]
+Exodus output interval: 1
+CSV output interval: 1
+initial time: 0.0
+final time: 0.01
+time step: 0.01
+minimum iterations: 1
+maximum iterations: 1
+relative tolerance: 1.0e-12
+absolute tolerance: 1.0e-08

examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboids.yaml

@@ -0,0 +1,11 @@
+type: multi
+domains: ["cuboid-1.yaml", "cuboid-2.yaml"]
+Exodus output interval: 1
+CSV output interval: 1
+initial time: 0.0
+final time: 1.0
+time step: 0.01
+minimum iterations: 1
+maximum iterations: 16
+relative tolerance: 1.0e-12
+absolute tolerance: 1.0e-08

examples/ahead/overlap/torsion/dynamic-opinf-fom/make_op_inf_models.py

@@ -1,117 +1,23 @@
-import numpy as np
-import opinf
-import os
-from matplotlib import pyplot as plt
 import normaopinf
-import normaopinf.readers
-import normaopinf.calculus
-import romtools
-
-if __name__ == "__main__":
-    # Load in snapshots
-    cur_dir = os.getcwd()
-    solution_id = 2
-    displacement_snapshots,times = normaopinf.readers.load_displacement_csv_files(solution_directory=cur_dir,solution_id=solution_id,skip_files=1)
-
-    # Identify which DOFs are free
-    free_dofs = normaopinf.readers.get_free_dofs(solution_directory=cur_dir,solution_id=solution_id)
-
-    # Set values = 0 if DOFs are fixed
-    displacement_snapshots[free_dofs[:,:]==False] = 0.
-
-    #Get sideset snapshots
-    sidesets = ["-Z_ss"]
-    sideset_snapshots = normaopinf.readers.load_sideset_displacement_csv_files(solution_directory=cur_dir,sidesets=sidesets,solution_id=solution_id,skip_files=1)
-
-    # Create an energy-based truncater 
-    #tolerance = 1.e-5
-    #my_energy_truncater = romtools.vector_space.utils.EnergyBasedTruncater(1. - tolerance)
-    my_energy_truncater = romtools.vector_space.utils.BasisSizeTruncater(4)
-
-    # Now load in sidesets and create reduced spaces
-    # Note that I construct a separate basis for each x,y,z component. This isn't necessary 
-    ss_tspace = {}
-    reduced_sideset_snapshots = {}
-    for sideset in sidesets:
-        if sideset_snapshots[sideset].shape[0] ==  1:
-          ss_tspace[sideset] = romtools.VectorSpaceFromPOD(snapshots=sideset_snapshots[sideset],
-                                              truncater=my_energy_truncater,
-                                              shifter = None,
-                                              orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
-                                              scaler = romtools.vector_space.utils.NoOpScaler())
-          # Compute L2 orthogonal projection onto trial spaces
-          reduced_sideset_snapshots[sideset] = romtools.rom.optimal_l2_projection(sideset_snapshots[sideset],ss_tspace[sideset]) 
-        else:
-          comp_trial_space = []
-          for i in range(0,3):
-            tspace = romtools.VectorSpaceFromPOD(snapshots=sideset_snapshots[sideset][i:i+1],
-                                              truncater=my_energy_truncater,
-                                              shifter = None,
-                                              orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
-                                              scaler = romtools.vector_space.utils.NoOpScaler())
-            comp_trial_space.append(tspace)
-          ss_tspace[sideset] = romtools.CompositeVectorSpace(comp_trial_space)
-          reduced_sideset_snapshots[sideset] = romtools.rom.optimal_l2_projection(sideset_snapshots[sideset],ss_tspace[sideset]) 
-
-
-    ## Stack sidesets into one matrix for OpInf 
-    reduced_stacked_sideset_snapshots = None
-    for sideset in sidesets:
-        if reduced_stacked_sideset_snapshots is None:
-            reduced_stacked_sideset_snapshots = reduced_sideset_snapshots[sideset]*1.
-        else: 
-            reduced_stacked_sideset_snapshots = np.append(reduced_stacked_sideset_snapshots,reduced_sideset_snapshots[sideset],axis=0)
-
-    # Create trial space for displacement vector
-    # Note again that I construct a separate basis for each x,y,z component. This isn't necessary 
-    trial_spaces = []
-
-    my_energy_truncater = romtools.vector_space.utils.BasisSizeTruncater(60)
-
-    for i in range(0,3):
-      trial_space = romtools.VectorSpaceFromPOD(snapshots=displacement_snapshots[i:i+1],
-                                              truncater=my_energy_truncater,
-                                              shifter = None,
-                                              orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
-                                              scaler = romtools.vector_space.utils.NoOpScaler())
-      trial_spaces.append(trial_space)
-
-    trial_space = romtools.CompositeVectorSpace(trial_spaces)
-
-    # Compute L2 orthogonal projection onto trial spaces 
-    uhat = romtools.rom.optimal_l2_projection(displacement_snapshots,trial_space)
-    u_ddots,times_dummy = normaopinf.readers.load_acceleration_csv_files(solution_directory=cur_dir,solution_id=solution_id,skip_files=1)
-    # Set values = 0 if DOFs are fixed
-    u_ddots[free_dofs[:,:]==False] = 0.
-    uhat_ddots = romtools.rom.optimal_l2_projection(u_ddots*1.,trial_space)
-
-    # Construct an opinf "AB" model (linear in the state and linear in the exogenous inputs)
-    #   Note: I don't construct a cAB ROM in this example since I know there is no forcing vector
-    l2solver = opinf.lstsq.L2Solver(regularizer=1e-11)
-    opinf_model = opinf.models.ContinuousModel("AB",solver=l2solver)
-    opinf_model.fit(states=uhat, ddts=uhat_ddots,inputs=reduced_stacked_sideset_snapshots)
-
-    # Flip signs to match convention of K on the LHS
-    K = -opinf_model.A_.entries
-    B = opinf_model.B_.entries
-
-    ## Now extract boundary operators and create dictionary to save
-    col_start = 0
-    sideset_operators = {}
-    for sideset in sidesets:
-      num_dofs = reduced_sideset_snapshots[sideset].shape[0]
-      val = np.einsum('kr,vnr->vkn',B[:,col_start:col_start + num_dofs] , ss_tspace[sideset].get_basis() )
-      shape2 = B[:,col_start:col_start + num_dofs] @ ss_tspace[sideset].get_basis()[0].transpose()
-      sideset_operators["B_" + sideset] = val#
-      col_start += num_dofs
-
-    f = np.zeros(K.shape[0])
-    vals_to_save = sideset_operators 
-    vals_to_save["basis"] = trial_space.get_basis() 
-    vals_to_save["K"] = K
-    vals_to_save["f"] = f
-
-    np.savez('opinf-operator',**vals_to_save)
-
-
+import normaopinf.opinf
+import os
+import numpy as np
 
+if __name__ == '__main__':
+  settings = {}
+  settings['fom-yaml-file'] = "torsion.yaml"
+  settings['training-data-directories'] = [os.getcwd()]
+  settings['solution-id'] = 1
+  settings['model-type'] = 'quadratic'
+  settings['forcing'] =  False
+  settings['truncation-type'] = 'size'
+  settings['boundary-truncation-type'] =  'energy' 
+  settings['regularization-parameter'] =  'automatic'
+  settings['trial-space-splitting-type'] = 'split'
+  settings['acceleration-computation-type'] = 'finite-difference'
+  sizes = np.array([2,4,6,8,10,12,14,16,18,20,40,60],dtype=int)
+  for i,size in enumerate(sizes):
+    settings['truncation-value'] = int(size)
+    settings['boundary-truncation-value'] = 1. - 1.e-5 
+    settings['operator-name'] = 'quadratic-opinf-operator-rom-dim-' + str(size)
+    normaopinf.opinf.make_opinf_model(settings)