Feature/implicit forcing

doc/pages/user-guide/case-file.md

@@ -30,8 +30,9 @@ The current high-level structure of the case file is shown below.
 }
 ~~~~~~~~~~~~~~~
 The `version` keywords is reserved to track changes in the format of the file.
-The the subsections below we list all the configuration options for each of the high-level objects.
-Some parameters will have default values, and are therefore optional.
+The the subsections below we list all the configuration options for each of the
+high-level objects. Some parameters will have default values, and are therefore
+optional.
 
 ## Output frequency control
 A common scheme for controlling the output frequency is applied for various
@@ -140,11 +141,11 @@ arbitrary parameters read from the case file. Additionally, this allows to
 change the material properties in time.
 
 ### Boundary types {#case-file_boundary-types}
-The optional `boundary_types` keyword can be used to specify boundary conditions.
-The reason for it being optional, is that some conditions can be specified
-directly inside the mesh file.
-In particular, this happens when Nek5000 `.re2` files are converted to `.nmsh`.
-Periodic boundary conditions are *always* defined inside the mesh file.
+The optional `boundary_types` keyword can be used to specify boundary
+conditions. The reason for it being optional, is that some conditions can be
+specified directly inside the mesh file. In particular, this happens when
+Nek5000 `.re2` files are converted to `.nmsh`. Periodic boundary conditions are
+*always* defined inside the mesh file.
 
 The value of the keyword is an array of strings, with the following possible
 values:
@@ -154,8 +155,8 @@ values:
   * `v`, a velocity Dirichlet boundary.
   * `sym`, a symmetry boundary.
   * `o`, outlet boundary.
-  * `on`, Dirichlet for the boundary-parallel velocity and homogeneous Neumann for
-   the wall-normal. The wall-parallel velocity is defined by the initial
+  * `on`, Dirichlet for the boundary-parallel velocity and homogeneous Neumann
+   for the wall-normal. The wall-parallel velocity is defined by the initial
    condition.
 
 * Advanced boundary conditions
@@ -214,10 +215,9 @@ The means of prescribing the values are controlled via the `type` keyword:
    `case.point_zones` object. See more about point zones @ref point-zones.md.
 
 ### Blasius profile
-The `blasius` object is used to specify the Blasius profile that can be used for the
-initial and inflow condition.
-The boundary cannot be tilted with respect to the coordinate axes.
-It requires  the following parameters:
+The `blasius` object is used to specify the Blasius profile that can be used for
+the initial and inflow condition. The boundary cannot be tilted with respect to
+the coordinate axes. It requires  the following parameters:
 
 1. `delta`, the thickness of the boundary layer.
 2. `freestream_velocity`, the velocity value in the free stream.
@@ -268,18 +268,19 @@ The following types are currently implemented.
 
 #### Brinkman
 The Brinkman source term introduces regions of resistance in the fluid domain.
-The volume force \f$ f_i \f$ applied in the selected regions are proportional to the
-fluid velocity component \f$ u_i \f$.
+The volume force \f$ f_i \f$ applied in the selected regions are proportional to
+the fluid velocity component \f$ u_i \f$.
 
 \f{eqnarray*}{
    f_i(x) &=& - B(x) u_i(x), \\
    B(x) &=& \kappa_0 + (\kappa_1 - \kappa_0) \xi(x) \frac{q + 1}{q + \xi(x)},
  \f}
 
-where, \f$ x \f$ is the current location in the domain, \f$ \xi: x \mapsto [0,1] \f$
-represent an indicator function for the resistance where \f$ \xi(x) = 0 \f$ is a free
-flow. \f$ \kappa_i \f$ describes the limits for the force application at \f$ \xi(x)=0 \f$
-and \f$ \xi(x)=1 \f$. A penalty parameter \f$ q \f$ help us to reduce numerical problems.
+where, \f$ x \f$ is the current location in the domain, \f$ \xi: x \mapsto [0,1]
+\f$ represent an indicator function for the resistance where \f$ \xi(x) = 0 \f$
+is a free flow. \f$ \kappa_i \f$ describes the limits for the force application
+at \f$ \xi(x)=0 \f$ and \f$ \xi(x)=1 \f$. A penalty parameter \f$ q \f$ help us
+to reduce numerical problems.
 
 The indicator function will be defined based on the object type. The following
 types are currently implemented.
