Merge pull request #75 from hugary1995/stokes

Fix mms
Argonne-National-Laboratory · Apr 5, 2024 · 9c827c4 · 9c827c4
2 parents 5d5fa91 + c758bb9
commit 9c827c4
Show file tree

Hide file tree

Showing 3 changed files with 120 additions and 128 deletions.
diff --git a/src/kernels/ADStokesIncompressibility.C b/src/kernels/ADStokesIncompressibility.C
@@ -30,5 +30,5 @@ ADStokesIncompressibility::ADStokesIncompressibility(const InputParameters & par
 ADReal
 ADStokesIncompressibility::computeQpResidual()
 {
-  return -(_grad_vel[_qp](0, 0) + _grad_vel[_qp](1, 1) + _grad_vel[_qp](2, 2)) * _test[_i][_qp];
-}
+  return -_grad_vel[_qp].tr() * _test[_i][_qp];
+}
diff --git a/src/kernels/ADStokesStressDivergence.C b/src/kernels/ADStokesStressDivergence.C
@@ -32,7 +32,5 @@ ADStokesStressDivergence::ADStokesStressDivergence(const InputParameters & param
 ADReal
 ADStokesStressDivergence::computeQpResidual()
 {
-  return _stress[_qp].contract(_grad_test[_i][_qp]) -
-         (_grad_test[_i][_qp](0, 0) + _grad_test[_i][_qp](1, 1) + _grad_test[_i][_qp](2, 2)) *
-             _pressure[_qp];
-}
+  return _stress[_qp].contract(_grad_test[_i][_qp]) - _grad_test[_i][_qp].tr() * _pressure[_qp];
+}
diff --git a/tests/stokes/manufactured.i b/tests/stokes/manufactured.i
@@ -8,59 +8,59 @@
     ymax = 1
     zmin = 0
     zmax = 1
-    nx = 6
-    ny = 6
-    nz = 6
+    nx = 2
+    ny = 2
+    nz = 2
     elem_type = HEX20
   []
 []
 
 [Variables]
-    [p]
-        order = FIRST
-        family = LAGRANGE
-    []
-    [u]
-        order = SECOND
-        family = LAGRANGE_VEC
-    []
+  [p]
+    order = FIRST
+    family = LAGRANGE
+  []
+  [u]
+    order = SECOND
+    family = LAGRANGE_VEC
+  []
 []
 
 [AuxVariables]
-    [u_exact]
-        order = FIRST
-        family = LAGRANGE_VEC
-    []
-    [p_exact]
-        order = FIRST
-        family = LAGRANGE
-    []
+  [u_exact]
+    order = FIRST
+    family = LAGRANGE_VEC
+  []
+  [p_exact]
+    order = FIRST
+    family = LAGRANGE
+  []
 []
 
 [AuxKernels]
-    [u_exact]
-        type = VectorFunctionAux
-        variable = u_exact
-        function = u_exact
-    []
-    [p_exact]
-        type = FunctionAux
-        variable = p_exact
-        function = p_exact
-    []
+  [u_exact]
+    type = VectorFunctionAux
+    variable = u_exact
+    function = u_exact
+  []
+  [p_exact]
+    type = FunctionAux
+    variable = p_exact
+    function = p_exact
+  []
 []
 
 [Functions]
-    [u_exact]
-        type = ParsedVectorFunction
-        expression_x = '2 * x^2 + y^2 + z^2'
-        expression_y = '2 * x^2 - 2 * x * y'
-        expression_z = '2 * x^2 - 2 * x * z'
-    []
-    [p_exact]
-        type = ParsedFunction
-        expression = 'x + y + z - 3.0/2.0'
-    []
+  [u_exact]
+    type = ParsedVectorFunction
+    expression_x = '2 * x^2 + y^2 + z^2'
+    expression_y = '2 * x^2 - 2 * x * y'
+    expression_z = '2 * x^2 - 2 * x * z'
+  []
+  [p_exact]
+    type = ParsedFunction
+    expression = 'x + y + z - 3.0/2.0'
+  []
 []
 
 [ICs]
@@ -74,111 +74,105 @@
 []
 
 [Kernels]
-    [equil]
-        type = ADStokesStressDivergence
-        variable = u
-        pressure = p
-    []
-    [incompressible]
-        type = ADStokesIncompressibility
-        variable = p
-        velocity = u
-    []
-    [body]
-        type = VectorBodyForce
-        variable = u
-        function_x = 7
-        function_y = 3
-        function_z = 3
-    []
+  [equil]
+    type = ADStokesStressDivergence
+    variable = u
+    pressure = p
+  []
+  [incompressible]
+    type = ADStokesIncompressibility
+    variable = p
+    velocity = u
+  []
+  [body]
+    type = VectorBodyForce
+    variable = u
+    function_x = -7
+    function_y = -3
+    function_z = -3
+  []
 []
 
 [Materials]
