Add a curl operator.

src/solver.f90

@@ -47,6 +47,7 @@ module m_solver
       procedure :: divergence_v2p
       procedure :: divergence_v2v
       procedure :: gradient_p2v
+      procedure :: curl_v2v
       procedure :: run
    end type solver_t
 
@@ -442,6 +443,67 @@ subroutine gradient_p2v(self, dpdx, dpdy, dpdz, pressure)
 
    end subroutine gradient_p2v
 
+   subroutine curl_v2v(self, o_x, o_y, o_z, u, v, w)
+      implicit none
+
+      class(solver_t) :: self
+      class(field_t), intent(inout) :: o_x, o_y, o_z !! omega_x/_y/_z
+      class(field_t), intent(in) :: u, v, w
+
+      class(field_t), pointer :: dv_x, dw_x, du_y, dw_y, du_z, dv_z
+      class(field_t), pointer :: u_y, w_y, u_z, v_z
+
+      ! o_x = dw/dy - dv/dz
+      ! o_y = du/dz - dw/dx
+      ! o_z = dv/dx - du/dy
+
+      ! obtain dw/dx, dv/dx and store them directly in omega_y, omega_z
+      call self%backend%tds_solve(o_y, w, self%xdirps, self%xdirps%der1st)
+      call self%backend%tds_solve(o_z, v, self%xdirps, self%xdirps%der1st)
+
+      u_y => self%backend%allocator%get_block()
+      w_y => self%backend%allocator%get_block()
+
+      call self%backend%reorder(u_y, u, RDR_X2Y)
+      call self%backend%reorder(w_y, w, RDR_X2Y)
+
+      du_y => self%backend%allocator%get_block()
+
+      ! obtain du/dy, dw/dy
+      ! store du/dy in a temporary field to add into omega_z later
+      ! dw/dy can be stored directly in omega_x as it is empty
+      call self%backend%tds_solve(du_y, u_y, self%ydirps, self%ydirps%der1st)
+      call self%backend%tds_solve(o_x, w_y, self%ydirps, self%ydirps%der1st)
+
+      call self%backend%allocator%release_block(u_y)
+      call self%backend%allocator%release_block(w_y)
+
+      ! omega_z = dv/dz - du/dy
+      call self%backend%vecadd(-1._dp, du_y, 1._dp, o_z)
+
+      call self%backend%allocator%release_block(du_y)
+
+      u_z => self%backend%allocator%get_block()
+      v_z => self%backend%allocator%get_block()
+      du_z => self%backend%allocator%get_block()
+      dv_z => self%backend%allocator%get_block()
+
+      ! obtain du/dz, dv/dz and store them in temporary fields
+      call self%backend%tds_solve(du_z, u_z, self%zdirps, self%zdirps%der1st)
+      call self%backend%tds_solve(dv_z, v_z, self%zdirps, self%zdirps%der1st)
+
+      ! omega_x = dw/dy - dv/dz
+      call self%backend%vecadd(-1._dp, dv_z, 1._dp, o_x)
+      ! omega_y = du/dz - dw/dx
+      call self%backend%vecadd(1._dp, du_z, -1._dp, o_y)
+
+      call self%backend%allocator%release_block(u_z)
+      call self%backend%allocator%release_block(v_z)
+      call self%backend%allocator%release_block(du_z)
+      call self%backend%allocator%release_block(dv_z)
+
+   end subroutine curl_v2v
+
    subroutine run(self, n_iter, u_out, v_out, w_out)
       implicit none