Eparish1/op inf

Manifest.toml

@@ -1049,3 +1049,23 @@ version = "1.52.0+1"
 deps = ["Artifacts", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
 version = "17.4.0+2"
+
+[[deps.ZipFile]]
+deps = ["Libdl", "Printf", "Zlib_jll"]
+git-tree-sha1 = "f492b7fe1698e623024e873244f10d89c95c340a"
+uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
+version = "0.10.1"
+
+[[deps.FileIO]]
+deps = ["Pkg", "Requires", "UUIDs"]
+git-tree-sha1 = "91e0e5c68d02bcdaae76d3c8ceb4361e8f28d2e9"
+uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+version = "1.16.5"
+
+
+[[deps.NPZ]]
+deps = ["FileIO", "ZipFile"]
+git-tree-sha1 = "60a8e272fe0c5079363b28b0953831e2dd7b7e6f"
+uuid = "15e1cf62-19b3-5cfa-8e77-841668bca605"
+version = "0.4.3"
+

Project.toml

@@ -17,3 +17,4 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
+NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605"

examples/ahead/overlap/cuboid/dynamic-opinf-fom/cuboid-1.yaml

@@ -0,0 +1,41 @@
+type: single
+input mesh file: ../cuboid-1.g
+output mesh file: cuboid-1.e
+CSV write sidesets: true
+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-fom/cuboid-2.yaml

@@ -0,0 +1,41 @@
+type: single
+input mesh file: ../cuboid-2.g
+output mesh file: cuboid-2.e
+CSV write sidesets: true
+model:
+  type: solid mechanics
+  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:
+  Dirichlet:
+    - node set: nsx-
+      component: x
+      function: "0.0"
+    - node set: nsy-
+      component: y
+      function: "0.0"
+    - node set: nsz+
+      component: z
+      function: "1.0 * t"
+  Schwarz overlap:
+    - side set: ssz-
+      source: cuboid-1.yaml
+      source block: fine
+      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-fom/cuboids.yaml

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

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

@@ -0,0 +1,109 @@
+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 = ["nsx-","nsy-","ssz-","nsz+"]
+    sideset_snapshots = normaopinf.readers.load_sideset_displacement_csv_files(solution_directory=cur_dir,sidesets=sidesets,solution_id=solution_id,skip_files=1)
+
+
+    # Create bases for sidesets 
+    my_energy_truncater = romtools.vector_space.utils.EnergyBasedTruncater(1. - 1.e-5)
+
+    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())
+          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]) 
+
+
+    ##accumulate sideset_snapshots into one large data matrix for opinf
+    #stacked_sideset_snapshots = None
+    #for sideset in sidesets:
+    #    if stacked_sideset_snapshots is None:
+    #        stacked_sideset_snapshots = sideset_snapshots[sideset]*1.
+    #    else: 
+    #        stacked_sideset_snapshots = np.append(stacked_sideset_snapshots,sideset_snapshots[sideset],axis=0)
+
+    ## Repeat for reducuced sideset snapshots
+    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)
+
+    trial_spaces = []
+    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)
+
+    uhat = romtools.rom.optimal_l2_projection(displacement_snapshots,trial_space)
+    u_ddots = normaopinf.calculus.d2dx2(displacement_snapshots,times)
+    uhat_ddots = romtools.rom.optimal_l2_projection(u_ddots*1.,trial_space)
+    l2solver = opinf.lstsq.L2Solver(regularizer=1e-8)
+    opinf_model = opinf.models.ContinuousModel("AB",solver=l2solver)
+
+    opinf_model.fit(states=uhat, ddts=uhat_ddots,inputs=reduced_stacked_sideset_snapshots)
+    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)
+
+
+

examples/ahead/overlap/cuboid/dynamic-opinf-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-rom/cuboid-2.yaml

@@ -0,0 +1,41 @@
+type: single
+input mesh file: cuboid-2.g
+output mesh file: cuboid-2.e
+model:
+  type: linear opinf rom
+  model-file: opinf-operator.npz
+  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:
+  Dirichlet:
+    - node set: nsx-
+      component: x
+      function: "0.0"
+    - node set: nsy-
+      component: y
+      function: "0.0"
+    - node set: nsz+
+      component: z
+      function: "1.0 * t"
+  Schwarz overlap:
+    - side set: ssz-
+      source: cuboid-1.yaml
+      source block: fine
+      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-rom/cuboids.yaml

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

examples/ahead/single/cuboid/dynamic-opinf-fom/cuboid.yaml

@@ -0,0 +1,45 @@
+type: single
+input mesh file: ../cuboid.g
+output mesh file: cuboid.e
+Exodus output interval: 1
+CSV output interval: 1
+CSV write sidesets: true
+
+model:
+  type: solid mechanics
+  material:
+    blocks:
+      cuboid: hyperelastic
+    hyperelastic:
+      model: linear elastic 
+      elastic modulus: 1.0e+09
+      Poisson's ratio: 0.25
+      density: 1000.0
+time integrator:
+  type: Newmark
+  β: 0.25
+  γ: 0.5
+  initial time: 0.0
+  final time: 1.0
+  time step: 0.005
+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"
+    - node set: nsz+
+      component: z
+      function: "1.0 * t"
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  minimum iterations: 1
+  maximum iterations: 16
+  relative tolerance: 1.0e-12
+  absolute tolerance: 1.0e-08