-    [strain]
-        type = StokesStrainRate
-        velocity = u
-    []
-    [stress]
-        type = StokesLinearViscous
-        mu = 1.0
-    []
+  [strain]
+    type = StokesStrainRate
+    velocity = u
+  []
+  [stress]
+    type = StokesLinearViscous
+    mu = 1.0
+  []
 []
 
 [BCs]
-    [fixed]
-        type = ADVectorFunctionDirichletBC
-        variable = u
-        boundary = 'bottom top back front'
-        function_x = '2 * x^2 + y^2 + z^2'
-        function_y = '2 * x^2 - 2 * x * y'
-        function_z = '2 * x^2 - 2 * x * z'
-    []
-#    [natural]
-#        type = ADVectorFunctionNeumannBC
-#        variable = u
-#        boundary = 'right'
-#        function_x = 7
-#        function_y = 3
-#        function_z = 3
-#    []
-#    [pressure]
-#        type = ADFunctionDirichletBC
-#        boundary = 'left'
-#        variable = p
-#        function = p_exact
-#    []
-[]
-
-[Preconditioning]
-  [FSP]
-    type = FSP
-    topsplit = 'up'
-    [up]
-      splitting = 'u p'
-      splitting_type  = schur
-      petsc_options_iname = '-pc_fieldsplit_schur_fact_type  -pc_fieldsplit_schur_precondition -ksp_gmres_restart -ksp_type -ksp_pc_side'
-      petsc_options_value = 'full                            self                              300                fgmres    right'
-    []
-    [u]
-        vars = 'u'
-        petsc_options_iname = '-pc_type'
-        petsc_options_value = 'lu' 
-    []
-    [p]
-        vars = 'p'
-        petsc_options = '-pc_lsc_scale_diag'
-        petsc_options_iname = '-ksp_type -ksp_gmres_restart -pc_type -ksp_pc_side -lsc_pc_type -lsc_ksp_type -lsc_ksp_pc_side -lsc_ksp_gmres_restart'
-        petsc_options_value = 'fgmres    300                lsc      right        lu            gmres        right            300'
-    []
+  [u_fixed]
+    type = ADVectorFunctionDirichletBC
+    variable = u
+    boundary = 'bottom top back front left right'
+    function_x = '2 * x^2 + y^2 + z^2'
+    function_y = '2 * x^2 - 2 * x * y'
+    function_z = '2 * x^2 - 2 * x * z'
+  []
+  [p_fixed]
+    type = FunctionDirichletBC
+    variable = p
+    boundary = 'bottom top back front left right'
+    function = 'p_exact'
   []
 []
 
+# [Preconditioning]
+#   [FSP]
+#     type = FSP
+#     topsplit = 'up'
+#     [up]
+#       splitting = 'u p'
+#       splitting_type = schur
+#       petsc_options_iname = '-pc_fieldsplit_schur_fact_type  -pc_fieldsplit_schur_precondition -ksp_gmres_restart -ksp_type -ksp_pc_side'
+#       petsc_options_value = 'full                            self                              300                fgmres    right'
+#     []
+#     [u]
+#       vars = 'u'
+#       petsc_options_iname = '-pc_type'
+#       petsc_options_value = 'lu'
+#     []
+#     [p]
+#       vars = 'p'
+#       petsc_options = '-pc_lsc_scale_diag'
+#       petsc_options_iname = '-ksp_type -ksp_gmres_restart -pc_type -ksp_pc_side -lsc_pc_type -lsc_ksp_type -lsc_ksp_pc_side -lsc_ksp_gmres_restart'
+#       petsc_options_value = 'fgmres    300                lsc      right        lu            gmres        right            300'
+#     []
+#   []
+# []
+
 [Executioner]
   type = Steady
 
   solve_type = 'newton'
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'lu'
   line_search = none
 
   l_max_its = 10
   l_tol = 1e-8
   nl_max_its = 15
-  nl_rel_tol = 1e-8
-  nl_abs_tol = 1e-6
+  nl_rel_tol = 1e-10
+  nl_abs_tol = 1e-12
 []
 
 [Postprocessors]
-    [u_error]
-        type = ElementVectorL2Error
-        function = u_exact
-        variable = u
-    []
-    [p_error]
-        type = ElementL2Error
-        function = p_exact
-        variable = p
-    []
+  [u_error]
+    type = ElementVectorL2Error
+    function = u_exact
+    variable = u
+  []
+  [p_error]
+    type = ElementL2Error
+    function = p_exact
+    variable = p
+  []
 []
 
 [Outputs]
-    exodus = true
+  exodus = true
 []