@@ -378,19 +379,22 @@ The following keywords are used, with the corresponding options.
   - `bicgstab`, a bi-conjugate gradient stabilized solver.
   - `cacg`, a communication-avoiding conjugate gradient solver.
   - `gmres`, a GMRES solver. Typically used for pressure.
-  - `fusedcg`, a conjugate gradient solver optimised for accelerators using kernel fusion.
+  - `fusedcg`, a conjugate gradient solver optimised for accelerators using
+    kernel fusion.
 * `preconditioner`, preconditioner type.
   - `jacobi`, a Jacobi preconditioner. Typically used for velocity.
-  - `hsmg`, a hybrid-Schwarz multigrid preconditioner. Typically used for pressure.
+  - `hsmg`, a hybrid-Schwarz multigrid preconditioner. Typically used for
+    pressure.
   - `ident`, an identity matrix (no preconditioner).
 * `absolute_tolerance`, tolerance criterion for convergence.
 * `max_iterations`, maximum number of iterations before giving up.
 * `projection_space_size`, size of the vector space used for accelerating the
    solution procedure. If 0, then the projection space is not used.
    More important for the pressure equation.
-* `projection_hold_steps`, steps for which the simulation does not use projection after starting
-   or time step changes. E.g. if 5, then the projection space will start to update at the 6th
-   time step and the space will be utilized at the 7th time step.
+* `projection_hold_steps`, steps for which the simulation does not use
+   projection after starting or time step changes. E.g. if 5, then the
+   projection space will start to update at the 6th time step and the space will
+   be utilized at the 7th time step.
 
 ### Flow rate forcing
 The optional `flow_rate_force` object can be used to force a particular flow

examples/immersed_bunny/immersed_bunny_implicit.case

@@ -0,0 +1,96 @@
+{
+    "version": 1.0,
+    "case": {
+        "mesh_file": "box.nmsh",
+        "output_directory": "fields",
+        "output_at_end": true,
+        "output_boundary": true,
+        "output_checkpoints": false,
+        "end_time": 20.0,
+        "timestep": 1e-4,
+        "numerics": {
+            "time_order": 3,
+            "polynomial_order": 8,
+            "dealias": true
+        },
+        "fluid": {
+            "scheme": "pnpn",
+            "Re": 1000.0,
+            "initial_condition": {
+                "type": "uniform",
+                "value": [
+                    1.0,
+                    0.0,
+                    0.0
+                ]
+            },
+            "inflow_condition": {
+                "type": "uniform",
+                "value": [
+                    1.0
+                ]
+            },
+            "velocity_solver": {
+                "type": "cg",
+                "preconditioner": "jacobi",
+                "projection_space_size": 0,
+                "absolute_tolerance": 1e-8,
+                "max_iterations": 800
+            },
+            "pressure_solver": {
+                "type": "gmres",
+                "preconditioner": "hsmg",
+                "projection_space_size": 0,
+                "absolute_tolerance": 1e-7,
+                "max_iterations": 800
+            },
+            "output_control": "nsamples",
+            "output_value": 250,
+            "boundary_types": [
+                "v",
+                "o",
+                "w",
+                "w",
+                "w",
+                "w"
+            ],
+            "source_terms": [
+                {
+                    "type": "brinkman",
+                    "objects": [
+                        {
+                            "type": "boundary_mesh",
+                            "name": "bunny.stl",
+                            "distance_transform": {
+                                "type": "smooth_step",
+                                "value": 0.05
+                            },
+                            "mesh_transform": {
+                                "type": "bounding_box",
+                                "box_min": [
+                                    0.75,
+                                    0.25,
+                                    0.00
+                                ],
+                                "box_max": [
+                                    1.25,
+                                    0.75,
+                                    0.50
+                                ],
+                                "keep_aspect_ratio": true
+                            }
+                        }
+                    ],
+                    "brinkman": {
+                        "limits": [
+                            0.0,
+                            100000000.0
+                        ],
+                        "penalty": 1.0,
+                        "implicit": true,
+                    },
+                }
+            ]
+        }
+    }
+}

src/fluid/bcknd/cpu/pnpn_res_cpu.f90

@@ -27,7 +27,7 @@ module pnpn_res_cpu
 contains
 
   subroutine pnpn_prs_res_cpu_compute(p, p_res, u, v, w, u_e, v_e, w_e, f_x, &
-       f_y, f_z, c_Xh, gs_Xh, bc_prs_surface,bc_sym_surface, Ax, bd, dt, mu, rho)
+                                      f_y, f_z, c_Xh, gs_Xh, bc_prs_surface,bc_sym_surface, Ax, bd, dt, mu, rho)
     type(field_t), intent(inout) :: p, u, v, w
     type(field_t), intent(inout) :: u_e, v_e, w_e
     type(field_t), intent(inout) :: p_res
