diff --git a/single-node-refactor/contact_test_inputs/block_sliding_test.yaml b/single-node-refactor/contact_test_inputs/block_sliding_test.yaml new file mode 100644 index 000000000..50853bdfc --- /dev/null +++ b/single-node-refactor/contact_test_inputs/block_sliding_test.yaml @@ -0,0 +1,100 @@ +# num_dims: 3 + +dynamic_options: + time_final: 3.0 + dt_min: 1e-6 + dt_max: 0.05 + dt_start: 1.e-5 + cycle_stop: 300000 + + +mesh_options: + source: file + file_path: contact_test_inputs/contact_test.geo + +#mesh_options: +# source: generate +# num_dims: 3 +# type: Box +# origin: [0.0, 0.0, 0.0] +# length: [1.2, 1.2, 1.2] +# num_elems: [12, 12, 12] + + +output_options: + timer_output_level: thorough + output_file_format: ensight + graphics_time_step: 0.02 + # graphics_iteration_step: 10 + +solver_options: + - solver: + method: SGH + # solver_vars: + # - blah + # - blah + # - blah + +boundary_conditions: + + # Tag Z plane + - boundary_condition: + solver: SGH + geometry: z_plane + direction: z_dir + value: 0.0 + type: reflected + + # global contact condition (this should be used when it is not clear what the contact surface(s) will be) + - boundary_condition: + solver: SGH + geometry: global + type: contact + +materials: + - material: + id: 0 + eos_model: ideal_gas + # strength_model: none + elastic_modulus: 10 + poisson_ratio: 0.3 + q1: 1.0 + q2: 1.333 + q1ex: 1.0 + q2ex: 1.333 + eos_global_vars: + - 10 + - 1.0E-14 + - 1.0 + +regions: + - fill_volume: + type: global + material_id: 0 + den: 1.0 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: 0.0 + w: 0.0 + + + + # top block initial velocity conditions + - fill_volume: + type: box + origin: [0.0, 0.0, 0.0] + x1: -0.25 + x2: 0.25 + y1: -0.25 + y2: 0.25 + z1: 1 + z2: 1.5 + material_id: 0 + den: 10.0 + # ie: 0.25833839995946534 + sie: 1.e-5 + velocity: cartesian + u: 0.0 + v: 0.2 + w: -0.001 \ No newline at end of file diff --git a/single-node-refactor/contact_test_inputs/contact_sie_test.yaml b/single-node-refactor/contact_test_inputs/contact_sie_test.yaml new file mode 100644 index 000000000..a005304e8 --- /dev/null +++ b/single-node-refactor/contact_test_inputs/contact_sie_test.yaml @@ -0,0 +1,100 @@ +# num_dims: 3 + +dynamic_options: + time_final: 3.0 + dt_min: 0.02 + dt_max: 0.05 + dt_start: 1.e-2 + cycle_stop: 300000 + + +mesh_options: + source: file + file_path: contact_test_inputs/contact_test.geo + +#mesh_options: +# source: generate +# num_dims: 3 +# type: Box +# origin: [0.0, 0.0, 0.0] +# length: [1.2, 1.2, 1.2] +# num_elems: [12, 12, 12] + + +output_options: + timer_output_level: thorough + output_file_format: ensight + graphics_time_step: 0.02 + # graphics_iteration_step: 10 + +solver_options: + - solver: + method: SGH + # solver_vars: + # - blah + # - blah + # - blah + +boundary_conditions: + + # Tag Z plane + - boundary_condition: + solver: SGH + geometry: z_plane + direction: z_dir + value: 0.0 + type: reflected + + # global contact condition (this should be used when it is not clear what the contact surface(s) will be) + - boundary_condition: + solver: SGH + geometry: global + type: contact + +materials: + - material: + id: 0 + eos_model: ideal_gas + # strength_model: none + elastic_modulus: 10 + poisson_ratio: 0.3 + q1: 1.0 + q2: 1.333 + q1ex: 1.0 + q2ex: 1.333 + eos_global_vars: + - 3 + - 1.0E-14 + - 1.0 + +regions: + - fill_volume: + type: global + material_id: 0 + den: 10.0 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: 0.0 + w: 0.0 + + + + # top block initial velocity conditions + - fill_volume: + type: box + origin: [0.0, 0.0, 0.0] + x1: -1.0 + x2: 1.0 + y1: -1.0 + y2: 1.0 + z1: 0.0 + z2: 1.5 + material_id: 0 + den: 1.0 + # ie: 0.25833839995946534 + sie: 0.5 + velocity: cartesian + u: 0.0 + v: 0.0 + w: 0.0 \ No newline at end of file diff --git a/single-node-refactor/contact_test_inputs/contact_test.geo b/single-node-refactor/contact_test_inputs/contact_test.geo new file mode 100644 index 000000000..b647d563a --- /dev/null +++ b/single-node-refactor/contact_test_inputs/contact_test.geo @@ -0,0 +1,94 @@ +Description line 1 +Description line 2 +node id assign +element id assign +part + 1 +Mesh +coordinates + 26 +1.00000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +-1.00000e+00 +-1.00000e+00 +-1.00000e+00 +-1.00000e+00 +-1.00000e+00 +-1.00000e+00 +-2.50000e-01 +-2.50000e-01 +-2.50000e-01 +-2.50000e-01 +2.50000e-01 +2.50000e-01 +2.50000e-01 +2.50000e-01 +-1.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +-1.00000e+00 +-1.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +-1.00000e+00 +-1.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +-1.00000e+00 +-2.50000e-01 +2.50000e-01 +2.50000e-01 +-2.50000e-01 +-2.50000e-01 +2.50000e-01 +2.50000e-01 +-2.50000e-01 +0.00000e+00 +0.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +1.50000e+00 +1.50000e+00 +1.00000e+00 +1.00000e+00 +1.50000e+00 +1.50000e+00 +hexa8 + 5 + 7 8 11 12 1 2 5 6 + 8 9 10 11 2 3 4 5 + 13 14 17 18 7 8 11 12 + 14 15 16 17 8 9 10 11 + 19 20 21 22 23 24 25 26 diff --git a/single-node-refactor/contact_test_inputs/contact_test.yaml b/single-node-refactor/contact_test_inputs/contact_test.yaml new file mode 100644 index 000000000..320c4ba6b --- /dev/null +++ b/single-node-refactor/contact_test_inputs/contact_test.yaml @@ -0,0 +1,100 @@ +# num_dims: 3 + +dynamic_options: + time_final: 1.0 + dt_min: 0.5 + dt_max: 0.5 + dt_start: 1.e-2 + cycle_stop: 300000 + + +mesh_options: + source: file + file_path: contact_test_inputs/contact_test.geo + +#mesh_options: +# source: generate +# num_dims: 3 +# type: Box +# origin: [0.0, 0.0, 0.0] +# length: [1.2, 1.2, 1.2] +# num_elems: [12, 12, 12] + + +output_options: + timer_output_level: thorough + output_file_format: ensight + graphics_time_step: 0.5 + # graphics_iteration_step: 10 + +solver_options: + - solver: + method: SGH + # solver_vars: + # - blah + # - blah + # - blah + +boundary_conditions: + + # Tag Z plane + - boundary_condition: + solver: SGH + geometry: z_plane + direction: z_dir + value: 0.0 + type: reflected + + # global contact condition (this should be used when it is not clear what the contact surface(s) will be) + - boundary_condition: + solver: SGH + geometry: global + type: contact + +materials: + - material: + id: 0 + eos_model: ideal_gas + # strength_model: none + elastic_modulus: 10 + poisson_ratio: 0.3 + q1: 1.0 + q2: 1.333 + q1ex: 1.0 + q2ex: 1.333 + eos_global_vars: + - 3 + - 1.0E-14 + - 1.0 + +regions: + - fill_volume: + type: global + material_id: 0 + den: 1.0 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: 0.0 + w: 0.0 + + + + # top block initial velocity conditions + - fill_volume: + type: box + origin: [0.0, 0.0, 0.0] + x1: -0.25 + x2: 0.25 + y1: -0.25 + y2: 0.25 + z1: 1 + z2: 1.5 + material_id: 0 + den: 1.0 + # ie: 0.25833839995946534 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: 0.0 + w: -0.5 \ No newline at end of file diff --git a/single-node-refactor/contact_test_inputs/contact_test2.geo b/single-node-refactor/contact_test_inputs/contact_test2.geo new file mode 100644 index 000000000..bfeb2d19b --- /dev/null +++ b/single-node-refactor/contact_test_inputs/contact_test2.geo @@ -0,0 +1,189 @@ +Description line 1 +Description line 2 +node id assign +element id assign +part + 1 +Mesh +coordinates + 54 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +5.00000e-01 +1.00000e+00 +5.00000e-01 +0.00000e+00 +0.00000e+00 +1.00000e+00 +5.00000e-01 +1.00000e+00 +1.00000e+00 +0.00000e+00 +5.00000e-01 +0.00000e+00 +5.00000e-01 +5.00000e-01 +1.00000e+00 +5.00000e-01 +0.00000e+00 +5.00000e-01 +5.00000e-01 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +5.00000e-01 +1.00000e+00 +5.00000e-01 +0.00000e+00 +0.00000e+00 +1.00000e+00 +5.00000e-01 +1.00000e+00 +1.00000e+00 +0.00000e+00 +5.00000e-01 +0.00000e+00 +5.00000e-01 +5.00000e-01 +1.00000e+00 +5.00000e-01 +0.00000e+00 +5.00000e-01 +5.00000e-01 +0.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +5.00000e-01 +1.00000e+00 +5.00000e-01 +0.00000e+00 +0.00000e+00 +0.00000e+00 +1.00000e+00 +5.00000e-01 +1.00000e+00 +1.00000e+00 +5.00000e-01 +5.00000e-01 +0.00000e+00 +5.00000e-01 +1.00000e+00 +5.00000e-01 +5.00000e-01 +5.00000e-01 +1.00000e+00 +1.00000e+00 +2.00000e+00 +2.00000e+00 +1.00000e+00 +1.00000e+00 +2.00000e+00 +2.00000e+00 +1.00000e+00 +1.50000e+00 +2.00000e+00 +1.50000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +2.00000e+00 +1.50000e+00 +2.00000e+00 +2.00000e+00 +1.50000e+00 +1.50000e+00 +1.00000e+00 +1.50000e+00 +2.00000e+00 +1.50000e+00 +1.50000e+00 +1.50000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +5.00000e-01 +5.00000e-01 +1.00000e+00 +5.00000e-01 +1.00000e+00 +5.00000e-01 +1.00000e+00 +1.00000e+00 +0.00000e+00 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +1.00000e+00 +5.00000e-01 +0.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +5.00000e-01 +5.00000e-01 +1.00000e+00 +5.00000e-01 +1.00000e+00 +5.00000e-01 +1.00000e+00 +1.00000e+00 +0.00000e+00 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +1.00000e+00 +5.00000e-01 +hexa8 + 16 + 1 9 21 12 13 22 27 25 + 13 22 27 25 5 15 26 20 + 9 2 10 21 22 14 23 27 + 22 14 23 27 15 6 17 26 + 12 21 11 4 25 27 24 18 + 25 27 24 18 20 26 19 8 + 21 10 3 11 27 23 16 24 + 27 23 16 24 26 17 7 19 + 28 36 48 39 40 49 54 52 + 40 49 54 52 32 42 53 47 + 36 29 37 48 49 41 50 54 + 49 41 50 54 42 33 44 53 + 39 48 38 31 52 54 51 45 + 52 54 51 45 47 53 46 35 + 48 37 30 38 54 50 43 51 + 54 50 43 51 53 44 34 46 diff --git a/single-node-refactor/contact_test_inputs/contact_test2.yaml b/single-node-refactor/contact_test_inputs/contact_test2.yaml new file mode 100644 index 000000000..671c70b43 --- /dev/null +++ b/single-node-refactor/contact_test_inputs/contact_test2.yaml @@ -0,0 +1,100 @@ +# num_dims: 3 + +dynamic_options: + time_final: 1.0 + dt_min: 0.5 + dt_max: 0.5 + dt_start: 1.e-2 + cycle_stop: 300000 + + +mesh_options: + source: file + file_path: contact_test_inputs/contact_test2.geo + +#mesh_options: +# source: generate +# num_dims: 3 +# type: Box +# origin: [0.0, 0.0, 0.0] +# length: [1.2, 1.2, 1.2] +# num_elems: [12, 12, 12] + + +output_options: + timer_output_level: thorough + output_file_format: ensight + graphics_time_step: 0.5 + # graphics_iteration_step: 10 + +solver_options: + - solver: + method: SGH + # solver_vars: + # - blah + # - blah + # - blah + +boundary_conditions: + + # Tag Z plane + - boundary_condition: + solver: SGH + geometry: z_plane + direction: z_dir + value: 0.0 + type: reflected + + # global contact condition (this should be used when it is not clear what the contact surface(s) will be) + - boundary_condition: + solver: SGH + geometry: global + type: contact + +materials: + - material: + id: 0 + eos_model: ideal_gas + # strength_model: none + elastic_modulus: 10 + poisson_ratio: 0.3 + q1: 1.0 + q2: 1.333 + q1ex: 1.0 + q2ex: 1.333 + eos_global_vars: + - 3 + - 1.0E-14 + - 1.0 + +regions: + - fill_volume: + type: global + material_id: 0 + den: 1.0 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: 0.0 + w: 0.0 + + + + # top block initial velocity conditions + - fill_volume: + type: box + origin: [0.0, 0.0, 0.0] + x1: 0 + x2: 1 + y1: 1 + y2: 2 + z1: 0 + z2: 1 + material_id: 0 + den: 1.0 + # ie: 0.25833839995946534 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: -0.5 + w: 0.0 \ No newline at end of file diff --git a/single-node-refactor/contact_test_inputs/edge_case1.geo b/single-node-refactor/contact_test_inputs/edge_case1.geo new file mode 100644 index 000000000..55b349045 --- /dev/null +++ b/single-node-refactor/contact_test_inputs/edge_case1.geo @@ -0,0 +1,74 @@ +Description line 1 +Description line 2 +node id assign +element id assign +part + 1 +Mesh +coordinates + 20 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +0.00000e+00 +0.00000e+00 +0.00000e+00 +0.00000e+00 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +-2.00000e+00 +0.00000e+00 +2.00000e+00 +-2.00000e+00 +0.00000e+00 +2.00000e+00 +-2.00000e+00 +0.00000e+00 +2.00000e+00 +-2.00000e+00 +0.00000e+00 +2.00000e+00 +1.00000e+00 +2.00000e+00 +2.00000e+00 +1.00000e+00 +1.00000e+00 +2.00000e+00 +2.00000e+00 +1.00000e+00 +1.00000e+00 +2.00000e+00 +1.00000e+00 +0.00000e+00 +1.00000e+00 +0.00000e+00 +1.00000e+00 +2.00000e+00 +1.00000e+00 +0.00000e+00 +1.00000e+00 +0.00000e+00 +1.85000e+00 +1.85000e+00 +2.85000e+00 +2.85000e+00 +1.85000e+00 +1.85000e+00 +2.85000e+00 +2.85000e+00 +hexa8 + 3 + 7 8 11 10 1 2 5 4 + 8 9 12 11 2 3 6 5 + 13 14 15 16 17 18 19 20 diff --git a/single-node-refactor/contact_test_inputs/edge_case1.yaml b/single-node-refactor/contact_test_inputs/edge_case1.yaml new file mode 100644 index 000000000..d573a1116 --- /dev/null +++ b/single-node-refactor/contact_test_inputs/edge_case1.yaml @@ -0,0 +1,100 @@ +# num_dims: 3 + +dynamic_options: + time_final: 1.0 + dt_min: 0.5 + dt_max: 0.5 + dt_start: 1.e-2 + cycle_stop: 300000 + + +mesh_options: + source: file + file_path: contact_test_inputs/edge_case1.geo + +#mesh_options: +# source: generate +# num_dims: 3 +# type: Box +# origin: [0.0, 0.0, 0.0] +# length: [1.2, 1.2, 1.2] +# num_elems: [12, 12, 12] + + +output_options: + timer_output_level: thorough + output_file_format: ensight + graphics_time_step: 0.5 + # graphics_iteration_step: 10 + +solver_options: + - solver: + method: SGH + # solver_vars: + # - blah + # - blah + # - blah + +boundary_conditions: + + # Tag Z plane + - boundary_condition: + solver: SGH + geometry: z_plane + direction: z_dir + value: 0.0 + type: reflected + + # global contact condition (this should be used when it is not clear what the contact surface(s) will be) + - boundary_condition: + solver: SGH + geometry: global + type: contact + +materials: + - material: + id: 0 + eos_model: ideal_gas + # strength_model: none + elastic_modulus: 10 + poisson_ratio: 0.3 + q1: 1.0 + q2: 1.333 + q1ex: 1.0 + q2ex: 1.333 + eos_global_vars: + - 3 + - 1.0E-14 + - 1.0 + +regions: + - fill_volume: + type: global + material_id: 0 + den: 1.0 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: 0.0 + w: 0.0 + + + + # top block initial velocity conditions + - fill_volume: + type: box + origin: [0.0, 0.0, 0.0] + x1: 0.0 + x2: 0.5 + y1: 1.0 + y2: 2.0 + z1: 1.8 + z2: 2.9 + material_id: 0 + den: 1.0 + # ie: 0.25833839995946534 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: -1.5 + w: 0.0 \ No newline at end of file diff --git a/single-node-refactor/contact_test_inputs/edge_case2.geo b/single-node-refactor/contact_test_inputs/edge_case2.geo new file mode 100644 index 000000000..693a6cb5e --- /dev/null +++ b/single-node-refactor/contact_test_inputs/edge_case2.geo @@ -0,0 +1,74 @@ +Description line 1 +Description line 2 +node id assign +element id assign +part + 1 +Mesh +coordinates + 20 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +2.50000e-01 +2.50000e-01 +2.50000e-01 +2.50000e-01 +-2.00000e+00 +0.00000e+00 +2.00000e+00 +-2.00000e+00 +0.00000e+00 +2.00000e+00 +-2.00000e+00 +0.00000e+00 +2.00000e+00 +-2.00000e+00 +0.00000e+00 +2.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +1.00000e+00 +2.00000e+00 +1.00000e+00 +0.00000e+00 +1.00000e+00 +0.00000e+00 +1.00000e+00 +2.00000e+00 +1.00000e+00 +0.00000e+00 +1.00000e+00 +0.00000e+00 +2.50000e+00 +2.50000e+00 +3.50000e+00 +3.50000e+00 +2.50000e+00 +2.50000e+00 +3.50000e+00 +3.50000e+00 +hexa8 + 3 + 7 8 11 10 1 2 5 4 + 8 9 12 11 2 3 6 5 + 13 14 15 16 17 18 19 20 diff --git a/single-node-refactor/contact_test_inputs/edge_case2.yaml b/single-node-refactor/contact_test_inputs/edge_case2.yaml new file mode 100644 index 000000000..b9f08a10c --- /dev/null +++ b/single-node-refactor/contact_test_inputs/edge_case2.yaml @@ -0,0 +1,100 @@ +# num_dims: 3 + +dynamic_options: + time_final: 1.0 + dt_min: 0.5 + dt_max: 0.5 + dt_start: 1.e-2 + cycle_stop: 300000 + + +mesh_options: + source: file + file_path: contact_test_inputs/edge_case2.geo + +#mesh_options: +# source: generate +# num_dims: 3 +# type: Box +# origin: [0.0, 0.0, 0.0] +# length: [1.2, 1.2, 1.2] +# num_elems: [12, 12, 12] + + +output_options: + timer_output_level: thorough + output_file_format: ensight + graphics_time_step: 0.5 + # graphics_iteration_step: 10 + +solver_options: + - solver: + method: SGH + # solver_vars: + # - blah + # - blah + # - blah + +boundary_conditions: + + # Tag Z plane + - boundary_condition: + solver: SGH + geometry: z_plane + direction: z_dir + value: 0.0 + type: reflected + + # global contact condition (this should be used when it is not clear what the contact surface(s) will be) + - boundary_condition: + solver: SGH + geometry: global + type: contact + +materials: + - material: + id: 0 + eos_model: ideal_gas + # strength_model: none + elastic_modulus: 10 + poisson_ratio: 0.3 + q1: 1.0 + q2: 1.333 + q1ex: 1.0 + q2ex: 1.333 + eos_global_vars: + - 3 + - 1.0E-14 + - 1.0 + +regions: + - fill_volume: + type: global + material_id: 0 + den: 1.0 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: 0.0 + w: 0.0 + + + + # top block initial velocity conditions + - fill_volume: + type: box + origin: [0.0, 0.0, 0.0] + x1: -0.5 + x2: 0.0 + y1: 0.0 + y2: 1.0 + z1: 2.5 + z2: 3.5 + material_id: 0 + den: 1.0 + # ie: 0.25833839995946534 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: 0.0 + w: -1.0 \ No newline at end of file diff --git a/single-node-refactor/contact_test_inputs/edge_case3.geo b/single-node-refactor/contact_test_inputs/edge_case3.geo new file mode 100644 index 000000000..2ad6094d5 --- /dev/null +++ b/single-node-refactor/contact_test_inputs/edge_case3.geo @@ -0,0 +1,87 @@ +Description line 1 +Description line 2 +node id assign +element id assign +part + 1 +Mesh +coordinates + 24 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +-5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +5.00000e-01 +-2.00000e+00 +0.00000e+00 +2.00000e+00 +-2.00000e+00 +0.00000e+00 +2.00000e+00 +-2.00000e+00 +0.00000e+00 +2.00000e+00 +-2.00000e+00 +0.00000e+00 +2.00000e+00 +-1.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +-1.00000e+00 +-1.00000e+00 +0.00000e+00 +1.00000e+00 +1.00000e+00 +0.00000e+00 +-1.00000e+00 +1.00000e+00 +2.00000e+00 +1.00000e+00 +0.00000e+00 +1.00000e+00 +0.00000e+00 +1.00000e+00 +2.00000e+00 +1.00000e+00 +0.00000e+00 +1.00000e+00 +0.00000e+00 +2.50000e+00 +2.50000e+00 +2.50000e+00 +3.50000e+00 +3.50000e+00 +3.50000e+00 +2.50000e+00 +2.50000e+00 +2.50000e+00 +3.50000e+00 +3.50000e+00 +3.50000e+00 +hexa8 + 4 + 7 8 11 10 1 2 5 4 + 8 9 12 11 2 3 6 5 + 13 14 17 18 19 20 23 24 + 14 15 16 17 20 21 22 23 diff --git a/single-node-refactor/contact_test_inputs/edge_case3.yaml b/single-node-refactor/contact_test_inputs/edge_case3.yaml new file mode 100644 index 000000000..7152f5439 --- /dev/null +++ b/single-node-refactor/contact_test_inputs/edge_case3.yaml @@ -0,0 +1,100 @@ +# num_dims: 3 + +dynamic_options: + time_final: 1.0 + dt_min: 0.5 + dt_max: 0.5 + dt_start: 1.e-2 + cycle_stop: 300000 + + +mesh_options: + source: file + file_path: contact_test_inputs/edge_case3.geo + +#mesh_options: +# source: generate +# num_dims: 3 +# type: Box +# origin: [0.0, 0.0, 0.0] +# length: [1.2, 1.2, 1.2] +# num_elems: [12, 12, 12] + + +output_options: + timer_output_level: thorough + output_file_format: ensight + graphics_time_step: 0.5 + # graphics_iteration_step: 10 + +solver_options: + - solver: + method: SGH + # solver_vars: + # - blah + # - blah + # - blah + +boundary_conditions: + + # Tag Z plane + - boundary_condition: + solver: SGH + geometry: z_plane + direction: z_dir + value: 0.0 + type: reflected + + # global contact condition (this should be used when it is not clear what the contact surface(s) will be) + - boundary_condition: + solver: SGH + geometry: global + type: contact + +materials: + - material: + id: 0 + eos_model: ideal_gas + # strength_model: none + elastic_modulus: 10 + poisson_ratio: 0.3 + q1: 1.0 + q2: 1.333 + q1ex: 1.0 + q2ex: 1.333 + eos_global_vars: + - 3 + - 1.0E-14 + - 1.0 + +regions: + - fill_volume: + type: global + material_id: 0 + den: 1.0 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: 0.0 + w: 0.0 + + + + # top block initial velocity conditions + - fill_volume: + type: box + origin: [0.0, 0.0, 0.0] + x1: -0.5 + x2: 0.5 + y1: -1.0 + y2: 1.0 + z1: 2.5 + z2: 3.5 + material_id: 0 + den: 1.0 + # ie: 0.25833839995946534 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: 0.0 + w: -1.0 \ No newline at end of file diff --git a/single-node-refactor/contact_test_inputs/multi_hit_test.yaml b/single-node-refactor/contact_test_inputs/multi_hit_test.yaml new file mode 100644 index 000000000..cf15e4c49 --- /dev/null +++ b/single-node-refactor/contact_test_inputs/multi_hit_test.yaml @@ -0,0 +1,100 @@ +# num_dims: 3 + +dynamic_options: + time_final: 3.0 + dt_min: 1e-6 + dt_max: 0.05 + dt_start: 1.e-5 + cycle_stop: 300000 + + +mesh_options: + source: file + file_path: contact_test_inputs/contact_test.geo + +#mesh_options: +# source: generate +# num_dims: 3 +# type: Box +# origin: [0.0, 0.0, 0.0] +# length: [1.2, 1.2, 1.2] +# num_elems: [12, 12, 12] + + +output_options: + timer_output_level: thorough + output_file_format: ensight + graphics_time_step: 0.02 + # graphics_iteration_step: 10 + +solver_options: + - solver: + method: SGH + # solver_vars: + # - blah + # - blah + # - blah + +boundary_conditions: + + # Tag Z plane + - boundary_condition: + solver: SGH + geometry: z_plane + direction: z_dir + value: 0.0 + type: reflected + + # global contact condition (this should be used when it is not clear what the contact surface(s) will be) + - boundary_condition: + solver: SGH + geometry: global + type: contact + +materials: + - material: + id: 0 + eos_model: ideal_gas + # strength_model: none + elastic_modulus: 10 + poisson_ratio: 0.3 + q1: 1.0 + q2: 1.333 + q1ex: 1.0 + q2ex: 1.333 + eos_global_vars: + - 3 + - 1.0E-14 + - 1.0 + +regions: + - fill_volume: + type: global + material_id: 0 + den: 1.0 + sie: 1.e-10 + velocity: cartesian + u: 0.0 + v: 0.0 + w: 0.0 + + + + # top block initial velocity conditions + - fill_volume: + type: box + origin: [0.0, 0.0, 0.0] + x1: -0.25 + x2: 0.25 + y1: -0.25 + y2: 0.25 + z1: 1 + z2: 1.5 + material_id: 0 + den: 10.0 + # ie: 0.25833839995946534 + sie: 1.e-5 + velocity: cartesian + u: 0.0 + v: 0.0 + w: -0.5 \ No newline at end of file diff --git a/single-node-refactor/src/Solvers/SGH_solver/CMakeLists.txt b/single-node-refactor/src/Solvers/SGH_solver/CMakeLists.txt index 30ac8dcb3..afccd72a1 100755 --- a/single-node-refactor/src/Solvers/SGH_solver/CMakeLists.txt +++ b/single-node-refactor/src/Solvers/SGH_solver/CMakeLists.txt @@ -22,16 +22,15 @@ include_directories(src) set(SRC_Files -src/boundary.cpp +src/boundary.cpp +src/contact.cpp src/energy_sgh.cpp src/force_sgh.cpp src/position.cpp src/momentum.cpp src/properties.cpp src/sgh_solve.cpp -src/time_integration.cpp -src/user_mat.cpp -src/user_mat_init.cpp) +src/time_integration.cpp) # set(SGH_Solver_SRC src/sgh_solver.cpp ) diff --git a/single-node-refactor/src/Solvers/SGH_solver/include/_debug_tools.h b/single-node-refactor/src/Solvers/SGH_solver/include/_debug_tools.h new file mode 100644 index 000000000..6d91ac466 --- /dev/null +++ b/single-node-refactor/src/Solvers/SGH_solver/include/_debug_tools.h @@ -0,0 +1,53 @@ +#pragma once + +#include "matar.h" + +using namespace mtr; + +template +void matar_print(const T &array) +{ + int ord = array.order(); + + if (ord == 1) + { + int n = array.dims(0); + for (int i = 0; i < n; i++) + { + std::cout << array(i) << " "; + } + std::cout << std::endl; + } else if (ord == 2) + { + int m = array.dims(0); + int n = array.dims(1); + for (int i = 0; i < m; i++) + { + std::cout << i << " ---> "; + for (int j = 0; j < n; j++) + { + std::cout << array(i, j) << " "; + } + std::cout << std::endl; + } + } else + { + int l = array.dims(0); + int m = array.dims(1); + int n = array.dims(2); + for (int i = 0; i < l; i++) + { + std::cout << "Level " << i << std::endl; + for (int j = 0; j < m; j++) + { + std::cout << j << " ---> "; + for (int k = 0; k < n; k++) + { + std::cout << array(i, j, k) << " "; + } + std::cout << std::endl; + } + std::cout << std::endl; + } + } +} diff --git a/single-node-refactor/src/Solvers/SGH_solver/include/contact.h b/single-node-refactor/src/Solvers/SGH_solver/include/contact.h new file mode 100644 index 000000000..3398fa8e8 --- /dev/null +++ b/single-node-refactor/src/Solvers/SGH_solver/include/contact.h @@ -0,0 +1,563 @@ +#ifndef CONTACT_H +#define CONTACT_H + +#include "matar.h" +#include "mesh.h" +#include "simulation_parameters.h" +#include "_debug_tools.h" // Remove this file entirely once finished + +using namespace mtr; + +// solving options +static constexpr size_t max_iter = 100; // max number of iterations +static constexpr double tol = 1e-15; // tolerance for the things that are supposed to be zero +static constexpr double edge_tol = 1e-3; // tolerance for edge case solutions (see contact_check for more info) + +struct contact_node_t +{ + size_t gid; // global node id + double mass; // mass of the node + CArrayKokkos pos = CArrayKokkos(3); // position of the node + CArrayKokkos vel = CArrayKokkos(3); // velocity of the node + CArrayKokkos internal_force = CArrayKokkos(3); // any force that is not due to contact + CArrayKokkos contact_force = CArrayKokkos(3); // force due to contact + + contact_node_t(); + + contact_node_t(const ViewCArrayKokkos &pos, const ViewCArrayKokkos &vel, + const ViewCArrayKokkos &internal_force, const ViewCArrayKokkos &contact_force, + const double &mass); +}; + +struct contact_patch_t +{ + size_t gid; // global patch id + size_t lid; // local patch id (local to contact_patches_t::contact_patches); this is needed in contact_pairs_t + CArrayKokkos nodes_gid; // global node ids + CArrayKokkos nodes_obj; // contact node objects + + // Iso-parametric coordinates of the patch nodes (1D array of size mesh.num_nodes_in_patch) + // For a standard linear hex, xi = [-1.0, 1.0, 1.0, -1.0], eta = [-1.0, -1.0, 1.0, 1.0] + // For now, these are the same for all patch objects, but should they be different, then remove static and look to + // contact_patches_t::initialize for how to set these values + CArrayKokkos xi; // xi coordinates + CArrayKokkos eta; // eta coordinates + static size_t num_nodes_in_patch; // number of nodes in the patch (or surface) + static constexpr size_t max_nodes = 4; // max number of nodes in the patch (or surface); for allocating memory at compile time + + // members to be used in find_nodes and capture_box + // expected max number of nodes that could hit a patch + static constexpr size_t max_contacting_nodes_in_patch = 25; + static constexpr size_t max_number_buckets = 1000; // todo: this needs to be ridden of and determined in initialize + // bounds of the capture box (xc_max, yc_max, zc_max, xc_min, yc_min, zc_min) + CArrayKokkos bounds = CArrayKokkos(6); + // buckets that intersect the patch + CArrayKokkos buckets = CArrayKokkos(max_number_buckets); + // nodes that could potentially contact the patch + CArrayKokkos possible_nodes = CArrayKokkos(max_contacting_nodes_in_patch); + + contact_patch_t(); + + contact_patch_t(const ViewCArrayKokkos &points, const ViewCArrayKokkos &vel_points, + const ViewCArrayKokkos &internal_force_points, + const ViewCArrayKokkos &contact_force_points, const ViewCArrayKokkos &mass_points_); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn capture_box + /// + /// \brief Constructs the capture box for the patch + /// + /// The capture box is used to determine which buckets penetrate the surface/patch. The nodes in the intersecting + /// buckets are considered for potential contact. The capture box is constructed from the maximum absolute value + /// of velocity and acceleration by considering the position at time dt, which is equal to + /// + /// position + velocity_max*dt + 0.5*acceleration_max*dt^2 and + /// position - velocity_max*dt - 0.5*acceleration_max*dt^2 + /// + /// The maximum and minimum components of the capture box are recorded and will be used in + /// contact_patches_t::find_nodes. + /// + /// \param vx_max absolute maximum x velocity across all nodes in the patch + /// \param vy_max absolute maximum y velocity across all nodes in the patch + /// \param vz_max absolute maximum z velocity across all nodes in the patch + /// \param ax_max absolute maximum x acceleration across all nodes in the patch + /// \param ay_max absolute maximum y acceleration across all nodes in the patch + /// \param az_max absolute maximum z acceleration across all nodes in the patch + /// \param dt time step + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + void capture_box(const double &vx_max, const double &vy_max, const double &vz_max, + const double &ax_max, const double &ay_max, const double &az_max, + const double &dt); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn get_contact_point + /// + /// \brief Finds the contact point in the reference space with the given contact node + /// + /// The row node_lid of det_sol is taken as the guess which is of the order (xi, eta, del_tc) where del_tc is the + /// time it takes for the node to penetrate the patch/surface. This will iteratively solve using a Newton-Raphson + /// scheme and will change det_sol in place. + /// + /// \param node Contact node object that is potentially penetrating this patch/surface + /// \param xi_val xi value to change in place + /// \param eta_val eta value to change in place + /// \param del_tc del_tc value to change in place + /// + /// \return true if a solution was found in less than max_iter iterations; false if the solution took up to max_iter + /// iterations or if a singularity was encountered + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + KOKKOS_FUNCTION // will be called inside a macro + bool get_contact_point(const contact_node_t &node, double &xi_val, double &eta_val, double &del_tc) const; + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn contact_check + /// + /// \brief Determines if a contact pair should be formed + /// + /// This is responsible for getting a guess value to feed into the get_contact_point method as well as constructing + /// the logic to determine if a contact pair should be formed. If a solution is found from that scheme, the + /// reference coordinates of xi and eta are between -1 and 1, and the calculated del_tc is between 0 and the current + /// time step (del_t), then this will return true and a contact pair should be formed. The exception to not adding a + /// contact pair between 'this' and 'node' is if the solution is on the edge. This behavior is handled in + /// contact_patches_t::get_contact_pairs. + /// + /// As for the significance of the `tol` and `edge_tol` parameters, `tol` is used to determine the convergence of + /// the Newton-Raphson scheme as well as the edge case for the time bound. The del_tc value is then considered true + /// at `0 - tol` and `del_t + tol`. Similarly, `edge_tol` is used for the solution of `xi` and `eta` to determine if + /// the solution is within the bounds of the patch/surface. If the absolute value of `xi` and `eta` are less than + /// or equal to `1 + edge_tol`, then the node is passing through the patch/surface. + /// + /// \param node Contact node object that is being checked for contact with 'this' + /// \param del_t time step + /// \param xi_val xi value to change in place + /// \param eta_val eta value to change in place + /// \param del_tc del_tc value to change in place + /// + /// \return true if a contact pair should be formed; false otherwise + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + KOKKOS_FUNCTION + bool contact_check(const contact_node_t &node, const double &del_t, double &xi_val, double &eta_val, + double &del_tc) const; + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn construct_basis + /// + /// \brief Constructs the basis matrix for the patch + /// + /// The columns of A are defined as the position of each patch/surface node at time del_t. This is a 3x4 matrix for + /// a standard linear hex element and its columns are constructed like this: + /// + /// ⎡p_{nx} + v_{nx}*del_t + 0.5*a_{nx}*del_t^2 ... for each n⎤ + /// ⎢ ⎥ + /// ⎡p_{ny} + v_{ny}*del_t + 0.5*a_{ny}*del_t^2 ... for each n⎤ + /// ⎢ ⎥ + /// ⎣p_{nz} + v_{nz}*del_t + 0.5*a_{nz}*del_t^2 ... for each n⎦ + /// + /// \param A basis matrix as defined above (will be changed in place) + /// \param del_t time step to construct the basis matrix + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + KOKKOS_FUNCTION + void construct_basis(ViewCArrayKokkos &A, const double &del_t) const; + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn ref_to_physical + /// + /// \brief Converts the reference coordinates to physical coordinates + /// + /// This method will convert the reference coordinates defined by 'this->xi' and 'this->eta' to the physical/global + /// coordinates. + /// + /// \param ref 1D reference coordinates (xi, eta) + /// \param A basis matrix as defined in construct_basis + /// \param phys 1D physical coordinates (x, y, z) + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + KOKKOS_FUNCTION + void ref_to_physical(const ViewCArrayKokkos &ref, const ViewCArrayKokkos &A, + ViewCArrayKokkos &phys) const; + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn phi + /// + /// \brief Modifies the phi_k array to contain the basis function values at the given xi and eta values + /// + /// \param phi_k basis function values that correspond to the `this->xi` and `this->eta` values + /// \param xi_value xi value + /// \param eta_value eta value + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + KOKKOS_FUNCTION + void phi(ViewCArrayKokkos &phi_k, const double &xi_value, const double &eta_value) const; + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn d_phi_d_xi + /// + /// \brief Modifies the d_phi_k_d_xi array to contain the basis function derivatives with respect to xi at the given + /// xi and eta values + /// + /// \param d_phi_k_d_xi basis function values that correspond to the `this->xi` and `this->eta` values + /// \param xi_value xi value + /// \param eta_value eta value + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + KOKKOS_FUNCTION + void d_phi_d_xi(ViewCArrayKokkos &d_phi_k_d_xi, const double &xi_value, const double &eta_value) const; + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn d_phi_d_eta + /// + /// \brief Modifies the d_phi_k_d_eta array to contain the basis function derivatives with respect to eta at the + /// given xi and eta values + /// + /// \param d_phi_k_d_eta basis function values that correspond to the `this->xi` and `this->eta` values + /// \param xi_value xi value + /// \param eta_value eta value + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + KOKKOS_FUNCTION + void d_phi_d_eta(ViewCArrayKokkos &d_phi_k_d_eta, const double &xi_value, const double &eta_value) const; + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn get_normal + /// + /// \brief Computes the normal vector of the patch/surface at the given xi and eta values + /// + /// \param xi_val xi value + /// \param eta_val eta value + /// \param del_t time step to compute the normal at + /// \param normal kokkos view that will be changed in place + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + KOKKOS_FUNCTION + void get_normal(const double &xi_val, const double &eta_val, const double &del_t, + ViewCArrayKokkos &normal) const; +}; + +// forward declaration +struct contact_patches_t; + +struct contact_pair_t +{ + contact_patch_t patch; // patch object (or surface) + contact_node_t node; // node object + double xi; // xi coordinate of the contact point + double eta; // eta coordinate of the contact point + double del_tc; // time it takes for the node to penetrate the patch/surface (only useful for initial contact) + CArrayKokkos normal = CArrayKokkos(3); // normal vector of the patch/surface at the contact point + + bool active = false; // if the pair is active or not + + // force members + double fc_inc = 0.0; // force increment to be added to contact_node_t::contact_force + double fc_inc_total = 0.0; // all previous force increments get summed to this member (see contact_patches_t::force_resolution()) + + enum contact_types + { + frictionless, // no friction; only normal force + glue // contact point stays constant + }; + + // todo: frictionless is the only contact type implemented, but in the future, changing this member before the + // force resolution call will allow for different contact types + contact_types contact_type = frictionless; // default contact type + + contact_pair_t(); + + KOKKOS_FUNCTION + contact_pair_t(contact_patches_t &contact_patches_obj, const contact_patch_t &patch_obj, + const contact_node_t &node_obj, const double &xi_val, const double &eta_val, + const double &del_tc_val, const ViewCArrayKokkos &normal_view); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn frictionless_increment + /// + /// \brief Computes the force increment for the contact pair with no friction + /// + /// This method will compute the force increment between a contact patch and node with no friction. The force is + /// strictly calculated in the normal direction only. The xi, eta, and fc_inc members will be changed in place. The + /// force increment value is determined by kinematic conditions and will result in the position of the node being + /// on the patch/surface at time del_t. + /// + /// \param contact_patches contact_patches object + /// \param del_t current time step in the analysis + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + KOKKOS_FUNCTION + void frictionless_increment(const double &del_t); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn distribute_frictionless_force + /// + /// \brief Distributes the force increment to the penetrating node and the patch nodes + /// + /// This method will distribute an incremental force value (fc_inc member) to the penetrating node and the patch + /// nodes. For the penetrating node, it's N*fc_inc and for the patch nodes, a value of -N*fc_inc*phi_k where N is + /// the unit normal and phi_k is the basis function array values as defined in contact_patch_t::phi. This method + /// will also add the force increment to fc_inc_total. If fc_inc_total is less than zero, then this means that a + /// tensile force is required to keep the node on the patch, but this is not possible since contact is always + /// compressive when there is no adhesive phenomena. When fc_inc_total goes below zero, fc_inc is set to zero, + /// and the left over fc_inc_total will be subtracted from the penetrating node and the patch nodes, then + /// fc_inc_total is set to zero. + /// + /// \param force_scale instead of distributing the full fc_inc, a fraction of it can be distributed to prevent large + /// shocks to the solving scheme + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + KOKKOS_FUNCTION + void distribute_frictionless_force(const double &force_scale); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn should_remove + /// + /// \brief Determines if the contact pair should be removed + /// + /// This method will determine if the contact pair should be removed. If the fc_inc_total value is zero or if the + /// xi and eta values are out of the bounds of the patch (not between -1 - edge_tol and 1 + edge_tol), then the pair + /// will be removed. This method is not used for all contact types as some (i.e. glue) require different conditions. + /// Additionally, this method will update the unit normal to the current xi and eta values. At the moment, this + /// method is used for frictionless contact only. + /// + /// \param del_t current time step in the analysis + /// + /// \return true if the contact pair should be removed; false otherwise + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + bool should_remove(const double &del_t); +}; + +struct contact_patches_t +{ + CArrayKokkos contact_patches; // patches that will be checked for contact + CArrayKokkos contact_nodes; // all nodes that are in contact_patches (accessed through node gid) + CArrayKokkos patches_gid; // global patch ids + CArrayKokkos nodes_gid; // global node ids + size_t num_contact_patches; // total number of patches that will be checked for contact + RaggedRightArrayKokkos patches_in_node; // each row is the node gid and the columns are the patches (local to contact_patches) that the node is in + CArrayKokkos num_patches_in_node; // the strides for patches_in_node + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn initialize + /// + /// \brief Initializes the contact_patches array + /// + /// \param mesh mesh object + /// \param bdy_contact_patches global ids of patches that will be checked for contact + /// \param nodes node object that contains coordinates and velocities of all nodes + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + void initialize(const mesh_t &mesh, const CArrayKokkos &bdy_contact_patches, const node_t &nodes, + const corner_t &corners); + + /* + * Here is a description of each array below: + * + * nbox: An array consisting of the number of nodes that a bucket contains. For example, nbox[0] = 2 would mean + * that bucket 0 has 2 nodes. + * lbox: An array consisting of the bucket id for each node. For example, lbox[5] = 12 would mean that node 5 is in + * bucket 12. + * nsort: An array consisting of sorted nodes based off the bucket id. + * npoint: An array consisting of indices into nsort where each index is the starting location. For example, + * nsort[npoint[5]] would return the starting node in bucket 5. This is a means for finding the nodes given + * a bucket id. + * + * With the above data structure, you could easily get the nodes in a bucket by the following: + * nsort[npoint[bucket_id]:npoint[bucket_id] + nbox[bucket_id]] + * + * Buckets are ordered by propagating first in the x direction, then in the y direction, and finally in the z. + */ + CArrayKokkos nbox; // Size nb buckets + CArrayKokkos lbox; // Size num_contact_nodes nodes (num_contact_nodes is the total number of nodes being checked for penetration) + CArrayKokkos nsort; // Size num_contact_nodes nodes + CArrayKokkos npoint; // Size nb buckets + + static double bucket_size; // bucket size (defined as 1.001*min_node_distance) + static size_t num_contact_nodes; // total number of contact nodes (always less than or equal to mesh.num_bdy_nodes) + double x_max = 0.0; // maximum x coordinate + double y_max = 0.0; // maximum y coordinate + double z_max = 0.0; // maximum z coordinate + double x_min = 0.0; // minimum x coordinate + double y_min = 0.0; // minimum y coordinate + double z_min = 0.0; // minimum z coordinate + double vx_max = 0.0; // maximum x velocity + double vy_max = 0.0; // maximum y velocity + double vz_max = 0.0; // maximum z velocity + double ax_max = 0.0; // maximum x acceleration + double ay_max = 0.0; // maximum y acceleration + double az_max = 0.0; // maximum z acceleration + size_t Sx = 0; // number of buckets in the x direction + size_t Sy = 0; // number of buckets in the y direction + size_t Sz = 0; // number of buckets in the z direction + + CArrayKokkos contact_pairs; // contact pairs (accessed through node gid) + DynamicRaggedRightArrayKokkos contact_pairs_access; // each row is the patch gid and the columns represent the node in contact with the patch; used for quick access and iterating + CArrayKokkos is_patch_node; // container for determining if a node is a patch node for a contact pair + CArrayKokkos is_pen_node; // container for determining if a node is a penetrating node for a contact pair + CArrayKokkos active_pairs; // array of only the active pairs (accessed through node gid) + CArrayKokkos forces; // member to store contact force increments (only used to check convergence) + size_t num_active_pairs = 0; // number of active pairs + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn sort + /// + /// \brief Constructs nbox, lbox, nsort, and npoint according to the Sandia Algorithm + /// + /// \param mesh mesh object + /// \param nodes node object that contains coordinates and velocities of all nodes + /// \param corner corner object that contains corner forces + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + void sort(); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn update_nodes + /// + /// \brief Updates the coordinates, velocities, internal forces, and zeros contact force for all contact nodes + /// + /// \param mesh mesh object + /// \param nodes node object that contains coordinates and velocities of all nodes + /// \param corner corner object that contains corner forces + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + void update_nodes(const mesh_t &mesh, const node_t &nodes, const corner_t &corner); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn find_nodes + /// + /// \brief Finds the nodes that could potentially contact a surface/patch + /// + /// \param contact_patch patch object of interest + /// \param del_t current time step in the analysis + /// \param num_nodes_found number of nodes that could potentially contact the patch (used to access possible_nodes) + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + void find_nodes(contact_patch_t &contact_patch, const double &del_t, size_t &num_nodes_found); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn get_contact_pairs + /// + /// \brief Constructs the contact pairs + /// + /// This will construct this->contact_pairs and this->contact_pairs_access. This member will be called once before + /// force resolution, then it will be called iteratively up to a certain max or until no new contact pairs are found + /// after the force resolution. The algorithm presented in this member does not have a master and slave hierarchy, + /// and the pairs are determined by whichever node is penetrating first. An important characteristic of the + /// datastructure is that a contact patch can have multiple nodes, but a contact node is only associated with one + /// contact patch. + /// + /// \param del_t current time step in the analysis + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + void get_contact_pairs(const double &del_t); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn remove_pair + /// + /// \brief Removes a contact pair from the contact_pairs_access array + /// + /// \param pair Contact pair object to remove + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + KOKKOS_FUNCTION + void remove_pair(contact_pair_t &pair); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn remove_pairs + /// + /// \brief Loops through all active pairs and removes the pairs that don't meet the criteria + /// + /// This method will walk through all the active contact pairs and remove the pairs that don't meet the criteria of + /// the corresponding `should_remove` function. + /// + /// \param del_t current time step in the analysis + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + void remove_pairs(const double &del_t); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn get_edge_pair + /// + /// \brief Determines the more dominant pair for the case when a penetrating node contacts an edge + /// + /// This method will determine the best pair to use based off the most opposing normal. For each pair, the normal + /// at the contact point is dotted with the surface normal of the penetrating node. The most negative dot product + /// value indicates the superior patch to pair to. If the dot products are the same, then the second patch is used + /// with a normal being the average of the two. + /// + /// \param normal1 normal of the already existing pair + /// \param normal2 normal of the current pair in the iterations + /// \param node_gid global node id of the penetrating node + /// \param del_t current time step in the analysis + /// \param new_normal modified normal to be used in the case where a new pair should be added + /// + /// \return true if a new pair should be added; false otherwise + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + KOKKOS_FUNCTION + bool get_edge_pair(const ViewCArrayKokkos &normal1, const ViewCArrayKokkos &normal2, + const size_t &node_gid, const double &del_t, ViewCArrayKokkos &new_normal) const; + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// \fn force_resolution + /// + /// \brief Resolves the contact forces + /// + /// todo: add more information here + /// + /// \param del_t current time step in the analysis + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + void force_resolution(const double &del_t); +}; + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +/// \fn mat_mul +/// +/// \brief Matrix multiplication with A*x = b +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +void mat_mul(const ViewCArrayKokkos &A, const ViewCArrayKokkos &x, ViewCArrayKokkos &b); + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +/// \fn norm +/// +/// \brief Computes the norm (sqrt(x1^2 + x2^2 + ...)) of a 1D array +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +double norm(const ViewCArrayKokkos &x); + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +/// \fn det +/// +/// \brief Finds the determinant of a 3x3 matrix +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +KOKKOS_INLINE_FUNCTION +double det(const ViewCArrayKokkos &A); + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +/// \fn inv +/// +/// \brief Finds the inverse of a 3x3 matrix +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +void inv(const ViewCArrayKokkos &A, ViewCArrayKokkos &A_inv, const double &A_det); + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +/// \fn dot +/// +/// \brief Computes the dot product of two 1D arrays +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +double dot(const ViewCArrayKokkos &a, const ViewCArrayKokkos &b); + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +/// \fn outer +/// +/// \brief Computes the outer product of two 1D arrays +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +void outer(const ViewCArrayKokkos &a, const ViewCArrayKokkos &b, ViewCArrayKokkos &c); + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +/// \fn all +/// +/// \brief Checks if all elements in the array are true up to the provided size +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +bool all(const ViewCArrayKokkos &a, const size_t &size); + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +/// \fn any +/// +/// \brief Checks if any elements in the array are true up to the provided size +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +KOKKOS_FUNCTION +bool any(const ViewCArrayKokkos &a, const size_t &size); + +// run tests +void run_contact_tests(contact_patches_t &contact_patches_obj, const mesh_t &mesh, const node_t &nodes, + const corner_t &corner, const simulation_parameters_t &sim_params); + +#endif // CONTACT_H \ No newline at end of file diff --git a/single-node-refactor/src/Solvers/SGH_solver/include/sgh_solver.h b/single-node-refactor/src/Solvers/SGH_solver/include/sgh_solver.h index a2124f656..9f0735316 100644 --- a/single-node-refactor/src/Solvers/SGH_solver/include/sgh_solver.h +++ b/single-node-refactor/src/Solvers/SGH_solver/include/sgh_solver.h @@ -38,6 +38,7 @@ #include "matar.h" #include "solver.h" #include "geometry_new.h" +#include "contact.h" // #include "io_utils.h" #include "simulation_parameters.h" @@ -75,6 +76,9 @@ class SGH : public Solver int rk_num_stages = 2; int cycle_stop = 1000000000; + contact_patches_t contact_bank; // keeps track of contact patches + bool doing_contact = false; // Condition used in SGH::execute + SGH() : Solver() { } @@ -113,8 +117,8 @@ class SGH : public Solver dt = sim_param.dynamic_options.dt_start; graphics_id = 0; - graphics_times(0) = 0.0; - graphics_time = 0.0; // the times for writing graphics dump + graphics_times(0) = sim_param.output_options.graphics_time_step; + graphics_time = sim_param.output_options.graphics_time_step; } ///////////////////////////////////////////////////////////////////////////// @@ -130,6 +134,23 @@ class SGH : public Solver std::cout << "Applying initial boundary conditions" << std::endl; boundary_velocity(mesh, sim_param.boundary_conditions, node.vel, time_value); + + // Setting up contact + for (size_t i = 0; i < mesh.num_bdy_sets; i++) { + boundary_condition_t bound = sim_param.boundary_conditions(i); + if (bound.type == boundary_conds::contact && bound.geometry == boundary_conds::global) { + std::cout << "Setting up global contact" << std::endl; + doing_contact = true; + + contact_bank.initialize(mesh, mesh.bdy_patches, node, corner); + // run_contact_tests(contact_bank, mesh, node, corner, sim_param); + break; + } else if (bound.type == boundary_conds:: contact && bound.geometry != boundary_conds::global) { + doing_contact = true; + std::cerr << "Contact boundary conditions are only supported for global at the moment." << std::endl; + exit(1); + } + } } ///////////////////////////////////////////////////////////////////////////// @@ -151,15 +172,11 @@ class SGH : public Solver // **** Functions defined in boundary.cpp **** // void boundary_velocity( const mesh_t& mesh, - const CArrayKokkos& boundary, + const DCArrayKokkos& boundary, DCArrayKokkos& node_vel, const double time_value); - void boundary_contact( - const mesh_t& mesh, - const CArrayKokkos& boundary, - DCArrayKokkos& node_vel, - const double time_value); + void boundary_contact(const double &del_t); // **** Functions defined in energy_sgh.cpp **** // void update_energy( @@ -174,7 +191,7 @@ class SGH : public Solver // **** Functions defined in force_sgh.cpp **** // void get_force( - const CArrayKokkos& material, + const DCArrayKokkos& material, const mesh_t& mesh, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, @@ -194,7 +211,7 @@ class SGH : public Solver const double rk_alpha); void get_force_2D( - const CArrayKokkos& material, + const DCArrayKokkos& material, const mesh_t& mesh, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, @@ -229,7 +246,8 @@ class SGH : public Solver const mesh_t& mesh, DCArrayKokkos& node_vel, const DCArrayKokkos& node_mass, - const DCArrayKokkos& corner_force); + const DCArrayKokkos& corner_force, + const CArrayKokkos &contact_nodes); KOKKOS_FUNCTION void get_velgrad( @@ -277,7 +295,7 @@ class SGH : public Solver // **** Functions defined in properties.cpp **** // void update_state( - const CArrayKokkos& material, + const DCArrayKokkos& material, const mesh_t& mesh, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, @@ -294,7 +312,7 @@ class SGH : public Solver const double rk_alpha); void update_state2D( - const CArrayKokkos& material, + const DCArrayKokkos& material, const mesh_t& mesh, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, @@ -351,36 +369,6 @@ class SGH : public Solver double& dt, const double fuzz); - // **** Functions defined in user_mat.cpp **** // - // NOTE: Pull up into high level - KOKKOS_FUNCTION - void user_eos_model( - const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const size_t elem_gid, - const size_t mat_id, - const DCArrayKokkos& elem_state_vars, - const DCArrayKokkos& elem_sspd, - const double den, - const double sie); - - KOKKOS_FUNCTION - void user_strength_model( - const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const size_t elem_gid, - const size_t mat_id, - const DCArrayKokkos& elem_state_vars, - const DCArrayKokkos& elem_sspd, - const double den, - const double sie, - const ViewCArrayKokkos& vel_grad, - const ViewCArrayKokkos& elem_node_gids, - const DCArrayKokkos& node_coords, - const DCArrayKokkos& node_vel, - const double vol, - const double dt, - const double rk_alpha); }; #endif // end HEADER_H diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/boundary.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/boundary.cpp index 6d557a819..1c784c54f 100644 --- a/single-node-refactor/src/Solvers/SGH_solver/src/boundary.cpp +++ b/single-node-refactor/src/Solvers/SGH_solver/src/boundary.cpp @@ -1,5 +1,5 @@ /********************************************************************************************** - © 2020. Triad National Security, LLC. All rights reserved. + � 2020. Triad National Security, LLC. All rights reserved. This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. Department of Energy/National Nuclear Security Administration. All rights in the program are @@ -47,7 +47,7 @@ /// ///////////////////////////////////////////////////////////////////////////// void SGH::boundary_velocity(const mesh_t& mesh, - const CArrayKokkos& boundary, + const DCArrayKokkos& boundary, DCArrayKokkos& node_vel, const double time_value) { @@ -116,12 +116,10 @@ void SGH::boundary_velocity(const mesh_t& mesh, /// \param The current simulation time /// ///////////////////////////////////////////////////////////////////////////// -void SGH::boundary_contact(const mesh_t& mesh, - const CArrayKokkos& boundary, - DCArrayKokkos& node_vel, - const double time_value) +void SGH::boundary_contact(const double &del_t) { - - - return; + contact_bank.sort(); + contact_bank.get_contact_pairs(del_t); + contact_bank.force_resolution(del_t); + contact_bank.remove_pairs(del_t); } // end boundary_contact function \ No newline at end of file diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/contact.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/contact.cpp new file mode 100644 index 000000000..91747ab50 --- /dev/null +++ b/single-node-refactor/src/Solvers/SGH_solver/src/contact.cpp @@ -0,0 +1,2034 @@ +#include "contact.h" + +// Definition of static member variables +size_t contact_patch_t::num_nodes_in_patch; +double contact_patches_t::bucket_size; +size_t contact_patches_t::num_contact_nodes; + +/// beginning of global, linear algebra functions ////////////////////////////////////////////////////////////////////// +#pragma clang diagnostic push +#pragma ide diagnostic ignored "OCDFAInspection" +KOKKOS_FUNCTION +void mat_mul(const ViewCArrayKokkos &A, const ViewCArrayKokkos &x, ViewCArrayKokkos &b) +{ + size_t x_ord = x.order(); + size_t m = A.dims(0); + size_t n = A.dims(1); + + if (x_ord == 1) + { + for (size_t i = 0; i < m; i++) + { + b(i) = 0.0; + for (size_t k = 0; k < n; k++) + { + b(i) += A(i, k)*x(k); + } + } + } else + { + size_t p = x.dims(1); + for (size_t i = 0; i < m; i++) + { + for (size_t j = 0; j < p; j++) + { + b(i, j) = 0.0; + for (size_t k = 0; k < n; k++) + { + b(i, j) += A(i, k)*x(k, j); + } + } + } + } +} // end mat_mul +#pragma clang diagnostic pop + +KOKKOS_FUNCTION +double norm(const ViewCArrayKokkos &x) +{ + double sum = 0.0; + for (size_t i = 0; i < x.size(); i++) + { + sum += pow(x(i), 2); + } + return sqrt(sum); +} // end norm + +KOKKOS_INLINE_FUNCTION +double det(const ViewCArrayKokkos &A) +{ + return A(0, 0)*(A(1, 1)*A(2, 2) - A(1, 2)*A(2, 1)) - A(0, 1)*(A(1, 0)*A(2, 2) - A(1, 2)*A(2, 0)) + + A(0, 2)*(A(1, 0)*A(2, 1) - A(1, 1)*A(2, 0)); +} // end det + +KOKKOS_FUNCTION +void inv(const ViewCArrayKokkos &A, ViewCArrayKokkos &A_inv, const double &A_det) +{ + // A_inv = 1/det(A)*adj(A) + A_inv(0, 0) = (A(1, 1)*A(2, 2) - A(1, 2)*A(2, 1))/A_det; + A_inv(0, 1) = (A(0, 2)*A(2, 1) - A(0, 1)*A(2, 2))/A_det; + A_inv(0, 2) = (A(0, 1)*A(1, 2) - A(0, 2)*A(1, 1))/A_det; + A_inv(1, 0) = (A(1, 2)*A(2, 0) - A(1, 0)*A(2, 2))/A_det; + A_inv(1, 1) = (A(0, 0)*A(2, 2) - A(0, 2)*A(2, 0))/A_det; + A_inv(1, 2) = (A(0, 2)*A(1, 0) - A(0, 0)*A(1, 2))/A_det; + A_inv(2, 0) = (A(1, 0)*A(2, 1) - A(1, 1)*A(2, 0))/A_det; + A_inv(2, 1) = (A(0, 1)*A(2, 0) - A(0, 0)*A(2, 1))/A_det; + A_inv(2, 2) = (A(0, 0)*A(1, 1) - A(0, 1)*A(1, 0))/A_det; +} // end inv + +KOKKOS_FUNCTION +double dot(const ViewCArrayKokkos &a, const ViewCArrayKokkos &b) +{ + double sum = 0.0; + for (size_t i = 0; i < a.size(); i++) + { + sum += a(i)*b(i); + } + return sum; +} // end dot + +KOKKOS_FUNCTION +void outer(const ViewCArrayKokkos &a, const ViewCArrayKokkos &b, ViewCArrayKokkos &c) +{ + for (size_t i = 0; i < a.size(); i++) + { + for (size_t j = 0; j < b.size(); j++) + { + c(i, j) = a(i)*b(j); + } + } +} // end outer + +KOKKOS_FUNCTION +bool all(const ViewCArrayKokkos &a, const size_t &size) +{ + for (size_t i = 0; i < size; i++) + { + if (!a(i)) + { + return false; + } + } + return true; +} // end all + +KOKKOS_FUNCTION +bool any(const ViewCArrayKokkos &a, const size_t &size) +{ + for (size_t i = 0; i < size; i++) + { + if (a(i)) + { + return true; + } + } + return false; +} // end any +/// end of global, linear algebra functions //////////////////////////////////////////////////////////////////////////// + +/// beginning of contact_patch_t member functions ////////////////////////////////////////////////////////////////////// +void contact_patch_t::capture_box(const double &vx_max, const double &vy_max, const double &vz_max, + const double &ax_max, const double &ay_max, const double &az_max, + const double &dt) +{ + // collecting all the points that will be used to construct the bounding box into add_sub + // the bounding box is the maximum and minimum points for each dimension + // todo: add_sub needs to be moved to a member variable. Can't have this being allocated for every call. + CArrayKokkos add_sub(2, 3, contact_patch_t::num_nodes_in_patch); + double max_v_arr[3] = {vx_max, vy_max, vz_max}; + double max_a_arr[3] = {ax_max, ay_max, az_max}; + FOR_ALL_CLASS(i, 0, contact_patch_t::num_nodes_in_patch, { + const contact_node_t &node_obj = nodes_obj(i); + for (int j = 0; j < 3; j++) + { + add_sub(0, j, i) = node_obj.pos(j) + max_v_arr[j]*dt + 0.5*max_a_arr[j]*dt*dt; + add_sub(1, j, i) = node_obj.pos(j) - max_v_arr[j]*dt - 0.5*max_a_arr[j]*dt*dt; + } + }); + Kokkos::fence(); + + // todo: Does bounds(i) = result_max and bounds(i + 3) = result_min need to be inside Run()? + for (int i = 0; i < 3; i++) + { + // Find the max of dim i + double local_max; + double result_max; + REDUCE_MAX(j, 0, contact_patch_t::num_nodes_in_patch, local_max, { + if (local_max < add_sub(0, i, j)) + { + local_max = add_sub(0, i, j); + } + }, result_max); + bounds(i) = result_max; + + // Find the min of dim i + double local_min; + double result_min; + REDUCE_MIN(j, 0, contact_patch_t::num_nodes_in_patch, local_min, { + if (local_min > add_sub(1, i, j)) + { + local_min = add_sub(1, i, j); + } + }, result_min); + bounds(i + 3) = result_min; + } + Kokkos::fence(); +} // end capture_box + +KOKKOS_FUNCTION +void contact_patch_t::construct_basis(ViewCArrayKokkos &A, const double &del_t) const +{ + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < num_nodes_in_patch; j++) + { + const contact_node_t &node_obj = nodes_obj(j); + double a = (node_obj.internal_force(i) + node_obj.contact_force(i))/node_obj.mass; + A(i, j) = node_obj.pos(i) + node_obj.vel(i)*del_t + 0.5*a*del_t*del_t; + } + } +} // end construct_basis + +KOKKOS_FUNCTION // will be called inside a macro +bool contact_patch_t::get_contact_point(const contact_node_t &node, double &xi_val, double &eta_val, + double &del_tc) const +{ + // In order to understand this, just see this PDF: + // https://github.com/gabemorris12/contact_surfaces/blob/master/Finding%20the%20Contact%20Point.pdf + + // The python version of this is also found in contact.py from that same repo. + + // Using "max_nodes" for the array size to ensure that the array is large enough + double A_arr[3*contact_patch_t::max_nodes]; + ViewCArrayKokkos A(&A_arr[0], 3, num_nodes_in_patch); + + double phi_k_arr[contact_patch_t::max_nodes]; + ViewCArrayKokkos phi_k(&phi_k_arr[0], num_nodes_in_patch); + + double d_phi_d_xi_arr[contact_patch_t::max_nodes]; + ViewCArrayKokkos d_phi_d_xi_(&d_phi_d_xi_arr[0], num_nodes_in_patch); + + double d_phi_d_eta_arr[contact_patch_t::max_nodes]; + ViewCArrayKokkos d_phi_d_eta_(&d_phi_d_eta_arr[0], num_nodes_in_patch); + + double d_A_d_del_t_arr[3*contact_patch_t::max_nodes]; + ViewCArrayKokkos d_A_d_del_t(&d_A_d_del_t_arr[0], 3, num_nodes_in_patch); + + double rhs_arr[3]; // right hand side (A_arr*phi_k) + ViewCArrayKokkos rhs(&rhs_arr[0], 3); + + double lhs; // left hand side (node.pos + node.vel*del_t + 0.5*node.acc*del_t*del_t) + + double F_arr[3]; // defined as rhs - lhs + ViewCArrayKokkos F(&F_arr[0], 3); + + double J0_arr[3]; // column 1 of jacobian + ViewCArrayKokkos J0(&J0_arr[0], 3); + + double J1_arr[3]; // column 2 of jacobian + ViewCArrayKokkos J1(&J1_arr[0], 3); + + double J2_arr[3]; // column 3 of jacobian + ViewCArrayKokkos J2(&J2_arr[0], 3); + + double J_arr[9]; // jacobian + ViewCArrayKokkos J(&J_arr[0], 3, 3); + + double J_inv_arr[9]; // inverse of jacobian + ViewCArrayKokkos J_inv(&J_inv_arr[0], 3, 3); + + double J_det; // determinant of jacobian + double sol[3]; // solution containing (xi, eta, del_tc) + sol[0] = xi_val; + sol[1] = eta_val; + sol[2] = del_tc; + + double grad_arr[3]; // J_inv*F term + ViewCArrayKokkos grad(&grad_arr[0], 3); + + // begin Newton-Rasphson solver + size_t iters; // need to keep track of the number of iterations outside the scope of the loop + for (int i = 0; i < max_iter; i++) + { + iters = i; + construct_basis(A, del_tc); + phi(phi_k, xi_val, eta_val); + mat_mul(A, phi_k, rhs); + for (int j = 0; j < 3; j++) + { + double a = (node.internal_force(j) + node.contact_force(j))/node.mass; + lhs = node.pos(j) + node.vel(j)*del_tc + 0.5*a*del_tc*del_tc; + F(j) = rhs(j) - lhs; + } + + if (norm(F) <= tol) + { + break; + } + + d_phi_d_xi(d_phi_d_xi_, xi_val, eta_val); + d_phi_d_eta(d_phi_d_eta_, xi_val, eta_val); + + // Construct d_A_d_del_t + for (int j = 0; j < 3; j++) + { + for (int k = 0; k < num_nodes_in_patch; k++) + { + const contact_node_t &patch_node = nodes_obj(k); + double a = (patch_node.internal_force(j) + patch_node.contact_force(j))/patch_node.mass; + d_A_d_del_t(j, k) = patch_node.vel(j) + a*del_tc; + } + } + + mat_mul(A, d_phi_d_xi_, J0); + mat_mul(A, d_phi_d_eta_, J1); + mat_mul(d_A_d_del_t, phi_k, J2); + // lhs is a function of del_t, so have to subtract the derivative of lhs wrt del_t + for (int j = 0; j < 3; j++) + { + double a = (node.internal_force(j) + node.contact_force(j))/node.mass; + J2(j) = J2(j) - node.vel(j) - a*del_tc; + } + + // Construct the Jacobian + for (int j = 0; j < 3; j++) + { + J(j, 0) = J0(j); + J(j, 1) = J1(j); + J(j, 2) = J2(j); + } + + // Construct the inverse of the Jacobian and check for singularity. Singularities occur when the node and patch + // are travelling at the same velocity (same direction) or if the patch is perfectly planar and the node travels + // parallel to the plane. Both cases mean no force resolution is needed and will return a false condition. These + // conditions can be seen in singularity_detection_check.py in the python version. + J_det = det(J); + if (fabs(J_det) < tol) + { + return false; + } + + inv(J, J_inv, J_det); + mat_mul(J_inv, F, grad); + for (int j = 0; j < 3; j++) + { + sol[j] = sol[j] - grad(j); + } + // update xi, eta, and del_tc + xi_val = sol[0]; + eta_val = sol[1]; + del_tc = sol[2]; + } // end solver loop + + if (iters == max_iter - 1) + { + return false; + } else + { + return true; + } +} // end get_contact_point + +KOKKOS_FUNCTION +bool contact_patch_t::contact_check(const contact_node_t &node, const double &del_t, double &xi_val, double &eta_val, + double &del_tc) const +{ + // Constructing the guess value + // The guess is determined by projecting the node onto a plane formed by the patch at del_t/2. + // First, compute the centroid of the patch at time del_t/2 + double A_arr[3*contact_patch_t::max_nodes]; + ViewCArrayKokkos A(&A_arr[0], 3, num_nodes_in_patch); + construct_basis(A, del_t/2); + + double centroid[3]; + for (int i = 0; i < 3; i++) + { + centroid[i] = 0.0; + for (int j = 0; j < num_nodes_in_patch; j++) + { + centroid[i] += A(i, j); + } + centroid[i] /= num_nodes_in_patch; + } + + // Compute the position of the penetrating node at del_t/2 + double node_later[3]; + for (int i = 0; i < 3; i++) + { + double a = (node.internal_force(i) + node.contact_force(i))/node.mass; + node_later[i] = node.pos(i) + node.vel(i)*del_t/2 + 0.25*a*del_t*del_t; + } + + // Construct the basis vectors. The first row of the matrix is the vector from the centroid to the reference point + // (1, 0) on the patch. The second row of the matrix is the vector from the centroid to the reference point (0, 1). + // The basis matrix is a 2x3 matrix always. + double b1_arr[3]; + ViewCArrayKokkos b1(&b1_arr[0], 3); + double b2_arr[3]; + ViewCArrayKokkos b2(&b2_arr[0], 3); + double p1_arr[2]; + ViewCArrayKokkos p1(&p1_arr[0], 2); + p1(0) = 1.0; + p1(1) = 0.0; + double p2_arr[2]; + ViewCArrayKokkos p2(&p2_arr[0], 2); + p2(0) = 0.0; + p2(1) = 1.0; + ref_to_physical(p1, A, b1); + ref_to_physical(p2, A, b2); + + // Get b1, b2, and node_later relative to centroid + ViewCArrayKokkos v(&node_later[0], 3); // position of node relative to centroid + for (int i = 0; i < 3; i++) + { + b1(i) -= centroid[i]; + b2(i) -= centroid[i]; + v(i) -= centroid[i]; + } + + // b1 and b2 need to be unit vectors to ensure that the guess values are between -1 and 1. + double b1_norm = norm(b1); + double b2_norm = norm(b2); + for (int i = 0; i < 3; i++) + { + b1(i) /= b1_norm; + b2(i) /= b2_norm; + } + // v also needs to be a normal vector, but if its norm is zero, then we leave it as is. + double v_norm = norm(v); + if (v_norm != 0.0) + { + for (int i = 0; i < 3; i++) + { + v(i) /= v_norm; + } + } + + // Get A_basis, which is the basis vectors found above. + double A_basis_arr[2*3]; + ViewCArrayKokkos A_basis(&A_basis_arr[0], 2, 3); + for (int i = 0; i < 3; i++) + { + A_basis(0, i) = b1(i); + A_basis(1, i) = b2(i); + } + + // The matrix multiplication of A_basis*v is the projection of the node onto the plane formed by the patch. + double guess_arr[2]; + ViewCArrayKokkos guess(&guess_arr[0], 2); + mat_mul(A_basis, v, guess); + xi_val = guess(0); + eta_val = guess(1); + del_tc = del_t/2; + + // Get the solution + bool solution_found = get_contact_point(node, xi_val, eta_val, del_tc); + + if (solution_found && fabs(xi_val) <= 1.0 + edge_tol && fabs(eta_val) <= 1.0 + edge_tol && del_tc >= 0.0 - tol && + del_tc <= del_t + tol) + { + return true; + } else + { + return false; + } +} // end contact_check + +KOKKOS_FUNCTION +void contact_patch_t::ref_to_physical(const ViewCArrayKokkos &ref, const ViewCArrayKokkos &A, + ViewCArrayKokkos &phys) const +{ + const double &xi_ = ref(0); + const double &eta_ = ref(1); + + double phi_k_arr[contact_patch_t::max_nodes]; + ViewCArrayKokkos phi_k(&phi_k_arr[0], num_nodes_in_patch); + phi(phi_k, xi_, eta_); + + mat_mul(A, phi_k, phys); +} + +KOKKOS_FUNCTION +void contact_patch_t::phi(ViewCArrayKokkos &phi_k, const double &xi_value, const double &eta_value) const +{ + if (num_nodes_in_patch == 4) + { + for (int i = 0; i < 4; i++) + { + phi_k(i) = 0.25*(1.0 + xi(i)*xi_value)*(1.0 + eta(i)*eta_value); + } + } else + { + std::cerr << "Error: higher order elements are not yet tested for contact" << std::endl; + exit(1); + } +} // end phi + +#pragma clang diagnostic push +#pragma ide diagnostic ignored "UnusedParameter" +KOKKOS_FUNCTION +void contact_patch_t::d_phi_d_xi(ViewCArrayKokkos &d_phi_k_d_xi, const double &xi_value, + const double &eta_value) const +{ + // xi_value is used for higher order elements + if (num_nodes_in_patch == 4) + { + for (int i = 0; i < 4; i++) + { + d_phi_k_d_xi(i) = 0.25*xi(i)*(1.0 + eta(i)*eta_value); + } + } else + { + std::cerr << "Error: higher order elements are not yet tested for contact" << std::endl; + exit(1); + } +} // end d_phi_d_xi + +KOKKOS_FUNCTION +void contact_patch_t::d_phi_d_eta(ViewCArrayKokkos &d_phi_k_d_eta, const double &xi_value, + const double &eta_value) const +{ + // eta_value is used for higher order elements + if (num_nodes_in_patch == 4) + { + for (int i = 0; i < 4; i++) + { + d_phi_k_d_eta(i) = 0.25*(1.0 + xi(i)*xi_value)*eta(i); + } + } else + { + std::cerr << "Error: higher order elements are not yet tested for contact" << std::endl; + exit(1); + } +} // end d_phi_d_eta + +KOKKOS_FUNCTION +void contact_patch_t::get_normal(const double &xi_val, const double &eta_val, const double &del_t, + ViewCArrayKokkos &normal) const +{ + // The normal is defined as the cross product between dr/dxi and dr/deta where r is the position vector, that is + // r = A*phi_k. + + // Get the derivative arrays + double d_phi_d_xi_arr[contact_patch_t::max_nodes]; + ViewCArrayKokkos d_phi_d_xi_(&d_phi_d_xi_arr[0], num_nodes_in_patch); + d_phi_d_xi(d_phi_d_xi_, xi_val, eta_val); + + double d_phi_d_eta_arr[contact_patch_t::max_nodes]; + ViewCArrayKokkos d_phi_d_eta_(&d_phi_d_eta_arr[0], num_nodes_in_patch); + d_phi_d_eta(d_phi_d_eta_, xi_val, eta_val); + + // Construct the basis matrix A + double A_arr[3*contact_patch_t::max_nodes]; + ViewCArrayKokkos A(&A_arr[0], 3, num_nodes_in_patch); + construct_basis(A, del_t); + + // Get dr_dxi and dr_deta by performing the matrix multiplication A*d_phi_d_xi and A*d_phi_d_eta + double dr_dxi_arr[3]; + ViewCArrayKokkos dr_dxi(&dr_dxi_arr[0], 3); + mat_mul(A, d_phi_d_xi_, dr_dxi); + + double dr_deta_arr[3]; + ViewCArrayKokkos dr_deta(&dr_deta_arr[0], 3); + mat_mul(A, d_phi_d_eta_, dr_deta); + + // Get the normal by performing the cross product between dr_dxi and dr_deta + normal(0) = dr_dxi(1)*dr_deta(2) - dr_dxi(2)*dr_deta(1); + normal(1) = dr_dxi(2)*dr_deta(0) - dr_dxi(0)*dr_deta(2); + normal(2) = dr_dxi(0)*dr_deta(1) - dr_dxi(1)*dr_deta(0); + + // Make the normal a unit vector + double norm_val = norm(normal); + for (int i = 0; i < 3; i++) + { + normal(i) /= norm_val; + } +} // end get_normal + +#pragma clang diagnostic pop +/// end of contact_patch_t member functions //////////////////////////////////////////////////////////////////////////// + +/// beginning of contact_pair_t member functions /////////////////////////////////////////////////////////////////////// +contact_pair_t::contact_pair_t() +{ + // Default constructor +} + +KOKKOS_FUNCTION +contact_pair_t::contact_pair_t(contact_patches_t &contact_patches_obj, const contact_patch_t &patch_obj, + const contact_node_t &node_obj, const double &xi_val, const double &eta_val, + const double &del_tc_val, const ViewCArrayKokkos &normal_view) +{ + // Set the contact_pair_t members + // todo: patch and node might need to be pointers instead so that changes to the patch and node objects are + // reflected in contact_patches_t::contact_nodes and contact_patches_t::contact_patches; however, I think that + // the patch and node objects references in the call mean that this is already happening. + patch = patch_obj; + node = node_obj; + xi = xi_val; + eta = eta_val; + del_tc = del_tc_val; + normal(0) = normal_view(0); + normal(1) = normal_view(1); + normal(2) = normal_view(2); + + contact_patches_obj.is_pen_node(node.gid) = true; + for (int i = 0; i < contact_patch_t::num_nodes_in_patch; i++) + { + contact_patches_obj.is_patch_node(patch.nodes_gid(i)) = true; + } + + // Add the pair to the contact_pairs_access + size_t &patch_stride = contact_patches_obj.contact_pairs_access.stride(patch.lid); + patch_stride++; + contact_patches_obj.contact_pairs_access(patch.lid, patch_stride - 1) = node.gid; +} + +KOKKOS_FUNCTION +void contact_pair_t::frictionless_increment(const double &del_t) +{ + // In order to understand this, just see this PDF: + // https://github.com/gabemorris12/contact_surfaces/blob/master/Finding%20the%20Contact%20Force.pdf + + double A_arr[3*contact_patch_t::max_nodes]; + ViewCArrayKokkos A(&A_arr[0], 3, contact_patch_t::num_nodes_in_patch); + + double phi_k_arr[contact_patch_t::max_nodes]; + ViewCArrayKokkos phi_k(&phi_k_arr[0], contact_patch_t::num_nodes_in_patch); + + double d_phi_d_xi_arr[contact_patch_t::max_nodes]; + ViewCArrayKokkos d_phi_d_xi_(&d_phi_d_xi_arr[0], contact_patch_t::num_nodes_in_patch); + + double d_phi_d_eta_arr[contact_patch_t::max_nodes]; + ViewCArrayKokkos d_phi_d_eta_(&d_phi_d_eta_arr[0], contact_patch_t::num_nodes_in_patch); + + double ak; // place to store acceleration of a patch node + double as; // place to store acceleration of the penetrating node + + double rhs_arr[3]; // right hand side (A_arr*phi_k) + ViewCArrayKokkos rhs(&rhs_arr[0], 3); + + double lhs; // left hand side (node.pos + node.vel*del_t + 0.5*node.acc*del_t*del_t) + + double F_arr[3]; // defined as lhs - rhs + ViewCArrayKokkos F(&F_arr[0], 3); + + double d_A_d_xi_arr[3*contact_patch_t::max_nodes]; // derivative of A wrt xi + ViewCArrayKokkos d_A_d_xi(&d_A_d_xi_arr[0], 3, contact_patch_t::num_nodes_in_patch); + + double d_A_d_eta_arr[3*contact_patch_t::max_nodes]; // derivative of A wrt eta + ViewCArrayKokkos d_A_d_eta(&d_A_d_eta_arr[0], 3, contact_patch_t::num_nodes_in_patch); + + double d_A_d_fc_arr[3*contact_patch_t::max_nodes]; // derivative of A wrt fc_inc + ViewCArrayKokkos d_A_d_fc(&d_A_d_fc_arr[0], 3, contact_patch_t::num_nodes_in_patch); + + double neg_normal_arr[3]; + ViewCArrayKokkos neg_normal(&neg_normal_arr[0], 3); + for (int i = 0; i < 3; i++) + { + neg_normal(i) = -normal(i); + } + + double outer1_arr[contact_patch_t::max_nodes]; // right segment in first outer product + ViewCArrayKokkos outer1(&outer1_arr[0], contact_patch_t::num_nodes_in_patch); + + double outer2_arr[contact_patch_t::max_nodes]; // right segment in second outer product + ViewCArrayKokkos outer2(&outer2_arr[0], contact_patch_t::num_nodes_in_patch); + + double outer3_arr[contact_patch_t::max_nodes]; // right segment in third outer product + ViewCArrayKokkos outer3(&outer3_arr[0], contact_patch_t::num_nodes_in_patch); + + double J0_first_arr[3]; // first term in the J0 column calculation + ViewCArrayKokkos J0_first(&J0_first_arr[0], 3); + + double J0_second_arr[3]; // second term in the J0 column calculation + ViewCArrayKokkos J0_second(&J0_second_arr[0], 3); + + double J1_first_arr[3]; // first term in the J1 column calculation + ViewCArrayKokkos J1_first(&J1_first_arr[0], 3); + + double J1_second_arr[3]; // second term in the J1 column calculation + ViewCArrayKokkos J1_second(&J1_second_arr[0], 3); + + double J2_second_arr[3]; // second term in the J2 column calculation + ViewCArrayKokkos J2_second(&J2_second_arr[0], 3); + + double J0_arr[3]; // column 1 of jacobian + ViewCArrayKokkos J0(&J0_arr[0], 3); + + double J1_arr[3]; // column 2 of jacobian + ViewCArrayKokkos J1(&J1_arr[0], 3); + + double J2_arr[3]; // column 3 of jacobian + ViewCArrayKokkos J2(&J2_arr[0], 3); + + double J_arr[9]; // jacobian + ViewCArrayKokkos J(&J_arr[0], 3, 3); + + double J_inv_arr[9]; // inverse of jacobian + ViewCArrayKokkos J_inv(&J_inv_arr[0], 3, 3); + + double J_det; // determinant of jacobian + double sol[3]; + sol[0] = xi; + sol[1] = eta; + sol[2] = fc_inc; + + double grad_arr[3]; // J_inv*F term + ViewCArrayKokkos grad(&grad_arr[0], 3); + + for (int i = 0; i < max_iter; i++) + { + patch.phi(phi_k, xi, eta); + // construct A + for (int j = 0; j < 3; j++) + { + for (int k = 0; k < contact_patch_t::num_nodes_in_patch; k++) + { + const contact_node_t &patch_node = patch.nodes_obj(k); + ak = (patch_node.internal_force(j) - fc_inc*normal(j)*phi_k(k) + + patch_node.contact_force(j))/patch_node.mass; + A(j, k) = patch_node.pos(j) + patch_node.vel(j)*del_t + 0.5*ak*del_t*del_t; + } + } + + // construct F + mat_mul(A, phi_k, rhs); + for (int j = 0; j < 3; j++) + { + as = (node.internal_force(j) + fc_inc*normal(j) + node.contact_force(j))/node.mass; + lhs = node.pos(j) + node.vel(j)*del_t + 0.5*as*del_t*del_t; + F(j) = lhs - rhs(j); + } + + if (norm(F) <= tol) + { + break; + } + + // construct J + patch.d_phi_d_xi(d_phi_d_xi_, xi, eta); + patch.d_phi_d_eta(d_phi_d_eta_, xi, eta); + + for (int j = 0; j < contact_patch_t::num_nodes_in_patch; j++) + { + const contact_node_t &patch_node = patch.nodes_obj(j); + outer1(j) = (0.5*d_phi_d_xi_(j)*fc_inc*del_t*del_t)/patch_node.mass; + outer2(j) = (0.5*d_phi_d_eta_(j)*fc_inc*del_t*del_t)/patch_node.mass; + outer3(j) = (0.5*phi_k(j)*del_t*del_t)/patch_node.mass; + } + + outer(neg_normal, outer1, d_A_d_xi); + outer(neg_normal, outer2, d_A_d_eta); + outer(neg_normal, outer3, d_A_d_fc); + + mat_mul(A, d_phi_d_xi_, J0_first); + mat_mul(d_A_d_xi, phi_k, J0_second); + for (int j = 0; j < 3; j++) + { + J0(j) = -J0_first(j) - J0_second(j); + } + + mat_mul(A, d_phi_d_eta_, J1_first); + mat_mul(d_A_d_eta, phi_k, J1_second); + for (int j = 0; j < 3; j++) + { + J1(j) = -J1_first(j) - J1_second(j); + } + + mat_mul(d_A_d_fc, phi_k, J2_second); + for (int j = 0; j < 3; j++) + { + J2(j) = (0.5*del_t*del_t*normal(j))/node.mass - J2_second(j); + } + + for (int j = 0; j < 3; j++) + { + J(j, 0) = J0(j); + J(j, 1) = J1(j); + J(j, 2) = J2(j); + } + + J_det = det(J); // there should be no singularities in this calculation + if (J_det == 0.0) + { + fc_inc = 0.0; + printf("Error: Singularity detected in frictionless_increment\n"); + break; + } + + inv(J, J_inv, J_det); + mat_mul(J_inv, F, grad); + for (int j = 0; j < 3; j++) + { + sol[j] = sol[j] - grad(j); + } + xi = sol[0]; + eta = sol[1]; + fc_inc = sol[2]; + } +} + +KOKKOS_FUNCTION +void contact_pair_t::distribute_frictionless_force(const double &force_scale) +{ + double force_val = force_scale*fc_inc; + + // get phi_k + double phi_k_arr[contact_patch_t::max_nodes]; + ViewCArrayKokkos phi_k(&phi_k_arr[0], contact_patch_t::num_nodes_in_patch); + patch.phi(phi_k, xi, eta); + + // if tensile, then subtract left over fc_inc_total; if not, then distribute to nodes + if (force_val + fc_inc_total < 0.0) + { + // update penetrating node + for (int i = 0; i < 3; i++) + { + node.contact_force(i) -= fc_inc_total*normal(i); + } + + // update patch nodes + for (int k = 0; k < contact_patch_t::num_nodes_in_patch; k++) + { + contact_node_t &patch_node = patch.nodes_obj(k); + for (int i = 0; i < 3; i++) + { + patch_node.contact_force(i) += fc_inc_total*normal(i)*phi_k(k); + } + } + + fc_inc_total = 0.0; + fc_inc = 0.0; + } else + { + fc_inc_total += force_val; + + // update penetrating node + for (int i = 0; i < 3; i++) + { + node.contact_force(i) += force_val*normal(i); + } + + // update patch nodes + for (int k = 0; k < contact_patch_t::num_nodes_in_patch; k++) + { + contact_node_t &patch_node = patch.nodes_obj(k); + for (int i = 0; i < 3; i++) + { + patch_node.contact_force(i) -= force_val*normal(i)*phi_k(k); + } + } + } +} // end distribute_frictionless_force + +bool contact_pair_t::should_remove(const double &del_t) +{ + if (fc_inc_total == 0.0 || fabs(xi) > 1.0 + edge_tol || fabs(eta) > 1.0 + edge_tol) + { + fc_inc_total = 0.0; + return true; + } else + { + // update the normal and zero fc_inc_total for the next iteration + double new_normal_arr[3]; + ViewCArrayKokkos new_normal(&new_normal_arr[0], 3); + patch.get_normal(xi, eta, del_t, new_normal); + for (int i = 0; i < 3; i++) + { + normal(i) = new_normal(i); + } + fc_inc_total = 0.0; + + return false; + } +} // end should_remove +/// end of contact_pair_t member functions ///////////////////////////////////////////////////////////////////////////// + +/// beginning of contact_patches_t member functions //////////////////////////////////////////////////////////////////// +void contact_patches_t::initialize(const mesh_t &mesh, const CArrayKokkos &bdy_contact_patches, + const node_t &nodes, const corner_t &corners) +{ + // Contact is only supported in 3D + if (mesh.num_dims != 3) + { + std::cerr << "Error: contact is only supported in 3D" << std::endl; + exit(1); + } + + // Set up the patches that will be checked for contact + patches_gid = bdy_contact_patches; + num_contact_patches = patches_gid.size(); + + // Set up array of contact_patch_t + contact_patches = CArrayKokkos(num_contact_patches); + + CArrayKokkos nodes_in_patch(num_contact_patches, mesh.num_nodes_in_patch); + FOR_ALL_CLASS(i, 0, num_contact_patches, + j, 0, mesh.num_nodes_in_patch, { + nodes_in_patch(i, j) = mesh.nodes_in_patch(patches_gid(i), j); + }); // end parallel for + Kokkos::fence(); + + // todo: Kokkos arrays are being accessed and modified but this isn't inside a macro. Should this be inside Run()? + for (int i = 0; i < num_contact_patches; i++) + { + contact_patch_t &contact_patch = contact_patches(i); + contact_patch.gid = patches_gid(i); + contact_patch.lid = i; + + // Make contact_patch.nodes_gid equal to the row of nodes_in_patch(i) + // This line is what is limiting the parallelism + contact_patch.nodes_gid = CArrayKokkos(mesh.num_nodes_in_patch); + contact_patch.nodes_obj = CArrayKokkos(mesh.num_nodes_in_patch); + for (size_t j = 0; j < mesh.num_nodes_in_patch; j++) + { + contact_patch.nodes_gid(j) = nodes_in_patch(i, j); + } + } // end for + + // todo: This if statement might need a closer look + // Setting up the iso-parametric coordinates for all patch objects + if (mesh.num_nodes_in_patch == 4) + { + contact_patch_t::num_nodes_in_patch = 4; + double xi_temp[4] = {-1.0, 1.0, 1.0, -1.0}; + double eta_temp[4] = {-1.0, -1.0, 1.0, 1.0}; + + // Allocating memory for the points and vel_points arrays + for (int i = 0; i < num_contact_patches; i++) + { + contact_patch_t &contact_patch = contact_patches(i); + + contact_patch.xi = CArrayKokkos(4); + contact_patch.eta = CArrayKokkos(4); + + // todo: We have the Kokkos xi and eta members being modified on host. I think it may be that the xi_temp + // and eta_temp arrays need to be ridden and do a serial Run() for the construction of xi and eta. + for (int j = 0; j < 4; j++) + { + contact_patch.xi(j) = xi_temp[j]; + contact_patch.eta(j) = eta_temp[j]; + } + } + + } else + { + std::cerr << "Error: higher order elements are not yet tested for contact" << std::endl; + exit(1); + } // end if + + // Determine the bucket size. This is defined as 1.001*min_node_distance + CArrayKokkos node_distances(num_contact_patches, contact_patch_t::num_nodes_in_patch); + FOR_ALL_CLASS(i, 0, num_contact_patches, + j, 0, contact_patch_t::num_nodes_in_patch, { + const contact_patch_t &contact_patch = contact_patches(i); + if (j < contact_patch_t::num_nodes_in_patch - 1) + { + const size_t n1 = contact_patch.nodes_gid(j); // current node + const size_t n2 = contact_patch.nodes_gid(j + 1); // next node + + double sum_sq = 0.0; + + for (int k = 0; k < 3; k++) + { + sum_sq += pow(nodes.coords(0, n1, k) - nodes.coords(0, n2, k), 2); + } + + node_distances(i, j) = sqrt(sum_sq); + } else + { + const size_t n1 = contact_patch.nodes_gid(0); // current node + const size_t n2 = contact_patch.nodes_gid( + contact_patch_t::num_nodes_in_patch - 1); // next node + + double sum_sq = 0.0; + + for (int k = 0; k < 3; k++) + { + sum_sq += pow(nodes.coords(0, n1, k) - nodes.coords(0, n2, k), 2); + } + + node_distances(i, j) = sqrt(sum_sq); + } + }); + Kokkos::fence(); + + double result = 0.0; + double local_min = 1.0e10; + REDUCE_MIN(i, 0, num_contact_patches, + j, 0, contact_patch_t::num_nodes_in_patch, local_min, { + if (local_min > node_distances(i, j)) + { + local_min = node_distances(i, j); + } + }, result); + + contact_patches_t::bucket_size = 0.999*result; + + // Find the total number of nodes (this is should always be less than or equal to mesh.num_bdy_nodes) + size_t local_max_index = 0; + size_t max_index = 0; + REDUCE_MAX(i, 0, num_contact_patches, + j, 0, contact_patch_t::num_nodes_in_patch, local_max_index, { + if (local_max_index < contact_patches(i).nodes_gid(j)) + { + local_max_index = contact_patches(i).nodes_gid(j); + } + }, max_index); + Kokkos::fence(); + + CArrayKokkos node_count(max_index + 1); + + // zero node_count + FOR_ALL(i, 0, max_index + 1, { + node_count(i) = 0; + }); + Kokkos::fence(); + + for (int i = 0; i < num_contact_patches; i++) + { + contact_patch_t &contact_patch = contact_patches(i); + for (int j = 0; j < contact_patch_t::num_nodes_in_patch; j++) + { + size_t node_gid = contact_patch.nodes_gid(j); + if (node_count(node_gid) == 0) + { + node_count(node_gid) = 1; + contact_patches_t::num_contact_nodes += 1; + } + } + } + + // todo: instead of these arrays being accessed through the node gid directly, they can be made much smaller with a + // size of num_contact_nodes. The nsort member should simply be a sorted array of local indices (see the + // final portion of sort()). Then, use nodes_gid to get the global index to be used for mesh_t, node_t, and + // corner_t objects. + // Initialize the contact_nodes and contact_pairs arrays + contact_nodes = CArrayKokkos(max_index + 1); + contact_pairs = CArrayKokkos(max_index + 1); + contact_pairs_access = DynamicRaggedRightArrayKokkos(num_contact_patches, + contact_patch_t::max_contacting_nodes_in_patch); + is_patch_node = CArrayKokkos(max_index + 1); + is_pen_node = CArrayKokkos(max_index + 1); + active_pairs = CArrayKokkos(contact_patches_t::num_contact_nodes); + forces = CArrayKokkos(contact_patches_t::num_contact_nodes); + + // Construct nodes_gid and nodes_obj + nodes_gid = CArrayKokkos(contact_patches_t::num_contact_nodes); + size_t node_lid = 0; + for (int i = 0; i < num_contact_patches; i++) + { + contact_patch_t &contact_patch = contact_patches(i); + for (int j = 0; j < contact_patch_t::num_nodes_in_patch; j++) + { + size_t node_gid = contact_patch.nodes_gid(j); + contact_node_t &node_obj = contact_nodes(node_gid); + // todo: for some reason, these non-matar types do not update in the contact_patches_t::update_nodes func + // this is only a problem if the node mass is changing + node_obj.mass = nodes.mass(node_gid); + node_obj.gid = node_gid; + contact_patch.nodes_obj(j) = node_obj; + if (node_count(node_gid) == 1) + { + node_count(node_gid) = 2; + nodes_gid(node_lid) = node_gid; + node_lid += 1; + } + } + } + + // Construct patches_in_node and num_patches_in_node + num_patches_in_node = CArrayKokkos(max_index + 1); + for (int patch_lid = 0; patch_lid < num_contact_patches; patch_lid++) + { + const contact_patch_t &contact_patch = contact_patches(patch_lid); + // for each node in the patch, increment the counter in num_patches_in_node + FOR_ALL_CLASS(i, 0, contact_patch_t::num_nodes_in_patch, { + const size_t &node_gid = contact_patch.nodes_gid(i); + num_patches_in_node(node_gid) += 1; + }); + Kokkos::fence(); + } + + CArrayKokkos stride_index(max_index + 1); + patches_in_node = RaggedRightArrayKokkos(num_patches_in_node); + // Walk through the patches, and for each node in the patch, add the patch to the patches_in_node array + for (int patch_lid = 0; patch_lid < num_contact_patches; patch_lid++) + { + // todo: will the device have access to patch_lid as I have it here? + const contact_patch_t &contact_patch = contact_patches(patch_lid); + FOR_ALL_CLASS(i, 0, contact_patch_t::num_nodes_in_patch, { + const size_t &node_gid = contact_patch.nodes_gid(i); + const size_t &stride = stride_index(node_gid); + patches_in_node(node_gid, stride) = patch_lid; + stride_index(node_gid) += 1; + }); + Kokkos::fence(); + } + + // Update the node members + update_nodes(mesh, nodes, corners); + +} // end initialize + +void contact_patches_t::update_nodes(const mesh_t &mesh, const node_t &nodes, const corner_t &corner) +{ + // Update node objects + FOR_ALL_CLASS(i, 0, contact_patches_t::num_contact_nodes, { + const size_t &node_gid = nodes_gid(i); + contact_node_t &contact_node = contact_nodes(node_gid); + contact_node.gid = node_gid; + contact_node.mass = nodes.mass(node_gid); + + // Update pos, vel, acc, and internal force + for (int j = 0; j < 3; j++) + { + contact_node.pos(j) = nodes.coords(0, node_gid, j); + contact_node.vel(j) = nodes.vel(0, node_gid, j); + + // zero forces + contact_node.contact_force(j) = 0.0; + contact_node.internal_force(j) = 0.0; + + // Loop over the corners + for (size_t corner_lid = 0; corner_lid < mesh.num_corners_in_node(node_gid); corner_lid++) + { + size_t corner_gid = mesh.corners_in_node(node_gid, corner_lid); + contact_node.internal_force(j) += corner.force(corner_gid, j); + } + } + }); + + // todo: if node mass is changing, then this will need to be updated here by looping through contact patches as well + // I'm not sure why, but the non-matar members (gid and mass) of contact_node_t is not updating the + // nodes_obj member in contact_patch_t +} + +void contact_patches_t::sort() +{ + // todo: I don't think it's a good idea to have these allocated here. These should become member variables, and + // its allocation should be moved to initialize() because sort() is being called every step. + // Grouping all the coordinates, velocities, and accelerations + CArrayKokkos points(3, contact_patch_t::num_nodes_in_patch*num_contact_patches); + CArrayKokkos velocities(3, contact_patch_t::num_nodes_in_patch*num_contact_patches); + CArrayKokkos accelerations(3, contact_patch_t::num_nodes_in_patch*num_contact_patches); + + FOR_ALL(i, 0, num_contact_patches, { + const contact_patch_t &contact_patch = contact_patches(i); + for (int j = 0; j < contact_patch_t::num_nodes_in_patch; j++) + { + const contact_node_t &node = contact_patch.nodes_obj(j); + for (int k = 0; k < 3; k++) + { + points(k, i*contact_patch_t::num_nodes_in_patch + j) = node.pos(k); + velocities(k, i*contact_patch_t::num_nodes_in_patch + j) = fabs(node.vel(k)); + accelerations(k, i*contact_patch_t::num_nodes_in_patch + j) = fabs((node.internal_force(k) + + node.contact_force(k))/node.mass); + } + } + }); // end parallel for + Kokkos::fence(); + + // Find the max and min for each dimension + double local_x_max = 0.0; + REDUCE_MAX_CLASS(i, 0, contact_patch_t::num_nodes_in_patch*num_contact_patches, local_x_max, { + if (local_x_max < points(0, i)) + { + local_x_max = points(0, i); + } + }, x_max); + double local_y_max = 0.0; + REDUCE_MAX_CLASS(i, 0, contact_patch_t::num_nodes_in_patch*num_contact_patches, local_y_max, { + if (local_y_max < points(1, i)) + { + local_y_max = points(1, i); + } + }, y_max); + double local_z_max = 0.0; + REDUCE_MAX_CLASS(i, 0, contact_patch_t::num_nodes_in_patch*num_contact_patches, local_z_max, { + if (local_z_max < points(2, i)) + { + local_z_max = points(2, i); + } + }, z_max); + double local_x_min = 0.0; + REDUCE_MIN_CLASS(i, 0, contact_patch_t::num_nodes_in_patch*num_contact_patches, local_x_min, { + if (local_x_min > points(0, i)) + { + local_x_min = points(0, i); + } + }, x_min); + double local_y_min = 0.0; + REDUCE_MIN_CLASS(i, 0, contact_patch_t::num_nodes_in_patch*num_contact_patches, local_y_min, { + if (local_y_min > points(1, i)) + { + local_y_min = points(1, i); + } + }, y_min); + double local_z_min = 0.0; + REDUCE_MIN_CLASS(i, 0, contact_patch_t::num_nodes_in_patch*num_contact_patches, local_z_min, { + if (local_z_min > points(2, i)) + { + local_z_min = points(2, i); + } + }, z_min); + double local_vx_max = 0.0; + REDUCE_MAX_CLASS(i, 0, contact_patch_t::num_nodes_in_patch*num_contact_patches, local_vx_max, { + if (local_vx_max < velocities(0, i)) + { + local_vx_max = velocities(0, i); + } + }, vx_max); + double local_vy_max = 0.0; + REDUCE_MAX_CLASS(i, 0, contact_patch_t::num_nodes_in_patch*num_contact_patches, local_vy_max, { + if (local_vy_max < velocities(1, i)) + { + local_vy_max = velocities(1, i); + } + }, vy_max); + double local_vz_max = 0.0; + REDUCE_MAX_CLASS(i, 0, contact_patch_t::num_nodes_in_patch*num_contact_patches, local_vz_max, { + if (local_vz_max < velocities(2, i)) + { + local_vz_max = velocities(2, i); + } + }, vz_max); + double local_ax_max = 0.0; + REDUCE_MAX_CLASS(i, 0, contact_patch_t::num_nodes_in_patch*num_contact_patches, local_ax_max, { + if (local_ax_max < accelerations(0, i)) + { + local_ax_max = accelerations(0, i); + } + }, ax_max); + double local_ay_max = 0.0; + REDUCE_MAX_CLASS(i, 0, contact_patch_t::num_nodes_in_patch*num_contact_patches, local_ay_max, { + if (local_ay_max < accelerations(1, i)) + { + local_ay_max = accelerations(1, i); + } + }, ay_max); + double local_az_max = 0.0; + REDUCE_MAX_CLASS(i, 0, contact_patch_t::num_nodes_in_patch*num_contact_patches, local_az_max, { + if (local_az_max < accelerations(2, i)) + { + local_az_max = accelerations(2, i); + } + }, az_max); + Kokkos::fence(); + + // If the max velocity is zero, then we want to set it to a small value. The max velocity and acceleration are used + // for creating a capture box around the contact patch. We want there to be at least some thickness to the box. + double* vel_max[3] = {&vx_max, &vy_max, &vz_max}; + for (auto &i: vel_max) + { + if (*i == 0.0) + { + *i = 1.0e-3; // Set to a small value + } + } + + // Define Sx, Sy, and Sz + Sx = floor((x_max - x_min)/bucket_size) + 1; // NOLINT(*-narrowing-conversions) + Sy = floor((y_max - y_min)/bucket_size) + 1; // NOLINT(*-narrowing-conversions) + Sz = floor((z_max - z_min)/bucket_size) + 1; // NOLINT(*-narrowing-conversions) + + // todo: Something similar to the points, velocities, and accelerations arrays above should be done here. The issue + // with this one is that nb changes through the iterations. lbox and npoint might need to change to Views. + // Initializing the nbox, lbox, nsort, and npoint arrays + size_t nb = Sx*Sy*Sz; // total number of buckets + nbox = CArrayKokkos(nb); + lbox = CArrayKokkos(contact_patches_t::num_contact_nodes); + nsort = CArrayKokkos(contact_patches_t::num_contact_nodes); + npoint = CArrayKokkos(nb); + CArrayKokkos nsort_lid(contact_patches_t::num_contact_nodes); + + // Find the bucket id for each node by constructing lbox + FOR_ALL_CLASS(i, 0, contact_patches_t::num_contact_nodes, { + size_t node_gid = nodes_gid(i); + const contact_node_t &node = contact_nodes(node_gid); + double x = node.pos(0); + double y = node.pos(1); + double z = node.pos(2); + + size_t Si_x = floor((x - x_min)/bucket_size); + size_t Si_y = floor((y - y_min)/bucket_size); + size_t Si_z = floor((z - z_min)/bucket_size); + + lbox(i) = Si_z*Sx*Sy + Si_y*Sx + Si_x; + nbox(lbox(i)) += 1; // increment nbox + }); + Kokkos::fence(); + + // Calculate the pointer for each bucket into a sorted list of nodes + for (size_t i = 1; i < nb; i++) + { + npoint(i) = npoint(i - 1) + nbox(i - 1); + } + + // Zero nbox + FOR_ALL_CLASS(i, 0, nb, { + nbox(i) = 0; + }); + Kokkos::fence(); + + // Sort the slave nodes according to their bucket id into nsort + for (int i = 0; i < contact_patches_t::num_contact_nodes; i++) + { + nsort_lid(nbox(lbox(i)) + npoint(lbox(i))) = i; + nbox(lbox(i)) += 1; + } + + // Change nsort to reflect the global node id's + FOR_ALL_CLASS(i, 0, contact_patches_t::num_contact_nodes, { + nsort(i) = nodes_gid(nsort_lid(i)); + }); + Kokkos::fence(); +} // end sort + +void contact_patches_t::find_nodes(contact_patch_t &contact_patch, const double &del_t, + size_t &num_nodes_found) +{ + // Get capture box + contact_patch.capture_box(vx_max, vy_max, vz_max, ax_max, ay_max, az_max, del_t); + const CArrayKokkos &bounds = contact_patch.bounds; + + // Determine the buckets that intersect with the capture box + size_t ibox_max = fmax(0, fmin(Sx - 1, floor((bounds(0) - x_min)/bucket_size))); // NOLINT(*-narrowing-conversions) + size_t jbox_max = fmax(0, fmin(Sy - 1, floor((bounds(1) - y_min)/bucket_size))); // NOLINT(*-narrowing-conversions) + size_t kbox_max = fmax(0, fmin(Sz - 1, floor((bounds(2) - z_min)/bucket_size))); // NOLINT(*-narrowing-conversions) + size_t ibox_min = fmax(0, fmin(Sx - 1, floor((bounds(3) - x_min)/bucket_size))); // NOLINT(*-narrowing-conversions) + size_t jbox_min = fmax(0, fmin(Sy - 1, floor((bounds(4) - y_min)/bucket_size))); // NOLINT(*-narrowing-conversions) + size_t kbox_min = fmax(0, fmin(Sz - 1, floor((bounds(5) - z_min)/bucket_size))); // NOLINT(*-narrowing-conversions) + + size_t bucket_index = 0; + for (size_t i = ibox_min; i < ibox_max + 1; i++) + { + for (size_t j = jbox_min; j < jbox_max + 1; j++) + { + for (size_t k = kbox_min; k < kbox_max + 1; k++) + { + // todo: If there is a segfault here, then that means that we have more than max_number_buckets. + // increase that value as this means that there was so much strain that the number of buckets + // exceeded the value. + contact_patch.buckets(bucket_index) = k*Sx*Sy + j*Sx + i; + bucket_index += 1; + } + } + } + + // Get all nodes in each bucket + num_nodes_found = 0; // local node index for possible_nodes member + for (int bucket_lid = 0; bucket_lid < bucket_index; bucket_lid++) + { + size_t b = contact_patch.buckets(bucket_lid); + for (size_t i = 0; i < nbox(b); i++) + { + size_t node_gid = nsort(npoint(b) + i); + bool add_node = true; + // If the node is in the current contact_patch, then continue; else, add it to possible_nodes + for (int j = 0; j < contact_patch_t::num_nodes_in_patch; j++) + { + if (node_gid == contact_patch.nodes_gid(j)) + { + add_node = false; + break; + } + } + + if (add_node) + { + contact_patch.possible_nodes(num_nodes_found) = node_gid; + num_nodes_found += 1; + } + } + } +} // end find_nodes + +void contact_patches_t::get_contact_pairs(const double &del_t) +{ + // clear the is_patch_node and is_pen_node arrays to be false + FOR_ALL_CLASS(i, 0, is_patch_node.size(), { + is_patch_node(i) = false; + is_pen_node(i) = false; + }); + Kokkos::fence(); + + for (int patch_lid = 0; patch_lid < num_contact_patches; patch_lid++) + { + contact_patch_t &contact_patch = contact_patches(patch_lid); + size_t num_nodes_found; + find_nodes(contact_patch, del_t, num_nodes_found); + // todo: original plan was for this to be a parallel loop, but there are collisions in the contact_pairs_access + for (int node_lid = 0; node_lid < num_nodes_found; node_lid++) + { + const size_t &node_gid = contact_patch.possible_nodes(node_lid); + const contact_node_t &node = contact_nodes(node_gid); + contact_pair_t ¤t_pair = contact_pairs(node_gid); + // If this is already an active pair, then back out of the thread/continue + if (current_pair.active) + { + is_pen_node(node_gid) = true; + for (int i = 0; i < contact_patch_t::num_nodes_in_patch; i++) + { + is_patch_node(contact_patch.nodes_gid(i)) = true; + } + // todo: ending the thread here has performance consequences for cuda; though it may be fine to do this + // because the number of possible nodes is small + // return; + continue; + } + + double xi_val, eta_val, del_tc; + bool is_hitting = contact_patch.contact_check(node, del_t, xi_val, eta_val, del_tc); + + if (is_hitting && !is_pen_node(node_gid) && !is_patch_node(node_gid)) + { + double normal_arr[3]; + ViewCArrayKokkos normal(&normal_arr[0], 3); + contact_patch.get_normal(xi_val, eta_val, del_t, normal); + current_pair = contact_pair_t(*this, contact_patch, node, xi_val, eta_val, del_tc, normal); + } else if (is_hitting && is_pen_node(node_gid)) + { + // This means that the current_pair has already been initialized. We need to compare the pair stored in + // `current_pair` to the parameters in this iteration. If `current_pair.del_tc` is larger than `del_tc` + // from this iteration, then that means that this iteration will intersect the patch first. For this, + // we need to update the parameters in `current_pair` to reflect the ones here. The only thing that + // stays constant is the node, while the patch, xi, eta, del_tc, etc. are updated. + if (del_tc + tol < current_pair.del_tc) + { + // see edge_cases.py for testing this branch (first test in file) + // remove all the patch nodes from is_patch_node + const contact_patch_t &original_patch = current_pair.patch; + for (int i = 0; i < contact_patch_t::num_nodes_in_patch; i++) + { + is_patch_node(original_patch.nodes_gid(i)) = false; + } + + // todo: remove_pair might be a reason why this whole loop should be serial + remove_pair(current_pair); + double normal_arr[3]; + ViewCArrayKokkos normal(&normal_arr[0], 3); + contact_patch.get_normal(xi_val, eta_val, del_t, normal); + current_pair = contact_pair_t(*this, contact_patch, node, xi_val, eta_val, del_tc, normal); + } else if (current_pair.del_tc - tol <= del_tc && current_pair.del_tc + tol >= del_tc) + { + // see edge_cases.py second and fourth test in file + + // This means that the node is hitting an edge. The dominant pair to be selected is based off the + // normal at the penetrating node's surface and the normal at the contact point. + double normal1_arr[3]; + ViewCArrayKokkos normal1(&normal1_arr[0], 3); + // current pair stores normal1 + for (int i = 0; i < 3; i++) + { + normal1(i) = current_pair.normal(i); + } + + double normal2_arr[3]; + ViewCArrayKokkos normal2(&normal2_arr[0], 3); + contact_patch.get_normal(xi_val, eta_val, del_t, normal2); + double new_normal_arr[3]; + ViewCArrayKokkos new_normal(&new_normal_arr[0], 3); + bool add_new_pair = get_edge_pair(normal1, normal2, node_gid, del_t, new_normal); + + if (add_new_pair) + { + // remove all the patch nodes from is_patch_node + const contact_patch_t &original_patch = current_pair.patch; + for (int i = 0; i < contact_patch_t::num_nodes_in_patch; i++) + { + is_patch_node(original_patch.nodes_gid(i)) = false; + } + + remove_pair(current_pair); + current_pair = contact_pair_t(*this, contact_patch, node, xi_val, eta_val, del_tc, new_normal); + } + } + } else if (is_hitting && is_patch_node(node_gid)) + { + // This means that the current node is a patch node in a previous contact pair but is also being + // considered as a penetrating node in this iteration. The logic with this case is loop through the + // nodes in the patch of this iteration and see if the node_gid exists as a penetrating node in a + // contact pair. If it does, then a comparison must be made with all the found pairs. + bool add_current_pair = false; + bool hitting_before_arr[contact_patch_t::max_nodes]; // stores true or false if the node of this iteration is hitting before its adjacent pairs + ViewCArrayKokkos hitting_before(&hitting_before_arr[0], contact_patch_t::num_nodes_in_patch); + size_t hitting_index = 0; + + for (int i = 0; i < contact_patch_t::num_nodes_in_patch; i++) + { + const size_t &patch_node_gid = contact_patch.nodes_gid(i); + if (is_pen_node(patch_node_gid)) + { + const contact_pair_t &pair = contact_pairs(patch_node_gid); // adjacent pair + if (del_tc + tol < pair.del_tc) + { + // this means that the current node is hitting the patch object before the patch object + // node is hitting its patch + hitting_before(hitting_index) = true; + // todo: If the contact pair is a glue condition, then we need to remove the adjacent pair + // here. It might be best to remove it for any case. + } else if (pair.del_tc - tol <= del_tc && pair.del_tc + tol >= del_tc) + { + // this means we are hitting at the same time + // if the contact point is not on the edge, then we hitting before is true + if (fabs(xi_val) < 1.0 - tol && fabs(eta_val) < 1.0 - tol) + { + hitting_before(hitting_index) = true; + } else + { + hitting_before(hitting_index) = false; + } + } else + { + // this means that the adjacent pair is hitting before the current node + hitting_before(hitting_index) = false; + } + hitting_index += 1; + } + } + + if (hitting_index == 0) + { + add_current_pair = true; + } else if (all(hitting_before, hitting_index)) + { + add_current_pair = true; + } else if (any(hitting_before, hitting_index) && fabs(xi_val) < 1.0 - tol && + fabs(eta_val) < 1.0 - tol) + { + add_current_pair = true; + } + + if (add_current_pair) + { + double normal_arr[3]; + ViewCArrayKokkos normal(&normal_arr[0], 3); + contact_patch.get_normal(xi_val, eta_val, del_t, normal); + current_pair = contact_pair_t(*this, contact_patch, node, xi_val, eta_val, del_tc, normal); + } + } + } + } + + // set the active pairs + num_active_pairs = 0; + for (int patch_lid = 0; patch_lid < num_contact_patches; patch_lid++) + { + for (int patch_stride = 0; patch_stride < contact_pairs_access.stride(patch_lid); patch_stride++) + { + const size_t &node_gid = contact_pairs_access(patch_lid, patch_stride); + contact_pair_t &pair = contact_pairs(node_gid); + pair.active = true; + active_pairs(num_active_pairs) = node_gid; + num_active_pairs += 1; + } + } +} // end get_contact_pairs + +KOKKOS_FUNCTION +void contact_patches_t::remove_pair(contact_pair_t &pair) +{ + pair.active = false; + pair.fc_inc_total = 0.0; + pair.fc_inc = 0.0; + + // modify the contact_pairs_access array + // keep iterating in the row until the pair is found and shift all the elements to the left + // then decrement the stride + size_t &patch_stride = contact_pairs_access.stride(pair.patch.lid); + bool found_node = false; + for (size_t i = 0; i < patch_stride; i++) + { + const size_t &node_gid = contact_pairs_access(pair.patch.lid, i); + if (node_gid == pair.node.gid) + { + found_node = true; + } else if (found_node) + { + contact_pairs_access(pair.patch.lid, i - 1) = node_gid; + } + } + patch_stride -= 1; + assert(found_node && "Error: attempted to remove pair that doesn't exist in contact_pairs_access"); +} // end remove_pair + +KOKKOS_FUNCTION +bool contact_patches_t::get_edge_pair(const ViewCArrayKokkos &normal1, const ViewCArrayKokkos &normal2, + const size_t &node_gid, const double &del_t, + ViewCArrayKokkos &new_normal) const +{ + // Get the surface normal of the penetrating node by averaging all the normals at that node + // we do this by looping through all the patches that the node is in + double node_normal_arr[3]; + ViewCArrayKokkos node_normal(&node_normal_arr[0], 3); + // zero node normal + for (int i = 0; i < 3; i++) + { + node_normal(i) = 0.0; + } + + const size_t &num_patches = num_patches_in_node(node_gid); + double local_normal_arr[3]; + ViewCArrayKokkos local_normal(&local_normal_arr[0], 3); + for (size_t i = 0; i < num_patches; i++) + { + const contact_patch_t &patch = contact_patches(patches_in_node(node_gid, i)); + // loop through the nodes in the patch until we find the node_gid the index matches with patch.xi and patch.eta + for (int j = 0; j < contact_patch_t::num_nodes_in_patch; j++) + { + if (patch.nodes_gid(j) == node_gid) + { + // get the normal at that node + patch.get_normal(patch.xi(j), patch.eta(j), del_t, local_normal); + for (int k = 0; k < 3; k++) + { + node_normal(k) += local_normal(k); + } + break; + } + } + } + // finish getting the average + for (int i = 0; i < 3; i++) + { + node_normal(i) /= num_patches; + } + // Make it a unit vector + double normal_norm = norm(node_normal); + for (int i = 0; i < 3; i++) + { + node_normal(i) /= normal_norm; + } + + // returning true means that a new pair should be formed + // the pair that should be selected is the one that has the most negative dot product with the node normal + // if normal1 is the most negative, then return false and make new_normal = normal1 + // if normal2 is the most negative, then return true and make new_normal = normal2 + // if the normals are the same, then return true and make new_normal the average between the two + double dot1 = dot(normal1, node_normal); + double dot2 = dot(normal2, node_normal); + + if (dot2 - tol <= dot1 && dot2 + tol >= dot1) + { + for (int i = 0; i < 3; i++) + { + new_normal(i) = (normal1(i) + normal2(i))/2.0; + } + + // make new_normal a unit vector + double new_normal_norm = norm(new_normal); + for (int i = 0; i < 3; i++) + { + new_normal(i) /= new_normal_norm; + } + return true; + } else if (dot1 < dot2) + { + for (int i = 0; i < 3; i++) + { + new_normal(i) = normal1(i); + } + return false; + } else + { + for (int i = 0; i < 3; i++) + { + new_normal(i) = normal2(i); + } + return true; + } +} // end get_edge_pair + +void contact_patches_t::force_resolution(const double &del_t) +{ + ViewCArrayKokkos forces_view(&forces(0), num_active_pairs); + for (int i = 0; i < max_iter; i++) + { + // find force increment for each pair + // todo: having trouble breaking this into two parts, find force and apply force. It's not giving the right + // results so keep it serial for now. + for (int j = 0; j < num_active_pairs; j++) + { + const size_t &node_gid = active_pairs(j); + contact_pair_t &pair = contact_pairs(node_gid); + + // The reason why we are doing if statements inside the loop instead of loops inside of if statements is to + // be able to have different contact types per pair. This makes it to where we can have a portion of the + // mesh do frictionless contact and another portion do glue contact. + if (pair.contact_type == contact_pair_t::contact_types::frictionless) + { + pair.frictionless_increment(del_t); + pair.distribute_frictionless_force(1.0); // if not doing serial, then this would be called in the second loop + forces_view(j) = pair.fc_inc; + } // else if (pair.contact_type == contact_pair_t::contact_types::glue) + } + // Kokkos::fence(); + + // check convergence (the force increments should be zero) + if (norm(forces_view) <= tol) + { + break; + } + + // distribute forces to nodes + // todo: this loop is a little more complicated since we are accessing patches at the same time + // keep this serial until determining a way to use kokkos atomics + // for (int j = 0; j < num_active_pairs; j++) + // { + // const size_t &node_gid = active_pairs(j); + // contact_pair_t &pair = contact_pairs(node_gid); + + // if (pair.contact_type == contact_pair_t::contact_types::frictionless) + // { + // pair.distribute_frictionless_force(*this, 0.5); + // } // else if (pair.contact_type == contact_pair_t::contact_types::glue) + // } + + // todo: multi-stage detection should get added here. To do this, wrap this whole function in a while loop and + // implement the multi-stage detection as is in the python version. + } +} // end force_resolution + +void contact_patches_t::remove_pairs(const double &del_t) +{ + for (int i = 0; i < num_active_pairs; i++) + { + const size_t &node_gid = active_pairs(i); + contact_pair_t &pair = contact_pairs(node_gid); + + bool should_remove = false; + if (pair.contact_type == contact_pair_t::contact_types::frictionless) + { + should_remove = pair.should_remove(del_t); + } // else if (pair.contact_type == contact_pair_t::contact_types::glue) + + if (should_remove) + { + remove_pair(pair); + } + } +} +/// end of contact_patches_t member functions ////////////////////////////////////////////////////////////////////////// + +/// beginning of internal, not to be used anywhere else tests ////////////////////////////////////////////////////////// +contact_patch_t::contact_patch_t() = default; + +contact_patch_t::contact_patch_t(const ViewCArrayKokkos &points, const ViewCArrayKokkos &vel_points, + const ViewCArrayKokkos &internal_force_points, + const ViewCArrayKokkos &contact_force_points, + const ViewCArrayKokkos &mass_points_) +{ + this->xi = CArrayKokkos(4); + this->eta = CArrayKokkos(4); + xi(0) = -1.0; + xi(1) = 1.0; + xi(2) = 1.0; + xi(3) = -1.0; + eta(0) = -1.0; + eta(1) = -1.0; + eta(2) = 1.0; + eta(3) = 1.0; + + this->nodes_obj = CArrayKokkos(4); + + for (int i = 0; i < num_nodes_in_patch; i++) + { + ViewCArrayKokkos pos_view = ViewCArrayKokkos(&points(i, 0), 3); + ViewCArrayKokkos vel_view = ViewCArrayKokkos(&vel_points(i, 0), 3); + ViewCArrayKokkos internal_force_view = ViewCArrayKokkos(&internal_force_points(i, 0), 3); + ViewCArrayKokkos contact_force_view = ViewCArrayKokkos(&contact_force_points(i, 0), 3); + + contact_node_t node = contact_node_t(pos_view, vel_view, internal_force_view, contact_force_view, + mass_points_(i)); + nodes_obj(i) = node; + } +} + +contact_node_t::contact_node_t() = default; + +contact_node_t::contact_node_t(const ViewCArrayKokkos &pos, const ViewCArrayKokkos &vel, + const ViewCArrayKokkos &internal_force, + const ViewCArrayKokkos &contact_force, const double &mass) +{ + this->mass = mass; + for (int i = 0; i < 3; i++) + { + this->pos(i) = pos(i); + this->vel(i) = vel(i); + this->internal_force(i) = internal_force(i); + this->contact_force(i) = contact_force(i); + } +} + +void run_contact_tests(contact_patches_t &contact_patches_obj, const mesh_t &mesh, const node_t &nodes, + const corner_t &corner, const simulation_parameters_t &sim_params) +{ + double err_tol = 1.0e-6; // error tolerance + + // Testing get_contact_point with abnormal patch velocities. See contact_check_visual_through_reference.py. + double test1_points_arr[4*3] = {1.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 1.0, + 1.0, 0.0, 1.0}; + double test1_vels_arr[4*3] = {0.0, 1.0, 0.0, + 0.0, 0.1, 0.0, + 0.0, 0.2, 0.0, + 0.0, 0.0, 0.0}; + double test1_internal_force_arr[4*3] = {0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0}; + double test1_contact_force_arr[4*3] = {0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0}; + double test1_mass_points_arr[4] = {1.0, 1.0, 1.0, 1.0}; + ViewCArrayKokkos test1_points(&test1_points_arr[0], 4, 3); + ViewCArrayKokkos test1_vels(&test1_vels_arr[0], 4, 3); + ViewCArrayKokkos test1_internal_force(&test1_internal_force_arr[0], 4, 3); + ViewCArrayKokkos test1_contact_force(&test1_contact_force_arr[0], 4, 3); + ViewCArrayKokkos test1_mass_points(&test1_mass_points_arr[0], 4); + contact_patch_t test1_patch(test1_points, test1_vels, test1_internal_force, test1_contact_force, + test1_mass_points); + + double test1_node_pos[3] = {0.25, 1.0, 0.2}; + double test1_node_vel[3] = {0.75, -1.0, 0.0}; + double test1_node_internal_force_arr[3] = {0.0, 0.0, 0.0}; + double test1_node_contact_force_arr[3] = {0.0, 0.0, 0.0}; + ViewCArrayKokkos test1_pos(&test1_node_pos[0], 3); + ViewCArrayKokkos test1_vel(&test1_node_vel[0], 3); + ViewCArrayKokkos test1_node_internal_force(&test1_node_internal_force_arr[0], 3); + ViewCArrayKokkos test1_node_contact_force(&test1_node_contact_force_arr[0], 3); + contact_node_t test1_node(test1_pos, test1_vel, test1_node_internal_force, test1_node_contact_force, + 1.0); + + double xi_val, eta_val, del_tc; + bool is_hitting = test1_patch.get_contact_point(test1_node, xi_val, eta_val, del_tc); + bool contact_check = test1_patch.contact_check(test1_node, 1.0, xi_val, eta_val, del_tc); + std::cout << "\nTesting get_contact_point and contact_check:" << std::endl; + std::cout << "-0.433241 -0.6 0.622161 vs. "; + std::cout << xi_val << " " << eta_val << " " << del_tc << std::endl; + assert(fabs(xi_val + 0.43324096) < err_tol); + assert(fabs(eta_val + 0.6) < err_tol); + assert(fabs(del_tc - 0.6221606424928471) < err_tol); + assert(is_hitting); + assert(contact_check); + + // Testing contact force calculation with the previous pair. See contact_check_visual_through_reference.py. + std::cout << "\nTesting frictionless_increment:" << std::endl; + + contact_pair_t test1_pair; + test1_pair.patch = test1_patch; + test1_pair.node = test1_node; + test1_pair.xi = xi_val; + test1_pair.eta = eta_val; + test1_pair.del_tc = del_tc; + + double force_normal[3]; + ViewCArrayKokkos force_n(&force_normal[0], 3); + // Normal is not taken at del_tc, but is taken at the current time step; this is just for testing purposes + test1_pair.patch.get_normal(test1_pair.xi, test1_pair.eta, test1_pair.del_tc, force_n); + for (int i = 0; i < 3; i++) + { + test1_pair.normal(i) = force_n(i); + } + + test1_pair.fc_inc = 0.5; + test1_pair.frictionless_increment(1.0); + std::cout << "-0.581465 -0.176368 0.858551 vs. "; + std::cout << test1_pair.xi << " " << test1_pair.eta << " " << test1_pair.fc_inc << std::endl; + assert(fabs(test1_pair.xi + 0.581465) < err_tol); + assert(fabs(test1_pair.eta + 0.176368) < err_tol); + assert(fabs(test1_pair.fc_inc - 0.858551) < err_tol); + + // Testing sort and get_contact_pairs + std::string file_name = sim_params.mesh_input.file_path; + std::string main_test = "contact_test.geo"; + std::string edge_case1 = "edge_case1.geo"; + std::string edge_case2 = "edge_case2.geo"; + std::string edge_case3 = "edge_case3.geo"; + + std::cout << "\nTesting sort and get_contact_pairs:" << std::endl; + if (file_name.find(main_test) != std::string::npos) + { + std::cout << "Patch with nodes 10 11 5 4 is paired with node 22" << std::endl; + std::cout << "Patch with nodes 9 10 4 3 is paired with node 23" << std::endl; + std::cout << "Patch with nodes 16 17 11 10 is paired with node 18" << std::endl; + std::cout << "Patch with nodes 15 16 10 9 is paired with node 19" << std::endl; + std::cout << "Patch with nodes 18 19 23 22 is paired with node 10" << std::endl; + std::cout << "vs." << std::endl; + contact_patches_obj.sort(); + contact_patches_obj.get_contact_pairs(0.1); + for (int i = 0; i < contact_patches_obj.num_contact_patches; i++) + { + for (int j = 0; j < contact_patches_obj.contact_pairs_access.stride(i); j++) + { + size_t node_gid = contact_patches_obj.contact_pairs_access(i, j); + contact_pair_t &pair = contact_patches_obj.contact_pairs(node_gid); + std::cout << "Patch with nodes "; + for (int k = 0; k < contact_patch_t::num_nodes_in_patch; k++) + { + std::cout << pair.patch.nodes_gid(k) << " "; + } + std::cout << "is paired with node " << pair.node.gid << std::endl; + } + } + assert(contact_patches_obj.contact_pairs_access(2, 0) == 22); + assert(contact_patches_obj.contact_pairs_access(6, 0) == 23); + assert(contact_patches_obj.contact_pairs_access(10, 0) == 18); + assert(contact_patches_obj.contact_pairs_access(14, 0) == 19); + assert(contact_patches_obj.contact_pairs_access(18, 0) == 10); + } else if (file_name.find(edge_case1) != std::string::npos) + { + std::cout << "Patch with nodes 7 8 2 1 is paired with node 12" << std::endl; + std::cout << "Patch with nodes 7 8 2 1 is paired with node 16" << std::endl; + std::cout << "vs." << std::endl; + contact_patches_obj.sort(); + contact_patches_obj.get_contact_pairs(1.0); + for (int i = 0; i < contact_patches_obj.num_contact_patches; i++) + { + for (int j = 0; j < contact_patches_obj.contact_pairs_access.stride(i); j++) + { + size_t node_gid = contact_patches_obj.contact_pairs_access(i, j); + contact_pair_t &pair = contact_patches_obj.contact_pairs(node_gid); + std::cout << "Patch with nodes "; + for (int k = 0; k < contact_patch_t::num_nodes_in_patch; k++) + { + std::cout << pair.patch.nodes_gid(k) << " "; + } + std::cout << "is paired with node " << pair.node.gid << std::endl; + } + } + assert(contact_patches_obj.contact_pairs_access(6, 0) == 12); + assert(contact_patches_obj.contact_pairs_access(6, 1) == 16); + } else if (file_name.find(edge_case2) != std::string::npos) + { + std::cout << "Patch with nodes 7 8 2 1 is paired with node 12" << std::endl; + std::cout << "Patch with nodes 7 8 2 1 is paired with node 13" << std::endl; + std::cout << "Patch with nodes 7 8 2 1 is paired with node 16" << std::endl; + std::cout << "Patch with nodes 7 8 2 1 is paired with node 17" << std::endl; + std::cout << "vs." << std::endl; + contact_patches_obj.sort(); + contact_patches_obj.get_contact_pairs(1.0); + for (int i = 0; i < contact_patches_obj.num_contact_patches; i++) + { + for (int j = 0; j < contact_patches_obj.contact_pairs_access.stride(i); j++) + { + size_t node_gid = contact_patches_obj.contact_pairs_access(i, j); + contact_pair_t &pair = contact_patches_obj.contact_pairs(node_gid); + std::cout << "Patch with nodes "; + for (int k = 0; k < contact_patch_t::num_nodes_in_patch; k++) + { + std::cout << pair.patch.nodes_gid(k) << " "; + } + std::cout << "is paired with node " << pair.node.gid << std::endl; + } + } + assert(contact_patches_obj.contact_pairs_access(6, 0) == 12); + assert(contact_patches_obj.contact_pairs_access(6, 1) == 13); + assert(contact_patches_obj.contact_pairs_access(6, 2) == 16); + assert(contact_patches_obj.contact_pairs_access(6, 3) == 17); + } else if (file_name.find(edge_case3) != std::string::npos) + { + std::cout << "Patch with nodes 6 7 1 0 is paired with node 12"; + std::cout << " ---> Pushback Direction: 0 -0.447214 0.894427" << std::endl; + std::cout << "Patch with nodes 6 7 1 0 is paired with node 18"; + std::cout << " ---> Pushback Direction: 0 -0.447214 0.894427" << std::endl; + std::cout << "Patch with nodes 7 8 2 1 is paired with node 13"; + std::cout << " ---> Pushback Direction: 0 0 1" << std::endl; + std::cout << "Patch with nodes 7 8 2 1 is paired with node 14"; + std::cout << " ---> Pushback Direction: 0 0.447214 0.894427" << std::endl; + std::cout << "Patch with nodes 7 8 2 1 is paired with node 19"; + std::cout << " ---> Pushback Direction: 0 0 1" << std::endl; + std::cout << "Patch with nodes 7 8 2 1 is paired with node 20"; + std::cout << " ---> Pushback Direction: 0 0.447214 0.894427" << std::endl; + std::cout << "vs." << std::endl; + contact_patches_obj.sort(); + contact_patches_obj.get_contact_pairs(1.0); + for (int i = 0; i < contact_patches_obj.num_contact_patches; i++) + { + for (int j = 0; j < contact_patches_obj.contact_pairs_access.stride(i); j++) + { + size_t node_gid = contact_patches_obj.contact_pairs_access(i, j); + contact_pair_t &pair = contact_patches_obj.contact_pairs(node_gid); + std::cout << "Patch with nodes "; + for (int k = 0; k < contact_patch_t::num_nodes_in_patch; k++) + { + std::cout << pair.patch.nodes_gid(k) << " "; + } + std::cout << "is paired with node " << pair.node.gid; + std::cout << " ---> Pushback Direction: "; + for (int k = 0; k < 3; k++) + { + std::cout << pair.normal(k) << " "; + } + std::cout << std::endl; + } + } + assert(contact_patches_obj.contact_pairs_access(1, 0) == 12); + assert(contact_patches_obj.contact_pairs_access(1, 1) == 18); + assert(contact_patches_obj.contact_pairs_access(6, 0) == 13); + assert(contact_patches_obj.contact_pairs_access(6, 1) == 14); + assert(contact_patches_obj.contact_pairs_access(6, 2) == 19); + assert(contact_patches_obj.contact_pairs_access(6, 3) == 20); + } + + // Testing force resolution + std::cout << "\nTesting force_resolution and remove_pairs:" << std::endl; + if (file_name.find(main_test) != std::string::npos) + { + std::cout << "Penetrating node 22 has a fc_inc_total value of 0.148976" << std::endl; + std::cout << "Penetrating node 23 has a fc_inc_total value of 0.148976" << std::endl; + std::cout << "Penetrating node 18 has a fc_inc_total value of 0.148976" << std::endl; + std::cout << "Penetrating node 19 has a fc_inc_total value of 0.148976" << std::endl; + std::cout << "Penetrating node 10 has a fc_inc_total value of 0" << std::endl; + std::cout << "vs." << std::endl; + contact_patches_obj.force_resolution(0.1); + + double pen_node_sum = 0.0; + double patch_node_sum = 0.0; + bool seen_patch_node[26]; + for (int i = 0; i < 26; i++) + { + seen_patch_node[i] = false; + } + for (int i = 0; i < contact_patches_obj.num_active_pairs; i++) + { + const size_t &node_gid = contact_patches_obj.active_pairs(i); + const contact_pair_t &pair = contact_patches_obj.contact_pairs(node_gid); + std::cout << "Penetrating node " << node_gid << " has a fc_inc_total_value of "; + std::cout << pair.fc_inc_total << std::endl; + + pen_node_sum += pair.fc_inc_total; + for (int j = 0; j < contact_patch_t::num_nodes_in_patch; j++) + { + const size_t &patch_node_gid = pair.patch.nodes_gid(j); + if (!seen_patch_node[patch_node_gid] && pair.fc_inc_total > 0.0) + { + seen_patch_node[patch_node_gid] = true; + + const contact_node_t &patch_node = contact_patches_obj.contact_nodes(patch_node_gid); + patch_node_sum += patch_node.contact_force(2); + } + } + } + std::cout << "Penetrating node sum: " << pen_node_sum << std::endl; + std::cout << "Patch node sum: " << patch_node_sum << std::endl; + + assert(fabs(pen_node_sum + patch_node_sum) < err_tol); + assert(fabs(contact_patches_obj.contact_pairs(22).fc_inc_total - 0.148976) < err_tol); + assert(fabs(contact_patches_obj.contact_pairs(23).fc_inc_total - 0.148976) < err_tol); + assert(fabs(contact_patches_obj.contact_pairs(18).fc_inc_total - 0.148976) < err_tol); + assert(fabs(contact_patches_obj.contact_pairs(19).fc_inc_total - 0.148976) < err_tol); + assert(contact_patches_obj.contact_pairs(10).fc_inc_total < err_tol); + + // Testing remove_pairs + contact_patches_obj.remove_pairs(0.1); + assert(!contact_patches_obj.contact_pairs(10).active); + } + + exit(0); +} +/// end of internal, not to be used anywhere else tests //////////////////////////////////////////////////////////////// diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/force_sgh.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/force_sgh.cpp index 2cdb89e9d..85b98312f 100644 --- a/single-node-refactor/src/Solvers/SGH_solver/src/force_sgh.cpp +++ b/single-node-refactor/src/Solvers/SGH_solver/src/force_sgh.cpp @@ -59,7 +59,7 @@ /// \param The current Runge Kutta integration alpha value /// ///////////////////////////////////////////////////////////////////////////// -void SGH::get_force(const CArrayKokkos& material, +void SGH::get_force(const DCArrayKokkos& material, const mesh_t& mesh, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, @@ -366,32 +366,32 @@ void SGH::get_force(const CArrayKokkos& material, } // end for loop over nodes in elem // --- Update Stress --- - // calculate the new stress at the next rk level, if it is a hypo model + // calculate the new stress at the next rk level, if it is a increment_based size_t mat_id = elem_mat_id(elem_gid); - // hypo elastic plastic model - if (material(mat_id).strength_type == model::hypo) { + // increment_based strength model + if (material(mat_id).strength_type == model::increment_based) { // cut out the node_gids for this element ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), 8); // --- call strength model --- - // material(mat_id).strength_model(elem_pres, - // elem_stress, - // elem_gid, - // mat_id, - // elem_statev, - // elem_sspd, - // elem_den(elem_gid), - // elem_sie(elem_gid), - // vel_grad, - // elem_node_gids, - // node_coords, - // node_vel, - // elem_vol(elem_gid), - // dt, - // rk_alpha); - } // end logical on hypo strength model + material(mat_id).strength_model(elem_pres, + elem_stress, + elem_gid, + mat_id, + elem_statev, + elem_sspd, + elem_den(elem_gid), + elem_sie(elem_gid), + vel_grad, + elem_node_gids, + node_coords, + node_vel, + elem_vol(elem_gid), + dt, + rk_alpha); + } // end logical on increment_based strength model }); // end parallel for loop over elements return; @@ -422,7 +422,7 @@ void SGH::get_force(const CArrayKokkos& material, /// \param The current Runge Kutta integration alpha value /// ///////////////////////////////////////////////////////////////////////////// -void SGH::get_force_2D(const CArrayKokkos& material, +void SGH::get_force_2D(const DCArrayKokkos& material, const mesh_t& mesh, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, @@ -740,32 +740,32 @@ void SGH::get_force_2D(const CArrayKokkos& material, } // end for loop over nodes in elem // --- Update Stress --- - // calculate the new stress at the next rk level, if it is a hypo model + // calculate the new stress at the next rk level, if it is a Increment based model size_t mat_id = elem_mat_id(elem_gid); - // hypo elastic plastic model - if (material(mat_id).strength_type == model::hypo) { + // Increment based elastic plastic model + if (material(mat_id).strength_type == model::increment_based) { // cut out the node_gids for this element - ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), 4); + ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), 8); // --- call strength model --- - // material(mat_id).strength_model(elem_pres, - // elem_stress, - // elem_gid, - // mat_id, - // elem_statev, - // elem_sspd, - // elem_den(elem_gid), - // elem_sie(elem_gid), - // vel_grad, - // elem_node_gids, - // node_coords, - // node_vel, - // elem_vol(elem_gid), - // dt, - // rk_alpha); - } // end logical on hypo strength model + material(mat_id).strength_model(elem_pres, + elem_stress, + elem_gid, + mat_id, + elem_statev, + elem_sspd, + elem_den(elem_gid), + elem_sie(elem_gid), + vel_grad, + elem_node_gids, + node_coords, + node_vel, + elem_vol(elem_gid), + dt, + rk_alpha); + } // end logical on increment_based strength model }); // end parallel for loop over elements return; diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/momentum.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/momentum.cpp index 7eddbd3fe..b69ad968a 100644 --- a/single-node-refactor/src/Solvers/SGH_solver/src/momentum.cpp +++ b/single-node-refactor/src/Solvers/SGH_solver/src/momentum.cpp @@ -1,5 +1,5 @@ /********************************************************************************************** - © 2020. Triad National Security, LLC. All rights reserved. + � 2020. Triad National Security, LLC. All rights reserved. This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. Department of Energy/National Nuclear Security Administration. All rights in the program are @@ -52,18 +52,30 @@ void SGH::update_velocity(double rk_alpha, const mesh_t& mesh, DCArrayKokkos& node_vel, const DCArrayKokkos& node_mass, - const DCArrayKokkos& corner_force + const DCArrayKokkos& corner_force, + const CArrayKokkos &contact_nodes ) { const size_t num_dims = mesh.num_dims; // walk over the nodes to update the velocity - FOR_ALL(node_gid, 0, mesh.num_nodes, { + //FOR_ALL(node_gid, 0, mesh.num_nodes, { + for(int node_gid = 0; node_gid& material, +void SGH::update_state(const DCArrayKokkos& material, const mesh_t& mesh, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, @@ -90,7 +90,7 @@ void SGH::update_state(const CArrayKokkos& material, // --- Stress --- // hyper elastic plastic model - if (material(mat_id).strength_type == model::hyper) { + if (material(mat_id).strength_type == model::state_based) { // cut out the node_gids for this element ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), num_nodes_in_elem); @@ -117,22 +117,22 @@ void SGH::update_state(const CArrayKokkos& material, elem_gid); // --- call strength model --- - // material(mat_id).strength_model(elem_pres, - // elem_stress, - // elem_gid, - // mat_id, - // elem_statev, - // elem_sspd, - // elem_den(elem_gid), - // elem_sie(elem_gid), - // vel_grad, - // elem_node_gids, - // node_coords, - // node_vel, - // elem_vol(elem_gid), - // dt, - // rk_alpha); - } // end logical on hyper strength model + material(mat_id).strength_model(elem_pres, + elem_stress, + elem_gid, + mat_id, + elem_statev, + elem_sspd, + elem_den(elem_gid), + elem_sie(elem_gid), + vel_grad, + elem_node_gids, + node_coords, + node_vel, + elem_vol(elem_gid), + dt, + rk_alpha); + } // end logical on state_based strength model // --- Pressure --- material(mat_id).eos_model(elem_pres, @@ -172,7 +172,7 @@ void SGH::update_state(const CArrayKokkos& material, /// \param The current Runge Kutta integration alpha value /// ///////////////////////////////////////////////////////////////////////////// -void SGH::update_state2D(const CArrayKokkos& material, +void SGH::update_state2D(const DCArrayKokkos& material, const mesh_t& mesh, const DCArrayKokkos& node_coords, const DCArrayKokkos& node_vel, @@ -203,8 +203,8 @@ void SGH::update_state2D(const CArrayKokkos& material, size_t mat_id = elem_mat_id(elem_gid); // --- Stress --- - // hyper elastic plastic model - if (material(mat_id).strength_type == model::hyper) { + // state_based elastic plastic model + if (material(mat_id).strength_type == model::state_based) { // cut out the node_gids for this element ViewCArrayKokkos elem_node_gids(&mesh.nodes_in_elem(elem_gid, 0), num_nodes_in_elem); @@ -246,7 +246,7 @@ void SGH::update_state2D(const CArrayKokkos& material, // elem_vol(elem_gid), // dt, // rk_alpha); - } // end logical on hyper strength model + } // end logical on state_based strength model // --- Pressure --- // material(mat_id).eos_model(elem_pres, diff --git a/single-node-refactor/src/Solvers/SGH_solver/src/sgh_solve.cpp b/single-node-refactor/src/Solvers/SGH_solver/src/sgh_solve.cpp index a725ee39c..e9b688c39 100644 --- a/single-node-refactor/src/Solvers/SGH_solver/src/sgh_solve.cpp +++ b/single-node-refactor/src/Solvers/SGH_solver/src/sgh_solve.cpp @@ -1,5 +1,5 @@ /********************************************************************************************** - © 2020. Triad National Security, LLC. All rights reserved. + � 2020. Triad National Security, LLC. All rights reserved. This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. Department of Energy/National Nuclear Security Administration. All rights in the program are @@ -49,8 +49,8 @@ void SGH::execute(simulation_parameters_t& sim_param, mesh_t& mesh, node_t& node { std::cout << "In execute function in sgh solver" << std::endl; - // printf("Writing outputs to file at %f \n", time_value); - // mesh_writer.write_mesh(mesh, elem, node, corner, sim_param, time_value, graphics_times); + printf("Writing outputs to file at %f \n", time_value); + mesh_writer.write_mesh(mesh, elem, node, corner, sim_param, time_value, graphics_times); CArrayKokkos node_extensive_mass(mesh.num_nodes); @@ -232,23 +232,42 @@ void SGH::execute(simulation_parameters_t& sim_param, mesh_t& mesh, node_t& node rk_alpha); } + // ---- apply contact boundary conditions to the boundary patches---- + if (doing_contact) // Structuring it like this to avoid having to sort() everytime. + { + contact_bank.update_nodes(mesh, node, corner); + boundary_contact(dt*rk_alpha); + } + + // Check to see if contact force all sums to zero + // double total_force[3] = {0.0, 0.0, 0.0}; + // for (size_t i = 0; i < mesh.num_nodes; i++) + // { + // const contact_node_t &contact_node = contact_bank.contact_nodes(i); + // for (size_t j = 0; j < 3; j++) + // { + // total_force[j] += contact_node.contact_force(j); + // } + // } + + // ViewCArrayKokkos total_force_view(&total_force[0], 3); + // print out total_force + // printf("Total force: %f %f %f\n", total_force[0], total_force[1], total_force[2]); + + // mpi_coms(); + // ---- Update nodal velocities ---- // update_velocity(rk_alpha, dt, mesh, node.vel, node.mass, - corner.force); + corner.force, + contact_bank.contact_nodes); // ---- apply velocity boundary conditions to the boundary patches---- boundary_velocity(mesh, sim_param.boundary_conditions, node.vel, time_value); - - // ---- apply contact boundary conditions to the boundary patches---- - boundary_contact(mesh, sim_param.boundary_conditions, node.vel, time_value); - - // mpi_coms(); - // ---- Update specific internal energy in the elements ---- update_energy(rk_alpha, dt, diff --git a/single-node-refactor/src/common/boundary_conditions.h b/single-node-refactor/src/common/boundary_conditions.h index 9abd6f2b4..dd7e8ec60 100644 --- a/single-node-refactor/src/common/boundary_conditions.h +++ b/single-node-refactor/src/common/boundary_conditions.h @@ -1,5 +1,5 @@ /********************************************************************************************** - © 2020. Triad National Security, LLC. All rights reserved. + � 2020. Triad National Security, LLC. All rights reserved. This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. Department of Energy/National Nuclear Security Administration. All rights in the program are @@ -47,7 +47,8 @@ enum bdy_tag y_plane = 1, // tag an y-plane z_plane = 2, // tag an z-plane cylinder = 3, // tag an cylindrical surface - sphere = 4 // tag a spherical surface + sphere = 4, // tag a spherical surface + global = 5 // tag all boundary patches (used for contact only) //read_file = 5 // read from a file currently unsupported }; @@ -80,7 +81,8 @@ static std::map bc_geometry_map { "y_plane", boundary_conds::y_plane }, { "z_plane", boundary_conds::z_plane }, { "cylinder", boundary_conds::cylinder }, - { "sphere", boundary_conds::sphere } + { "sphere", boundary_conds::sphere }, + { "global", boundary_conds::global } // { "read_file", boundary_conds::read_file } }; @@ -112,7 +114,7 @@ static std::map bc_direction_map ///////////////////////////////////////////////////////////////////////////// struct boundary_condition_t { - solver_input::method solver = solver_input::NONE; ///< Numerical solver method + solver_input::method solver = solver_input::NONE; ///< Numerical solver method boundary_conds::bdy_hydro_conds type; ///< Type of boundary condition boundary_conds::bdy_tag geometry; ///< Geometry boundary condition is applied to @@ -123,7 +125,7 @@ struct boundary_condition_t double v = 0.0; ///< WARNING: Currently unused double w = 0.0; ///< WARNING: Currently unused - std::vector origin = { 0.0, 0.0, 0.0 }; ///< Origin of boundary condition geometry + DCArrayKokkos origin; // // WARNING: CURRENTLY NOT PARSED OR USED double hydro_bc_vel_0 = 0.0; ///< Initial velocity for timed velocity boundary condition @@ -131,7 +133,6 @@ struct boundary_condition_t double hydro_bc_vel_t_start = 0.0; ///< Start time for velocity boundary condition double hydro_bc_vel_t_end = 0.0; ///< End time for velocity boundary condition }; // end boundary conditions - // ------------------------------------- // valid inputs for boundary conditions // ------------------------------------- diff --git a/single-node-refactor/src/common/geometry_new.h b/single-node-refactor/src/common/geometry_new.h index 15856920f..39b15f22f 100644 --- a/single-node-refactor/src/common/geometry_new.h +++ b/single-node-refactor/src/common/geometry_new.h @@ -708,7 +708,7 @@ size_t check_bdy(const size_t patch_gid, /// ///////////////////////////////////////////////////////////////////////////// KOKKOS_INLINE_FUNCTION -void tag_bdys(const CArrayKokkos& boundary, +void tag_bdys(const DCArrayKokkos& boundary, mesh_t& mesh, const DCArrayKokkos& node_coords) { diff --git a/single-node-refactor/src/common/io_utils.h b/single-node-refactor/src/common/io_utils.h index 9b343728c..f151a515a 100644 --- a/single-node-refactor/src/common/io_utils.h +++ b/single-node-refactor/src/common/io_utils.h @@ -1,36 +1,36 @@ /********************************************************************************************** - © 2020. Triad National Security, LLC. All rights reserved. - This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos - National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. - Department of Energy/National Nuclear Security Administration. All rights in the program are - reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear - Security Administration. The Government is granted for itself and others acting on its behalf a - nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare - derivative works, distribute copies to the public, perform publicly and display publicly, and - to permit others to do so. - This program is open source under the BSD-3 License. - Redistribution and use in source and binary forms, with or without modification, are permitted - provided that the following conditions are met: - 1. Redistributions of source code must retain the above copyright notice, this list of - conditions and the following disclaimer. - 2. Redistributions in binary form must reproduce the above copyright notice, this list of - conditions and the following disclaimer in the documentation and/or other materials - provided with the distribution. - 3. Neither the name of the copyright holder nor the names of its contributors may be used - to endorse or promote products derived from this software without specific prior - written permission. - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS - IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR - CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; - OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, - WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR - OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF - ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - **********************************************************************************************/ +© 2020. Triad National Security, LLC. All rights reserved. +This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos +National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. +Department of Energy/National Nuclear Security Administration. All rights in the program are +reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear +Security Administration. The Government is granted for itself and others acting on its behalf a +nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare +derivative works, distribute copies to the public, perform publicly and display publicly, and +to permit others to do so. +This program is open source under the BSD-3 License. +Redistribution and use in source and binary forms, with or without modification, are permitted +provided that the following conditions are met: +1. Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, this list of +conditions and the following disclaimer in the documentation and/or other materials +provided with the distribution. +3. Neither the name of the copyright holder nor the names of its contributors may be used +to endorse or promote products derived from this software without specific prior +written permission. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; +OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +**********************************************************************************************/ #ifndef FIERRO_IO_H #define FIERRO_IO_H @@ -50,8 +50,8 @@ /// /// \brief Class for simplifying reading meshes /// -/// This class contains the requisite functions required to read different -/// mesh formats. The idea is to set the mesh file name, and parse the +/// This class contains the requisite functions required to read different +/// mesh formats. The idea is to set the mesh file name, and parse the /// extension to decide which reader to use. Currently, only ensight .geo /// files are supported. /// @@ -175,8 +175,6 @@ class MeshReader else{ double dummy; fscanf(in, "%le", &dummy); - - } } // end for @@ -287,7 +285,6 @@ class MeshBuilder std::cout << "\t" << pair.first << std::endl; } throw std::runtime_error("**** 2D MESH TYPE NOT SUPPORTED ****"); - return; } } else if (sim_param.mesh_input.num_dims == 3) { @@ -311,17 +308,17 @@ class MeshBuilder /// \param Simulation parameters /// ///////////////////////////////////////////////////////////////////////////// - void build_2d_box(mesh_t& mesh, elem_t& elem, node_t& node, corner_t& corner, simulation_parameters_t& sim_param) + void build_2d_box(mesh_t& mesh, elem_t& elem, node_t& node, corner_t& corner, simulation_parameters_t& sim_param) const { printf(" Creating a 2D box mesh \n"); const int num_dim = 2; - const double lx = sim_param.mesh_input.length[0]; - const double ly = sim_param.mesh_input.length[1]; + const double lx = sim_param.mesh_input.length.host(0); + const double ly = sim_param.mesh_input.length.host(1); - const int num_elems_i = sim_param.mesh_input.num_elems[0]; - const int num_elems_j = sim_param.mesh_input.num_elems[1]; + const int num_elems_i = sim_param.mesh_input.num_elems.host(0); + const int num_elems_j = sim_param.mesh_input.num_elems.host(1); const int num_points_i = num_elems_i + 1; // num points in x const int num_points_j = num_elems_j + 1; // num points in y @@ -333,13 +330,15 @@ class MeshBuilder const int num_elems = num_elems_i * num_elems_j; - std::vector origin = sim_param.mesh_input.origin; + std::vector origin(num_dim); + sim_param.mesh_input.origin.update_host(); + for (int i = 0; i < num_dim; i++) { origin[i] = sim_param.mesh_input.origin.host(i); } // --- 2D parameters --- - const int num_faces_in_elem = 4; // number of faces in elem - const int num_points_in_elem = 4; // number of points in elem - const int num_points_in_face = 2; // number of points in a face - const int num_edges_in_elem = 4; // number of edges in a elem + // const int num_faces_in_elem = 4; // number of faces in elem + // const int num_points_in_elem = 4; // number of points in elem + // const int num_points_in_face = 2; // number of points in a face + // const int num_edges_in_elem = 4; // number of edges in a elem // --- mesh node ordering --- // Convert ijk index system to the finite element numbering convention @@ -365,17 +364,18 @@ class MeshBuilder int node_gid = get_id(i, j, 0, num_points_i, num_points_j); // store the point coordinates - node.coords(0, node_gid, 0) = origin[0] + (double)i * dx; - node.coords(0, node_gid, 1) = origin[1] + (double)j * dy; + node.coords.host(0, node_gid, 0) = origin[0] + (double)i * dx; + node.coords.host(0, node_gid, 1) = origin[1] + (double)j * dy; } // end for i } // end for j for (int rk_level = 1; rk_level < rk_num_bins; rk_level++) { for (int node_gid = 0; node_gid < num_nodes; node_gid++) { - node.coords(rk_level, node_gid, 0) = node.coords(0, node_gid, 0); - node.coords(rk_level, node_gid, 1) = node.coords(0, node_gid, 1); + node.coords.host(rk_level, node_gid, 0) = node.coords.host(0, node_gid, 0); + node.coords.host(rk_level, node_gid, 1) = node.coords.host(0, node_gid, 1); } } + node.coords.update_device(); // intialize elem variables mesh.initialize_elems(num_elems, num_dim); @@ -435,7 +435,7 @@ class MeshBuilder /// \param Simulation parameters /// ///////////////////////////////////////////////////////////////////////////// - void build_2d_polar(mesh_t& mesh, elem_t& elem, node_t& node, corner_t& corner, simulation_parameters_t& sim_param) + void build_2d_polar(mesh_t& mesh, elem_t& elem, node_t& node, corner_t& corner, simulation_parameters_t& sim_param) const { printf(" Creating a 2D polar mesh \n"); @@ -461,13 +461,15 @@ class MeshBuilder const int num_elems = num_elems_i * num_elems_j; - std::vector origin = sim_param.mesh_input.origin; + std::vector origin(num_dim); + sim_param.mesh_input.origin.update_host(); + for (int i = 0; i < num_dim; i++) { origin[i] = sim_param.mesh_input.origin.host(i); } // --- 2D parameters --- - const int num_faces_in_elem = 4; // number of faces in elem - const int num_points_in_elem = 4; // number of points in elem - const int num_points_in_face = 2; // number of points in a face - const int num_edges_in_elem = 4; // number of edges in a elem + // const int num_faces_in_elem = 4; // number of faces in elem + // const int num_points_in_elem = 4; // number of points in elem + // const int num_points_in_face = 2; // number of points in a face + // const int num_edges_in_elem = 4; // number of edges in a elem // --- mesh node ordering --- // Convert ijk index system to the finite element numbering convention @@ -492,17 +494,18 @@ class MeshBuilder double theta_j = start_angle + (double)j * dy; // store the point coordinates - node.coords(0, node_gid, 0) = origin[0] + r_i * cos(theta_j); - node.coords(0, node_gid, 1) = origin[1] + r_i * sin(theta_j); + node.coords.host(0, node_gid, 0) = origin[0] + r_i * cos(theta_j); + node.coords.host(0, node_gid, 1) = origin[1] + r_i * sin(theta_j); } // end for i } // end for j for (int rk_level = 1; rk_level < rk_num_bins; rk_level++) { for (int node_gid = 0; node_gid < num_nodes; node_gid++) { - node.coords(rk_level, node_gid, 0) = node.coords(0, node_gid, 0); - node.coords(rk_level, node_gid, 1) = node.coords(0, node_gid, 1); + node.coords.host(rk_level, node_gid, 0) = node.coords.host(0, node_gid, 0); + node.coords.host(rk_level, node_gid, 1) = node.coords.host(0, node_gid, 1); } } + node.coords.update_device(); // intialize elem variables mesh.initialize_elems(num_elems, num_dim); @@ -562,19 +565,21 @@ class MeshBuilder /// \param Simulation parameters /// ///////////////////////////////////////////////////////////////////////////// - void build_3d_box(mesh_t& mesh, elem_t& elem, node_t& node, corner_t& corner, simulation_parameters_t& sim_param) + void build_3d_box(mesh_t& mesh, elem_t& elem, node_t& node, corner_t& corner, simulation_parameters_t& sim_param) const { printf(" Creating a 3D box mesh \n"); const int num_dim = 3; - const double lx = sim_param.mesh_input.length[0]; - const double ly = sim_param.mesh_input.length[1]; - const double lz = sim_param.mesh_input.length[2]; + sim_param.mesh_input.length.update_host(); + const double lx = sim_param.mesh_input.length.host(0); + const double ly = sim_param.mesh_input.length.host(1); + const double lz = sim_param.mesh_input.length.host(2); - const int num_elems_i = sim_param.mesh_input.num_elems[0]; - const int num_elems_j = sim_param.mesh_input.num_elems[1]; - const int num_elems_k = sim_param.mesh_input.num_elems[2]; + sim_param.mesh_input.num_elems.update_host(); + const int num_elems_i = sim_param.mesh_input.num_elems.host(0); + const int num_elems_j = sim_param.mesh_input.num_elems.host(1); + const int num_elems_k = sim_param.mesh_input.num_elems.host(2); const int num_points_i = num_elems_i + 1; // num points in x const int num_points_j = num_elems_j + 1; // num points in y @@ -588,13 +593,15 @@ class MeshBuilder const int num_elems = num_elems_i * num_elems_j * num_elems_k; - std::vector origin = sim_param.mesh_input.origin; + std::vector origin(num_dim); + sim_param.mesh_input.origin.update_host(); + for (int i = 0; i < num_dim; i++) { origin[i] = sim_param.mesh_input.origin.host(i); } // --- 3D parameters --- - const int num_faces_in_elem = 6; // number of faces in elem - const int num_points_in_elem = 8; // number of points in elem - const int num_points_in_face = 4; // number of points in a face - const int num_edges_in_elem = 12; // number of edges in a elem + // const int num_faces_in_elem = 6; // number of faces in elem + // const int num_points_in_elem = 8; // number of points in elem + // const int num_points_in_face = 4; // number of points in a face + // const int num_edges_in_elem = 12; // number of edges in a elem // --- mesh node ordering --- // Convert ijk index system to the finite element numbering convention @@ -625,20 +632,21 @@ class MeshBuilder int node_gid = get_id(i, j, k, num_points_i, num_points_j); // store the point coordinates - node.coords(0, node_gid, 0) = origin[0] + (double)i * dx; - node.coords(0, node_gid, 1) = origin[1] + (double)j * dy; - node.coords(0, node_gid, 2) = origin[2] + (double)k * dz; + node.coords.host(0, node_gid, 0) = origin[0] + (double)i * dx; + node.coords.host(0, node_gid, 1) = origin[1] + (double)j * dy; + node.coords.host(0, node_gid, 2) = origin[2] + (double)k * dz; } // end for i } // end for j } // end for k for (int rk_level = 1; rk_level < rk_num_bins; rk_level++) { for (int node_gid = 0; node_gid < num_nodes; node_gid++) { - node.coords(rk_level, node_gid, 0) = node.coords(0, node_gid, 0); - node.coords(rk_level, node_gid, 1) = node.coords(0, node_gid, 1); - node.coords(rk_level, node_gid, 2) = node.coords(0, node_gid, 2); + node.coords.host(rk_level, node_gid, 0) = node.coords.host(0, node_gid, 0); + node.coords.host(rk_level, node_gid, 1) = node.coords.host(0, node_gid, 1); + node.coords.host(rk_level, node_gid, 2) = node.coords.host(0, node_gid, 2); } } + node.coords.update_device(); // intialize elem variables mesh.initialize_elems(num_elems, num_dim); @@ -695,7 +703,7 @@ class MeshBuilder /// /// \fn build_3d_HexN_box /// - /// \brief Builds an unstructured high order 3D rectilinear mesh + /// \brief Builds an unstructured high order 3D rectilinear mesh /// /// \param Simulation mesh that is built /// \param Element state data @@ -704,7 +712,7 @@ class MeshBuilder /// \param Simulation parameters /// ///////////////////////////////////////////////////////////////////////////// - void build_3d_HexN_box(mesh_t& mesh, elem_t& elem, node_t& node, corner_t& corner, simulation_parameters_t& sim_param) + void build_3d_HexN_box(mesh_t& mesh, elem_t& elem, node_t& node, corner_t& corner, simulation_parameters_t& sim_param) const { printf(" ***** WARNING:: build_3d_HexN_box not yet implemented\n"); } @@ -726,7 +734,7 @@ class MeshBuilder /// \param Number of j indices /// ///////////////////////////////////////////////////////////////////////////// - int get_id(int i, int j, int k, int num_i, int num_j) + int get_id(int i, int j, int k, int num_i, int num_j) const { return i + j * num_i + k * num_i * num_j; } @@ -735,7 +743,7 @@ class MeshBuilder /// /// \fn PointIndexFromIJK /// - /// \brief Given (i,j,k) coordinates within the Lagrange hex, return an + /// \brief Given (i,j,k) coordinates within the Lagrange hex, return an /// offset into the local connectivity (PointIds) array. /// /// Assumes that the grid has an i,j,k structure @@ -748,7 +756,7 @@ class MeshBuilder /// \param array of 3 integers specifying the order along each axis of the hexahedron /// ///////////////////////////////////////////////////////////////////////////// - int PointIndexFromIJK(int i, int j, int k, const int* order) + int PointIndexFromIJK(int i, int j, int k, const int* order) const { bool ibdy = (i == 0 || i == order[0]); bool jbdy = (j == 0 || j == order[1]); @@ -859,7 +867,7 @@ class MeshWriter std::cout << "\t" << pair.first << std::endl; } throw std::runtime_error("**** MESH OUTPUT TYPE NOT SUPPORTED ****"); - return; + // return; } } @@ -879,13 +887,13 @@ class MeshWriter /// ///////////////////////////////////////////////////////////////////////////// void write_ensight( - mesh_t& mesh, - elem_t& elem, - node_t& node, - corner_t& corner, - simulation_parameters_t& sim_param, - double time_value, - CArray graphics_times) + mesh_t& mesh, + elem_t& elem, + node_t& node, + corner_t& corner, + simulation_parameters_t& sim_param, + double time_value, + CArray graphics_times) { // Update host data elem.den.update_host(); @@ -1205,17 +1213,16 @@ class MeshWriter /// ///////////////////////////////////////////////////////////////////////////// void write_vtk( - mesh_t& mesh, - elem_t& elem, - node_t& node, - corner_t& corner, - simulation_parameters_t& sim_param, - double time_value, - CArray graphics_times) + mesh_t& mesh, + elem_t& elem, + node_t& node, + corner_t& corner, + simulation_parameters_t& sim_param, + double time_value, + CArray graphics_times) { // Not yet supported throw std::runtime_error("**** VTK OUTPUT TYPE NOT YET SUPPORTED ****"); - } }; diff --git a/single-node-refactor/src/common/material.h b/single-node-refactor/src/common/material.h index c1b270771..8220faace 100644 --- a/single-node-refactor/src/common/material.h +++ b/single-node-refactor/src/common/material.h @@ -1,36 +1,36 @@ /********************************************************************************************** - © 2020. Triad National Security, LLC. All rights reserved. - This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos - National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. - Department of Energy/National Nuclear Security Administration. All rights in the program are - reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear - Security Administration. The Government is granted for itself and others acting on its behalf a - nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare - derivative works, distribute copies to the public, perform publicly and display publicly, and - to permit others to do so. - This program is open source under the BSD-3 License. - Redistribution and use in source and binary forms, with or without modification, are permitted - provided that the following conditions are met: - 1. Redistributions of source code must retain the above copyright notice, this list of - conditions and the following disclaimer. - 2. Redistributions in binary form must reproduce the above copyright notice, this list of - conditions and the following disclaimer in the documentation and/or other materials - provided with the distribution. - 3. Neither the name of the copyright holder nor the names of its contributors may be used - to endorse or promote products derived from this software without specific prior - written permission. - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS - IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR - CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; - OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, - WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR - OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF - ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - **********************************************************************************************/ +© 2020. Triad National Security, LLC. All rights reserved. +This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos +National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. +Department of Energy/National Nuclear Security Administration. All rights in the program are +reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear +Security Administration. The Government is granted for itself and others acting on its behalf a +nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare +derivative works, distribute copies to the public, perform publicly and display publicly, and +to permit others to do so. +This program is open source under the BSD-3 License. +Redistribution and use in source and binary forms, with or without modification, are permitted +provided that the following conditions are met: +1. Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, this list of +conditions and the following disclaimer in the documentation and/or other materials +provided with the distribution. +3. Neither the name of the copyright holder nor the names of its contributors may be used +to endorse or promote products derived from this software without specific prior +written permission. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; +OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +**********************************************************************************************/ #ifndef FIERRO_MATERIAL_H #define FIERRO_MATERIAL_H @@ -41,14 +41,52 @@ namespace model { + // strength model types -enum strength_tag +enum strength_type { - none = 0, - hypo = 1, // hypoelastic plastic model - hyper = 2, // hyperelastic plastic model + no_strength_type, + increment_based, // Model evaluation is independent of time integration + state_based, // Model is dependent on time integration +}; + +// Specific strength models +enum strength_models +{ + no_strength_model, + user_defined_strength, +}; + +// EOS model types +enum eos_type +{ + no_eos, + ideal_gas, + user_defined_eos, +}; + +} // end model namespace + +static std::map eos_map +{ + { "no_eos", model::no_eos }, + { "ideal_gas", model::ideal_gas }, + { "user_defined", model::user_defined_eos}, +}; + +static std::map strength_type_map +{ + { "no_strength", model::no_strength_type }, + { "increment_based", model::increment_based }, + { "state_based", model::state_based }, }; -} // end of namespace + +static std::map strength_models_map +{ + { "no_strength", model::no_strength_model }, + { "user_defined_strength", model::user_defined_strength }, +}; + namespace model_init { @@ -60,7 +98,6 @@ enum strength_setup_tag }; } // end of namespace - ///////////////////////////////////////////////////////////////////////////// /// /// \struct material_t @@ -79,6 +116,9 @@ struct material_t // statev(4) = ref density // statev(5) = ref specific internal energy + // Type of EOS model used + model::eos_type eos_type = model::no_eos; + // Equation of state (EOS) function pointer void (*eos_model)(const DCArrayKokkos& elem_pres, const DCArrayKokkos& elem_stress, @@ -89,22 +129,43 @@ struct material_t const double den, const double sie); - // Strength model (EOS) function pointer - void (*strength_model)(double, double); // WARNING: a placeholder + // Strength model type + model::strength_type strength_type = model::no_strength_type; + + // Material strength model function pointer + void (*strength_model)(const DCArrayKokkos& elem_pres, + const DCArrayKokkos& elem_stress, + const size_t elem_gid, + const size_t mat_id, + const DCArrayKokkos& elem_state_vars, + const DCArrayKokkos& elem_sspd, + const double den, + const double sie, + const ViewCArrayKokkos &vel_grad, + const ViewCArrayKokkos &elem_node_gids, + const DCArrayKokkos &node_coords, + const DCArrayKokkos &node_vel, + const double vol, + const double dt, + const double rk_alpha); + + + - // hypo or hyper elastic plastic model - model::strength_tag strength_type; // setup the strength model via the input file for via a user_setup model_init::strength_setup_tag strength_setup = model_init::input; size_t num_eos_state_vars = 0; ///< Number of state variables for the EOS - size_t num_strength_state_vars = 0;///< Number of state variables for the strength model - size_t num_eos_global_vars = 0;///< Number of global variables for the EOS - size_t num_strength_global_vars = 0;///< Number of global variables for the strength model + size_t num_strength_state_vars = 0; ///< Number of state variables for the strength model + size_t num_eos_global_vars = 0; ///< Number of global variables for the EOS + size_t num_strength_global_vars = 0; ///< Number of global variables for the strength model DCArrayKokkos eos_global_vars; ///< Array of global variables for the EOS + + DCArrayKokkos strength_global_vars; ///< Array of global variables for the strength model + double q1 = 1.0; ///< acoustic coefficient in Riemann solver for compression double q1ex = 1.3333; ///< acoustic coefficient in Riemann solver for expansion double q2 = 1.0; ///< linear coefficient in Riemann solver for compression @@ -123,6 +184,7 @@ static std::vector str_material_inps { "id", "eos_model", + "strength_model_type" "strength_model", "q1", "q2", @@ -183,34 +245,111 @@ static void ideal_gas(const DCArrayKokkos& elem_pres, return; } // end of ideal_gas -// WARNING: placeholder -static void elastic_plastic(double stress, double strain) -{ - // do nothing - std::cout << "hello from elastic_plastic! Replace with actual strength model!" << std::endl; -} - -// EOS function pointer -typedef void (*eos_type)(const DCArrayKokkos& elem_pres, - const DCArrayKokkos& elem_stress, - const size_t elem_gid, - const size_t mat_id, - const DCArrayKokkos& elem_state_vars, - const DCArrayKokkos& elem_sspd, - const double den, - const double sie); - -static std::map eos_map + +KOKKOS_FUNCTION +static void no_eos(const DCArrayKokkos& elem_pres, + const DCArrayKokkos& elem_stress, + const size_t elem_gid, + const size_t mat_id, + const DCArrayKokkos& elem_state_vars, + const DCArrayKokkos& elem_sspd, + const double den, + const double sie) { - { "ideal_gas", ideal_gas } -}; + return; +} // end of no_eos + + +// ----------------------------------------------------------------------------- +// This is the user material model function for the equation of state +// An eos function must be supplied or the code will fail to run. +// The pressure and sound speed can be calculated from an analytic eos. +// The pressure can also be calculated using p = -1/3 Trace(Stress) +//------------------------------------------------------------------------------ +KOKKOS_FUNCTION +static void user_eos_model(const DCArrayKokkos& elem_pres, + const DCArrayKokkos& elem_stress, + const size_t elem_gid, + const size_t mat_id, + const DCArrayKokkos& elem_state_vars, + const DCArrayKokkos& elem_sspd, + const double den, + const double sie){ + + // ----------------------------------------------------------------------------- + // Required variables are here + //------------------------------------------------------------------------------ + std::cout<<"Inside user EOS model"<& elem_pres, + const DCArrayKokkos& elem_stress, + const size_t elem_gid, + const size_t mat_id, + const DCArrayKokkos& elem_state_vars, + const DCArrayKokkos& elem_sspd, + const double den, + const double sie, + const ViewCArrayKokkos &vel_grad, + const ViewCArrayKokkos &elem_node_gids, + const DCArrayKokkos &node_coords, + const DCArrayKokkos &node_vel, + const double vol, + const double dt, + const double rk_alpha){ + + // ----------------------------------------------------------------------------- + // Required variables are here + //------------------------------------------------------------------------------ + + std::cout<<"Inside user strength model"<& elem_pres, + const DCArrayKokkos& elem_stress, + const size_t elem_gid, + const size_t mat_id, + const DCArrayKokkos& elem_state_vars, + const DCArrayKokkos& elem_sspd, + const double den, + const double sie, + const ViewCArrayKokkos &vel_grad, + const ViewCArrayKokkos &elem_node_gids, + const DCArrayKokkos &node_coords, + const DCArrayKokkos &node_vel, + const double vol, + const double dt, + const double rk_alpha){ + + return; + +} // end of user mat -// add the strength models here -typedef void (*strength_ptr)(double, double); -static std::map strength_map -{ - { "elastic_plastic", elastic_plastic } -}; #endif // end Header Guard \ No newline at end of file diff --git a/single-node-refactor/src/common/mesh.h b/single-node-refactor/src/common/mesh.h index 7d0624684..96920b0a6 100644 --- a/single-node-refactor/src/common/mesh.h +++ b/single-node-refactor/src/common/mesh.h @@ -864,7 +864,7 @@ struct mesh_t /// REMOVE TO SETUP /// ///////////////////////////////////////////////////////////////////////////// - void build_boundry_node_sets(const CArrayKokkos& boundary, + void build_boundry_node_sets(const DCArrayKokkos& boundary, mesh_t& mesh) { // build boundary nodes in each boundary set diff --git a/single-node-refactor/src/common/mesh_inputs.h b/single-node-refactor/src/common/mesh_inputs.h index f0f272270..75cf4a03b 100644 --- a/single-node-refactor/src/common/mesh_inputs.h +++ b/single-node-refactor/src/common/mesh_inputs.h @@ -1,41 +1,42 @@ /********************************************************************************************** - © 2020. Triad National Security, LLC. All rights reserved. - This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos - National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. - Department of Energy/National Nuclear Security Administration. All rights in the program are - reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear - Security Administration. The Government is granted for itself and others acting on its behalf a - nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare - derivative works, distribute copies to the public, perform publicly and display publicly, and - to permit others to do so. - This program is open source under the BSD-3 License. - Redistribution and use in source and binary forms, with or without modification, are permitted - provided that the following conditions are met: - 1. Redistributions of source code must retain the above copyright notice, this list of - conditions and the following disclaimer. - 2. Redistributions in binary form must reproduce the above copyright notice, this list of - conditions and the following disclaimer in the documentation and/or other materials - provided with the distribution. - 3. Neither the name of the copyright holder nor the names of its contributors may be used - to endorse or promote products derived from this software without specific prior - written permission. - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS - IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR - CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; - OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, - WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR - OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF - ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - **********************************************************************************************/ +© 2020. Triad National Security, LLC. All rights reserved. +This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos +National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. +Department of Energy/National Nuclear Security Administration. All rights in the program are +reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear +Security Administration. The Government is granted for itself and others acting on its behalf a +nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare +derivative works, distribute copies to the public, perform publicly and display publicly, and +to permit others to do so. +This program is open source under the BSD-3 License. +Redistribution and use in source and binary forms, with or without modification, are permitted +provided that the following conditions are met: +1. Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, this list of +conditions and the following disclaimer in the documentation and/or other materials +provided with the distribution. +3. Neither the name of the copyright holder nor the names of its contributors may be used +to endorse or promote products derived from this software without specific prior +written permission. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; +OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +**********************************************************************************************/ #ifndef FIERRO_MESH_INPUT_OPTIONS_H #define FIERRO_MESH_INPUT_OPTIONS_H #include #include "matar.h" +#include "parse_yaml.h" namespace mesh_input { @@ -51,7 +52,7 @@ enum source enum type { Box = 0, // Create the mesh using the mesh builder - Cylinder = 1, // Read in the mesh from a file + Cylinder = 1, // Read in the mesh from a file }; } // end of namespace @@ -81,9 +82,10 @@ struct mesh_input_t std::string file_path = ""; ///< Absolute path of mesh file mesh_input::type type; ///< Type of mesh to generate if - std::vector origin { 0.0, 0.0, 0.0 }; ///< Mesh origin for generating a mesh - std::vector length { 1.0, 1.0, 1.0 }; ///< x,y,z length of generated mesh - std::vector num_elems { 2, 2, 2 }; ///< Number of elements along x,y, z for generating a mesh. + DCArrayKokkos origin; ///< Mesh origin for generating a mesh + DCArrayKokkos length; ///< x,y,z length of generated mesh + DCArrayKokkos num_elems; ///< Number of elements along x,y, z for generating a mesh. + size_t p_order = 1; // WARNING, NOT YET PARSED @@ -92,8 +94,8 @@ struct mesh_input_t double starting_angle = 0.0; ///< Starting angle in degrees for 2D RZ mesh double ending_angle = 90; ///< Ending angle in degrees for 2D RZ mesh - int num_radial_elems = 10; ///< Number of elements in the radial direction for 2DRZ mesh - int num_angular_elems = 10; ///< Number of elements in the radial direction for 2DRZ mesh + int num_radial_elems = 10; ///< Number of elements in the radial direction for 2DRZ mesh + int num_angular_elems = 10; ///< Number of elements in the radial direction for 2DRZ mesh }; // mesh_input_t // ---------------------------------- @@ -117,4 +119,5 @@ static std::vector str_mesh_inps "num_angular_elems" }; + #endif // end Header Guard \ No newline at end of file diff --git a/single-node-refactor/src/common/region.h b/single-node-refactor/src/common/region.h index 364fec26b..c5808b093 100644 --- a/single-node-refactor/src/common/region.h +++ b/single-node-refactor/src/common/region.h @@ -1,5 +1,5 @@ /********************************************************************************************** - © 2020. Triad National Security, LLC. All rights reserved. + � 2020. Triad National Security, LLC. All rights reserved. This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. Department of Energy/National Nuclear Security Administration. All rights in the program are @@ -47,7 +47,7 @@ namespace region // for tagging boundary faces enum vol_tag { - global = 0, // tag every elements in the mesh + global = 0, // tag every element in the mesh box = 1, // tag all elements inside a box cylinder = 2, // tag all elements inside a cylinder sphere = 3, // tag all elements inside a sphere @@ -86,8 +86,8 @@ struct reg_fill_t double x2 = 0.0; ///< Second X plane for creating a box double y1 = 0.0; ///< First Y plane for creating a box double y2 = 0.0; ///< Second Y plane for creating a box - double z1 = 0.0; ///< First Z plane for creating a box - double z2 = 0.0; ///< Second Z plane for creating a box + double z1 = 0.0; ///< First Z plane for creating a box + double z2 = 0.0; ///< Second Z plane for creating a box // radius double radius1 = 0.0; ///< Inner radius to fill for sphere @@ -101,14 +101,13 @@ struct reg_fill_t double v = 0.0; ///< V component of velocity double w = 0.0; ///< W component of velocity - double speed = 0.0; ///< velocity magnitude for radial velocity initialization double ie = 0.0; ///< extensive internal energy double sie = 0.0; ///< specific internal energy double den = 0.0; ///< density - std::vector origin = { 0.0, 0.0, 0.0 }; ///< Origin for region + DCArrayKokkos origin; ///< Origin for region }; // ---------------------------------- diff --git a/single-node-refactor/src/common/simulation_parameters.h b/single-node-refactor/src/common/simulation_parameters.h index 95fd89831..b913ff727 100644 --- a/single-node-refactor/src/common/simulation_parameters.h +++ b/single-node-refactor/src/common/simulation_parameters.h @@ -62,14 +62,11 @@ struct simulation_parameters_t std::vector solver_inputs; ///< Solvers to use during the simulation - CArrayKokkos boundary_conditions; ///< Simulation boundary conditions + DCArrayKokkos boundary_conditions; ///< Simulation boundary conditions - CArrayKokkos region_fills; ///< Region data for simulation mesh - - CArrayKokkos materials; ///< Material data for simulation - - std::vector> eos_global_vars; ///< EOS data for simulation + DCArrayKokkos region_fills; ///< Region data for simulation mesh + DCArrayKokkos materials; ///< Material data for simulation }; // simulation_parameters_t #endif // end Header Guard \ No newline at end of file diff --git a/single-node-refactor/src/driver.h b/single-node-refactor/src/driver.h index 83748a6a7..7526788aa 100644 --- a/single-node-refactor/src/driver.h +++ b/single-node-refactor/src/driver.h @@ -248,9 +248,9 @@ class Driver elem_coords[2] = 0.0; } } // end loop over nodes in element (NOTE: Translating using origin so below check works) - elem_coords[0] = (elem_coords[0] / mesh.num_nodes_in_elem) - sim_param.region_fills(f_id).origin[0]; - elem_coords[1] = (elem_coords[1] / mesh.num_nodes_in_elem) - sim_param.region_fills(f_id).origin[1]; - elem_coords[2] = (elem_coords[2] / mesh.num_nodes_in_elem) - sim_param.region_fills(f_id).origin[2]; + elem_coords[0] = (elem_coords[0] / mesh.num_nodes_in_elem) - sim_param.region_fills(f_id).origin(0); + elem_coords[1] = (elem_coords[1] / mesh.num_nodes_in_elem) - sim_param.region_fills(f_id).origin(1); + elem_coords[2] = (elem_coords[2] / mesh.num_nodes_in_elem) - sim_param.region_fills(f_id).origin(2); // spherical radius double radius = sqrt(elem_coords[0] * elem_coords[0] + diff --git a/single-node-refactor/src/setup/parse_yaml.cpp b/single-node-refactor/src/setup/parse_yaml.cpp index 505e65884..c49e1ac8e 100644 --- a/single-node-refactor/src/setup/parse_yaml.cpp +++ b/single-node-refactor/src/setup/parse_yaml.cpp @@ -1,36 +1,36 @@ /********************************************************************************************** - © 2020. Triad National Security, LLC. All rights reserved. - This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos - National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. - Department of Energy/National Nuclear Security Administration. All rights in the program are - reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear - Security Administration. The Government is granted for itself and others acting on its behalf a - nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare - derivative works, distribute copies to the public, perform publicly and display publicly, and - to permit others to do so. - This program is open source under the BSD-3 License. - Redistribution and use in source and binary forms, with or without modification, are permitted - provided that the following conditions are met: - 1. Redistributions of source code must retain the above copyright notice, this list of - conditions and the following disclaimer. - 2. Redistributions in binary form must reproduce the above copyright notice, this list of - conditions and the following disclaimer in the documentation and/or other materials - provided with the distribution. - 3. Neither the name of the copyright holder nor the names of its contributors may be used - to endorse or promote products derived from this software without specific prior - written permission. - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS - IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR - CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; - OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, - WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR - OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF - ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - **********************************************************************************************/ +© 2020. Triad National Security, LLC. All rights reserved. +This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos +National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. +Department of Energy/National Nuclear Security Administration. All rights in the program are +reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear +Security Administration. The Government is granted for itself and others acting on its behalf a +nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare +derivative works, distribute copies to the public, perform publicly and display publicly, and +to permit others to do so. +This program is open source under the BSD-3 License. +Redistribution and use in source and binary forms, with or without modification, are permitted +provided that the following conditions are met: +1. Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, this list of +conditions and the following disclaimer in the documentation and/or other materials +provided with the distribution. +3. Neither the name of the copyright holder nor the names of its contributors may be used +to endorse or promote products derived from this software without specific prior +written permission. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; +OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +**********************************************************************************************/ #include #include #include @@ -54,7 +54,7 @@ bool VERBOSE = false; // ============================================================================== // Function Definitions // ============================================================================== -// modified code from stackover flow for string delimiter parsing +// modified code from stack over./flow for string delimiter parsing std::vector exact_array_values(std::string s, std::string delimiter) { size_t pos_start = 0, pos_end, delim_len = delimiter.length(); @@ -242,10 +242,10 @@ void print_yaml(Yaml::Node root) // ================================================================================= void parse_yaml(Yaml::Node& root, simulation_parameters_t& sim_param) { - if (VERBOSE) { - printf("\n"); - std::cout << "Printing YAML Input file:" << std::endl; - } + // if (VERBOSE) { + // printf("\n"); + // std::cout << "Printing YAML Input file:" << std::endl; + // } // print the input file print_yaml(root); @@ -291,7 +291,7 @@ void parse_yaml(Yaml::Node& root, simulation_parameters_t& sim_param) std::cout << "Parsing YAML materials:" << std::endl; } // parse the material yaml text into a vector of materials - parse_materials(root, sim_param.materials, sim_param.eos_global_vars); + parse_materials(root, sim_param.materials); } // ================================================================================= @@ -482,6 +482,11 @@ void parse_mesh_input(Yaml::Node& root, mesh_input_t& mesh_input) // get the mesh variables names set by the user std::vector user_mesh_inputs; + mesh_input.origin = DCArrayKokkos (3, "mesh_input.origin"); + mesh_input.length = DCArrayKokkos (3, "mesh_input.length"); + mesh_input.num_elems = DCArrayKokkos (3, "mesh_input.num_elems"); + + // extract words from the input file and validate they are correct validate_inputs(mesh_yaml, user_mesh_inputs, str_mesh_inps); @@ -523,7 +528,6 @@ void parse_mesh_input(Yaml::Node& root, mesh_input_t& mesh_input) // Number of dimensions for the mesh else if (a_word.compare("num_dims") == 0) { int num_dim = root["mesh_options"][a_word].As(); - std::cout << "\tNTEST TEST TEST NUM DIM = " << num_dim << std::endl; if (VERBOSE) { std::cout << "\tNum dimensions = " << num_dim << std::endl; } @@ -581,45 +585,75 @@ void parse_mesh_input(Yaml::Node& root, mesh_input_t& mesh_input) // get the origin numbers, values are words std::vector numbers = exact_array_values(origin, ","); - std::vector val; - for (int i = 0; i < 3; i++) { - val.push_back(std::stod(numbers[i])); + double x1 = std::stod(numbers[0]); + double y1 = std::stod(numbers[1]); + double z1 = std::stod(numbers[2]); + + if (VERBOSE) { + std::cout << "\tx1 = " << x1 << std::endl; + std::cout << "\ty1 = " << y1 << std::endl; + std::cout << "\tz1 = " << z1 << std::endl; } - mesh_input.origin = val; + // storing the origin values as + RUN({ + mesh_input.origin(0) = x1; + mesh_input.origin(1) = y1; + mesh_input.origin(2) = z1; + }); } // Extents of the mesh else if (a_word.compare("length") == 0) { - std::string origin = root["mesh_options"][a_word].As(); + std::string length = root["mesh_options"][a_word].As(); if (VERBOSE) { - std::cout << "\tlength = " << origin << std::endl; + std::cout << "\tlength = " << length << std::endl; } - std::vector numbers = exact_array_values(origin, ","); + std::vector numbers = exact_array_values(length, ","); + + double l1 = std::stod(numbers[0]); + double l2 = std::stod(numbers[1]); + double l3 = std::stod(numbers[2]); - std::vector val; - for (int i = 0; i < 3; i++) { - val.push_back(std::stod(numbers[i])); + if (VERBOSE) { + std::cout << "\tl1 = " << l1 << std::endl; + std::cout << "\tl2 = " << l2 << std::endl; + std::cout << "\tl3 = " << l3 << std::endl; } - mesh_input.length = val; + // storing the length values + RUN({ + mesh_input.length(0) = l1; + mesh_input.length(1) = l2; + mesh_input.length(2) = l3; + }); } // Number of elements per direction else if (a_word.compare("num_elems") == 0) { - std::string origin = root["mesh_options"][a_word].As(); + std::string elem_count = root["mesh_options"][a_word].As(); if (VERBOSE) { - std::cout << "\tnum_elems = " << origin << std::endl; + std::cout << "\tnum_elems = " << elem_count << std::endl; } - // get the origin numbers, values are words - std::vector numbers = exact_array_values(origin, ","); + // get the elem_count numbers, values are words + std::vector numbers = exact_array_values(elem_count, ","); + + int n1 = std::stod(numbers[0]); + int n2 = std::stod(numbers[1]); + int n3 = std::stod(numbers[2]); - std::vector val; - for (int i = 0; i < 3; i++) { - val.push_back(std::stoi(numbers[i])); + if (VERBOSE) { + std::cout << "\tn1 = " << n1 << std::endl; + std::cout << "\tn2 = " << n2 << std::endl; + std::cout << "\tn3 = " << n3 << std::endl; } - mesh_input.num_elems = val; + // storing the number of elements + RUN({ + mesh_input.num_elems(0) = n1; + mesh_input.num_elems(1) = n2; + mesh_input.num_elems(2) = n3; + }); } // Polynomial order for the mesh else if (a_word.compare("polynomial_order") == 0) { @@ -789,13 +823,19 @@ void parse_output_options(Yaml::Node& root, output_options_t& output_options) // ================================================================================= // Parse Fill regions // ================================================================================= -void parse_regions(Yaml::Node& root, CArrayKokkos& region_fills) +void parse_regions(Yaml::Node& root, DCArrayKokkos& region_fills) { Yaml::Node& region_yaml = root["regions"]; size_t num_regions = region_yaml.Size(); - region_fills = CArrayKokkos(num_regions); + region_fills = DCArrayKokkos(num_regions , "sim_param.region_fills"); + + for(int i=0; i< num_regions; i++){ + region_fills.host(i).origin = DCArrayKokkos (3, "region_fills.origin"); + region_fills.host(i).origin.update_device(); + } + region_fills.update_device(); // loop over the fill regions specified for (int reg_id = 0; reg_id < num_regions; reg_id++) { @@ -1067,9 +1107,9 @@ void parse_regions(Yaml::Node& root, CArrayKokkos& region_fills) // storing the origin values as (x1,y1,z1) RUN({ - region_fills(reg_id).origin[0] = x1; - region_fills(reg_id).origin[1] = y1; - region_fills(reg_id).origin[2] = z1; + region_fills(reg_id).origin(0) = x1; + region_fills(reg_id).origin(1) = y1; + region_fills(reg_id).origin(2) = z1; }); } // origin else { @@ -1086,17 +1126,15 @@ void parse_regions(Yaml::Node& root, CArrayKokkos& region_fills) // ================================================================================= // Parse Material Definitions // ================================================================================= -void parse_materials(Yaml::Node& root, CArrayKokkos& materials, - std::vector>& eos_global_vars) +void parse_materials(Yaml::Node& root, DCArrayKokkos& materials) { Yaml::Node& material_yaml = root["materials"]; size_t num_materials = material_yaml.Size(); - materials = CArrayKokkos(num_materials); + materials = DCArrayKokkos(num_materials, "sim_param.materials"); // allocate room for each material to store eos_global_vars - eos_global_vars = std::vector>(num_materials); // loop over the materials specified for (int mat_id = 0; mat_id < num_materials; mat_id++) { @@ -1160,7 +1198,6 @@ void parse_materials(Yaml::Node& root, CArrayKokkos& materials, if (VERBOSE) { std::cout << "\tq2ex = " << q2ex << std::endl; } - RUN({ materials(mat_id).q2ex = q2ex; }); @@ -1197,49 +1234,152 @@ void parse_materials(Yaml::Node& root, CArrayKokkos& materials, // set the EOS if (eos_map.find(eos) != eos_map.end()) { - auto eos_model = eos_map[eos]; + + switch(eos_map[eos]){ + + case model::no_eos: + RUN({ + materials(mat_id).eos_type = model::no_eos; + materials(mat_id).eos_model = no_eos; + }); + if (VERBOSE) { + std::cout << "\teos_model = " << eos << std::endl; + } + break; + + case model::ideal_gas: + RUN({ + materials(mat_id).eos_type = model::ideal_gas; + materials(mat_id).eos_model = ideal_gas; + }); + if (VERBOSE) { + std::cout << "\teos_model = " << eos << std::endl; + } + break; + case model::user_defined_eos: + RUN({ + materials(mat_id).eos_type = model::user_defined_eos; + materials(mat_id).eos_model = user_eos_model; + }); + if (VERBOSE) { + std::cout << "\teos_model = " << eos << std::endl; + } + break; + default: + std::cout << "ERROR: invalid input: " << eos << std::endl; + throw std::runtime_error("**** EOS Not Understood ****"); + break; + } // end switch on EOS type + } + else{ + std::cout << "ERROR: invalid EOS input: " << eos << std::endl; + throw std::runtime_error("**** EOS Not Understood ****"); + } // end if + } // EOS model - RUN({ - materials(mat_id).eos_model = eos_model; - // materials(mat_id).eos_model(0.,1.,2.); // WARNING BUG HERE, replace with real EOS model - }); + // Type of strength model + else if (a_word.compare("strength_model_type") == 0) { + std::string strength_model_type = root["materials"][mat_id]["material"]["strength_model_type"].As(); - if (VERBOSE) { - std::cout << "\teos_model = " << eos << std::endl; - } + // set the EOS + if (strength_type_map.find(strength_model_type) != strength_type_map.end()) { + + switch(strength_type_map[strength_model_type]){ + + case model::no_strength_type: + RUN({ + materials(mat_id).strength_type = model::no_strength_type; + }); + if (VERBOSE) { + std::cout << "\tstrength_model_type_type = " << strength_model_type << std::endl; + } + break; + + case model::increment_based: + RUN({ + materials(mat_id).strength_type = model::increment_based; + }); + + if (VERBOSE) { + std::cout << "\tstrength_model_type = " << strength_model_type << std::endl; + } + break; + case model::state_based: + RUN({ + materials(mat_id).strength_type = model::state_based; + }); + std::cout << "ERROR: state_based models not yet defined: " << std::endl; + throw std::runtime_error("**** ERROR: state_based models not yet defined ****"); + if (VERBOSE) { + std::cout << "\tstrength_model_type = " << strength_model_type << std::endl; + } + break; + default: + std::cout << "ERROR: invalid strength type input: " << strength_model_type << std::endl; + throw std::runtime_error("**** Strength Model Type Not Understood ****"); + break; + } // end switch on EOS type } else{ - std::cout << "ERROR: invalid input: " << eos << std::endl; + std::cout << "ERROR: Invalid strength model type input: " << strength_model_type << std::endl; + throw std::runtime_error("**** Strength model type not understood ****"); } // end if - } // EOS model + + } // Strength model type + + // Set specific strength model else if (a_word.compare("strength_model") == 0) { std::string strength_model = root["materials"][mat_id]["material"]["strength_model"].As(); - // set the strength_model - if (strength_map.find(strength_model) != strength_map.end()) { - auto strength = strength_map[strength_model]; - - RUN({ - materials(mat_id).strength_model = strength; - materials(mat_id).strength_model(0., 1.); // WARNING BUG HERE, replace with real strength model - }); - - if (VERBOSE) { - std::cout << "\tstrength_model = " << strength_model << std::endl; - } + // set the EOS + if (strength_models_map.find(strength_model) != strength_models_map.end()) { + + switch(strength_models_map[strength_model]){ + + case model::no_strength_model: + RUN({ + materials(mat_id).strength_model = no_strength; + }); + if (VERBOSE) { + std::cout << "\tstrength_model = " << strength_model << std::endl; + } + break; + + case model::user_defined_strength: + RUN({ + materials(mat_id).strength_model = user_strength_model; + }); + if (VERBOSE) { + std::cout << "\tstrength_model = " << strength_model << std::endl; + } + break; + default: + std::cout << "ERROR: invalid strength input: " << strength_model << std::endl; + throw std::runtime_error("**** Strength model Not Understood ****"); + break; + } // end switch on EOS type } else{ - std::cout << "ERROR: invalid input: " << strength_model << std::endl; + std::cout << "ERROR: invalid Strength model input: " << strength_model << std::endl; + throw std::runtime_error("**** Strength model Not Understood ****"); } // end if - } // EOS model + + } // Strength model // exact the eos_global_vars else if (a_word.compare("eos_global_vars") == 0) { - // Yaml::Node & mat_global_vars_yaml = root["materials"][mat_id]["material"][a_word]; + Yaml::Node & mat_global_vars_yaml = root["materials"][mat_id]["material"][a_word]; - size_t num_global_vars = material_inps_yaml.Size(); + size_t num_global_vars = mat_global_vars_yaml.Size(); std::cout << "*** parsing num global eos vars = " << num_global_vars << std::endl; - materials(mat_id).eos_global_vars = DCArrayKokkos(num_global_vars); + + + materials.host(mat_id).eos_global_vars = DCArrayKokkos(num_global_vars, "material.eos_global_vars"); + materials.host(mat_id).eos_global_vars.update_device(); + RUN({ + materials(mat_id).num_eos_global_vars = num_global_vars; + }); + materials.update_device(); if (VERBOSE) { std::cout << "num global eos vars = " << num_global_vars << std::endl; @@ -1247,7 +1387,7 @@ void parse_materials(Yaml::Node& root, CArrayKokkos& materials, for (int global_var_id = 0; global_var_id < num_global_vars; global_var_id++) { double eos_var = root["materials"][mat_id]["material"]["eos_global_vars"][global_var_id].As(); - eos_global_vars[mat_id].push_back(eos_var); + RUN({ materials(mat_id).eos_global_vars(global_var_id) = eos_var; @@ -1272,13 +1412,19 @@ void parse_materials(Yaml::Node& root, CArrayKokkos& materials, // ================================================================================= // Parse Boundary Conditions // ================================================================================= -void parse_bcs(Yaml::Node& root, CArrayKokkos& boundary_conditions) +void parse_bcs(Yaml::Node& root, DCArrayKokkos& boundary_conditions) { Yaml::Node& bc_yaml = root["boundary_conditions"]; size_t num_bcs = bc_yaml.Size(); - boundary_conditions = CArrayKokkos(num_bcs); + boundary_conditions = DCArrayKokkos(num_bcs, "sim_param.boundary_conditions"); + + for(int i=0; i< num_bcs; i++){ + boundary_conditions.host(i).origin = DCArrayKokkos (3, "boundary_conditions.origin"); + boundary_conditions.host(i).origin.update_device(); + } + boundary_conditions.update_device(); // loop over the fill regions specified for (int bc_id = 0; bc_id < num_bcs; bc_id++) { @@ -1307,7 +1453,7 @@ void parse_bcs(Yaml::Node& root, CArrayKokkos& boundary_co // set the solver if (map.find(solver) != map.end()) { - auto bc_solver = map[solver]; + solver_input::method bc_solver = map[solver]; RUN({ boundary_conditions(bc_id).solver = bc_solver; @@ -1451,9 +1597,9 @@ void parse_bcs(Yaml::Node& root, CArrayKokkos& boundary_co // storing the origin values as (x1,y1,z1) RUN({ - boundary_conditions(bc_id).origin[0] = x1; - boundary_conditions(bc_id).origin[1] = y1; - boundary_conditions(bc_id).origin[2] = z1; + boundary_conditions(bc_id).origin(0) = x1; + boundary_conditions(bc_id).origin(1) = y1; + boundary_conditions(bc_id).origin(2) = z1; }); } // origin else { diff --git a/single-node-refactor/src/setup/parse_yaml.h b/single-node-refactor/src/setup/parse_yaml.h index d68303e96..1b26c945b 100644 --- a/single-node-refactor/src/setup/parse_yaml.h +++ b/single-node-refactor/src/setup/parse_yaml.h @@ -1,36 +1,36 @@ /********************************************************************************************** - © 2020. Triad National Security, LLC. All rights reserved. - This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos - National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. - Department of Energy/National Nuclear Security Administration. All rights in the program are - reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear - Security Administration. The Government is granted for itself and others acting on its behalf a - nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare - derivative works, distribute copies to the public, perform publicly and display publicly, and - to permit others to do so. - This program is open source under the BSD-3 License. - Redistribution and use in source and binary forms, with or without modification, are permitted - provided that the following conditions are met: - 1. Redistributions of source code must retain the above copyright notice, this list of - conditions and the following disclaimer. - 2. Redistributions in binary form must reproduce the above copyright notice, this list of - conditions and the following disclaimer in the documentation and/or other materials - provided with the distribution. - 3. Neither the name of the copyright holder nor the names of its contributors may be used - to endorse or promote products derived from this software without specific prior - written permission. - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS - IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR - CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; - OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, - WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR - OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF - ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - **********************************************************************************************/ +© 2020. Triad National Security, LLC. All rights reserved. +This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos +National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. +Department of Energy/National Nuclear Security Administration. All rights in the program are +reserved by Triad National Security, LLC, and the U.S. Department of Energy/National Nuclear +Security Administration. The Government is granted for itself and others acting on its behalf a +nonexclusive, paid-up, irrevocable worldwide license in this material to reproduce, prepare +derivative works, distribute copies to the public, perform publicly and display publicly, and +to permit others to do so. +This program is open source under the BSD-3 License. +Redistribution and use in source and binary forms, with or without modification, are permitted +provided that the following conditions are met: +1. Redistributions of source code must retain the above copyright notice, this list of +conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, this list of +conditions and the following disclaimer in the documentation and/or other materials +provided with the distribution. +3. Neither the name of the copyright holder nor the names of its contributors may be used +to endorse or promote products derived from this software without specific prior +written permission. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; +OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +**********************************************************************************************/ #ifndef FIERRO_PARSE_YAML_H #define FIERRO_PARSE_YAML_H @@ -40,6 +40,7 @@ #include #include #include +#include #include "matar.h" @@ -91,14 +92,12 @@ void parse_mesh_input(Yaml::Node& root, mesh_input_t& mesh_input); void parse_output_options(Yaml::Node& root, output_options_t& output_options); // parse the region text -void parse_regions(Yaml::Node& root, CArrayKokkos& region_fills); +void parse_regions(Yaml::Node& root, DCArrayKokkos& region_fills); // parse the region text -void parse_materials(Yaml::Node& root, - CArrayKokkos& materials, - std::vector>& eos_global_vars); +void parse_materials(Yaml::Node& root, DCArrayKokkos& materials); // parse the boundary condition text -void parse_bcs(Yaml::Node& root, CArrayKokkos& boundary_conditions); +void parse_bcs(Yaml::Node& root, DCArrayKokkos& boundary_conditions); #endif // end Header Guard \ No newline at end of file