@@ -69,11 +69,11 @@ subroutine pnpn_prs_res_cpu_compute(p, p_res, u, v, w, u_e, v_e, w_e, f_x, &
 
     do concurrent (i = 1:n)
        ta1%x(i,1,1,1) = f_x%x(i,1,1,1) / rho &
-            - ((wa1%x(i,1,1,1) * (mu / rho)) * c_Xh%B(i,1,1,1))
+         - ((wa1%x(i,1,1,1) * (mu / rho)) * c_Xh%B(i,1,1,1))
        ta2%x(i,1,1,1) = f_y%x(i,1,1,1) / rho &
-            - ((wa2%x(i,1,1,1) * (mu / rho)) * c_Xh%B(i,1,1,1))
+         - ((wa2%x(i,1,1,1) * (mu / rho)) * c_Xh%B(i,1,1,1))
        ta3%x(i,1,1,1) = f_z%x(i,1,1,1) / rho &
-            - ((wa3%x(i,1,1,1) * (mu / rho)) * c_Xh%B(i,1,1,1))
+         - ((wa3%x(i,1,1,1) * (mu / rho)) * c_Xh%B(i,1,1,1))
     end do
 
     call gs_Xh%op(ta1, GS_OP_ADD)
@@ -94,7 +94,7 @@ subroutine pnpn_prs_res_cpu_compute(p, p_res, u, v, w, u_e, v_e, w_e, f_x, &
 
     do concurrent (i = 1:n)
        p_res%x(i,1,1,1) = (-p_res%x(i,1,1,1)) &
-                        + wa1%x(i,1,1,1) + wa2%x(i,1,1,1) + wa3%x(i,1,1,1)
+         + wa1%x(i,1,1,1) + wa2%x(i,1,1,1) + wa3%x(i,1,1,1)
     end do
 
     !
@@ -119,22 +119,23 @@ subroutine pnpn_prs_res_cpu_compute(p, p_res, u, v, w, u_e, v_e, w_e, f_x, &
 
     do concurrent (i = 1:n)
        p_res%x(i,1,1,1) = p_res%x(i,1,1,1) &
-            - (dtbd * (ta1%x(i,1,1,1) + ta2%x(i,1,1,1) + ta3%x(i,1,1,1)))&
-            - (wa1%x(i,1,1,1) + wa2%x(i,1,1,1) + wa3%x(i,1,1,1))
+         - (dtbd * (ta1%x(i,1,1,1) + ta2%x(i,1,1,1) + ta3%x(i,1,1,1)))&
+         - (wa1%x(i,1,1,1) + wa2%x(i,1,1,1) + wa3%x(i,1,1,1))
     end do
 
     call neko_scratch_registry%relinquish_field(temp_indices)
 
   end subroutine pnpn_prs_res_cpu_compute
 
   subroutine pnpn_vel_res_cpu_compute(Ax, u, v, w, u_res, v_res, w_res, &
-       p, f_x, f_y, f_z, c_Xh, msh, Xh, mu, rho, bd, dt, n)
+                                      p, f_x, f_y, f_z, chi, c_Xh, msh, Xh, &
+                                      mu, rho, bd, dt, n)
     class(ax_t), intent(in) :: Ax
     type(mesh_t), intent(inout) :: msh
     type(space_t), intent(inout) :: Xh
     type(field_t), intent(inout) :: p, u, v, w
     type(field_t), intent(inout) :: u_res, v_res, w_res
-    type(field_t), intent(inout) :: f_x, f_y, f_z
+    type(field_t), intent(inout) :: f_x, f_y, f_z, chi
     type(coef_t), intent(inout) :: c_Xh
     real(kind=rp), intent(in) :: mu
     real(kind=rp), intent(in) :: rho
@@ -147,7 +148,7 @@ subroutine pnpn_vel_res_cpu_compute(Ax, u, v, w, u_res, v_res, w_res, &
 
     do i = 1, n
        c_Xh%h1(i,1,1,1) = mu
-       c_Xh%h2(i,1,1,1) = rho * (bd / dt)
+       c_Xh%h2(i,1,1,1) = rho * (bd / dt) + rho * chi%x(i,1,1,1)
     end do
     c_Xh%ifh2 = .true.