Numerical advection of shape matrix using Lie integrator

src/2d/parcels/parcel_bc.f90

@@ -23,7 +23,7 @@ end subroutine apply_periodic_bc
         ! @param[inout] position vector of parcel
         ! @param[inout] B matrix of parcel
         subroutine apply_reflective_bc(position, B)
-            double precision, intent(inout) :: position(2), B(2)
+            double precision, intent(inout) :: position(2), B(3)
 
             if (position(2) > upper(2)) then
                 position(2) = two * upper(2) - position(2)

src/2d/parcels/parcel_container.f90

@@ -12,7 +12,7 @@ module parcel_container
     type parcel_container_type
         double precision, allocatable, dimension(:, :) :: &
             position,   &
-            B               ! B matrix entries; ordering B(:, 1) = B11, B(:, 2) = B12
+            B               ! B matrix entries; ordering B(1, :) = B11, B(2, :) = B12, B(3, :) = B33
 
         double precision, allocatable, dimension(:) :: &
             volume,     &
@@ -88,7 +88,7 @@ subroutine parcel_alloc(num)
 
             allocate(parcels%position(2, num))
             allocate(parcels%vorticity(num))
-            allocate(parcels%B(2, num))
+            allocate(parcels%B(3, num))
             allocate(parcels%volume(num))
             allocate(parcels%buoyancy(num))
 #ifndef ENABLE_DRY_MODE

src/2d/parcels/parcel_diagnostics.f90

@@ -88,7 +88,7 @@ subroutine calculate_parcel_diagnostics(velocity)
             double precision :: velocity(:, :)
             integer          :: n
             double precision :: b, z, vel(2), vol, zmin
-            double precision :: eval, lam, B22, lsum, l2sum, vsum, v2sum
+            double precision :: eval, lam, lsum, l2sum, vsum, v2sum
 
             call start_timer(parcel_stats_timer)
 
@@ -121,7 +121,7 @@ subroutine calculate_parcel_diagnostics(velocity)
             std_vol = zero
 
             !$omp parallel default(shared)
-            !$omp do private(n, vel, vol, b, z, eval, lam, B22) &
+            !$omp do private(n, vel, vol, b, z, eval, lam) &
             !$omp& reduction(+: ke, ape, lsum, l2sum, vsum, v2sum, n_small, rms_zeta) &
             !$omp& reduction(-: pe)
             do n = 1, n_parcels
@@ -141,8 +141,7 @@ subroutine calculate_parcel_diagnostics(velocity)
                     ape = ape + ape_den(b, z) * vol
                 endif
 
-                B22 = get_B22(parcels%B(1, n), parcels%B(2, n), vol)
-                eval = get_eigenvalue(parcels%B(1, n), parcels%B(2, n), B22)
+                eval = get_eigenvalue(parcels%B(1, n), parcels%B(2, n), parcels%B(3, n))
                 lam = get_aspect_ratio(eval, vol)
 
                 lsum = lsum + lam

src/2d/parcels/parcel_ellipse.f90

@@ -55,17 +55,14 @@ function get_eigenvector(a2, B11, B12, B22) result(evec)
         end function get_eigenvector
 
         ! used in unit tests only
-        function get_angle(B11, B12, volume) result(angle)
+        function get_angle(B11, B12, B22) result(angle)
             double precision, intent(in) :: B11
             double precision, intent(in) :: B12
-            double precision, intent(in) :: volume
-            double precision             :: B22
+            double precision, intent(in) :: B22
             double precision             :: a2
             double precision             :: evec(2)
             double precision             :: angle
 
-            B22 = get_B22(B11, B12, volume)
-
             a2 = get_eigenvalue(B11, B12, B22)
 
             evec = get_eigenvector(a2, B11, B12, B22)
@@ -119,24 +116,20 @@ elemental function get_aspect_ratio(a2, volume) result(lam)
         ! @param[in] volume of the parcel
         ! @param[in] B matrix elements of the parcel
         ! @returns the parcel support points
-        function get_ellipse_points(position, volume, B) result(points)
+        function get_ellipse_points(position, B) result(points)
             double precision, intent(in) :: position(2)
-            double precision, intent(in) :: volume
-            double precision, intent(in) :: B(2)        ! B11, B12
-            double precision             :: B22
+            double precision, intent(in) :: B(3)        ! B11, B12, B22
             double precision             :: c
             double precision             :: a2
             double precision             :: evec(2)
             double precision             :: h(2)
             double precision             :: points(2, 2)
 
-            B22 = get_B22(B(1), B(2), volume)
-
-            a2 = get_eigenvalue(B(1), B(2), B22)
+            a2 = get_eigenvalue(B(1), B(2), B(3))
 
-            c = dsqrt(dabs(two * a2 - B(1) - B22))
+            c = dsqrt(dabs(two * a2 - B(1) - B(3)))
 
-            evec = get_eigenvector(a2, B(1), B(2), B22)
+            evec = get_eigenvector(a2, B(1), B(2), B(3))
 
             h = f12 * c * evec
 

src/2d/parcels/parcel_init.f90

@@ -5,7 +5,7 @@ module parcel_init
     use options, only : parcel, output, verbose, field_tol
     use constants, only : zero, two, one, f12
     use parcel_container, only : parcels, n_parcels
-    use parcel_ellipse, only : get_ab, get_B22, get_eigenvalue
+    use parcel_ellipse, only : get_ab, get_eigenvalue, get_B22
     use parcel_split, only : split_ellipses
     use parcel_interpl, only : bilinear, ngp
     use parameters, only : dx, vcell, ncell,        &
@@ -83,6 +83,9 @@ subroutine init_parcels(fname, tol)
 
                 ! B12
                 parcels%B(2, n) = zero
+
+                ! B22
+                parcels%B(3, n) = get_B22(parcels%B(1, n), zero, parcels%volume(1))
             enddo
             !$omp end do
             !$omp end parallel
@@ -145,13 +148,12 @@ end subroutine init_regular_positions
 
         subroutine init_refine(lam)
             double precision, intent(inout) :: lam
-            double precision                :: B22, a2
+            double precision                :: a2
 
             ! do refining by splitting
             do while (lam >= parcel%lambda_max)
                 call split_ellipses(parcels, parcel%lambda_max)
-                B22 = get_B22(parcels%B(1, 1), zero, parcels%volume(1))
-                a2 = get_eigenvalue(parcels%B(1, 1), zero, B22)
+                a2 = get_eigenvalue(parcels%B(1, 1), parcels%B(2, 1), parcels%B(3, 1))
                 lam = a2 / get_ab(parcels%volume(1))
             end do
         end subroutine init_refine

src/2d/parcels/parcel_interpl.f90

@@ -50,7 +50,7 @@ subroutine vol2grid
                 pvol = parcels%volume(n)
 
                 points = get_ellipse_points(parcels%position(:, n), &
-                                            pvol, parcels%B(:, n))
+                                            parcels%B(:, n))
 
 
                 ! we have 2 points per ellipse
@@ -85,7 +85,7 @@ end subroutine vol2grid
 #ifndef NDEBUG
         ! Interpolate the parcel volume to the grid to check symmetry
         subroutine vol2grid_symmetry_error
-            double precision :: points(2, 2), V, B(2), pos(2)
+            double precision :: points(2, 2), V, B(3), pos(2)
             integer          :: n, p, l, m
             double precision :: pvol
 
@@ -107,7 +107,7 @@ subroutine vol2grid_symmetry_error
 
                     B(2) = dble(m) * B(2)
 
-                    points = get_ellipse_points(pos, V, B)
+                    points = get_ellipse_points(pos, B)
 
                     ! we have 2 points per ellipse
                     do p = 1, 2
@@ -180,7 +180,7 @@ subroutine par2grid
                 btot = parcels%buoyancy(n)
 #endif
                 points = get_ellipse_points(parcels%position(:, n), &
-                                            pvol, parcels%B(:, n))
+                                            parcels%B(:, n))
 
                 call get_index(parcels%position(:, n), i, j)
                 i = mod(i + nx, nx)
@@ -312,6 +312,7 @@ subroutine grid2par(vel, vor, vgrad, add)
                     do n = 1, n_parcels
                         vel(:, n) = zero
                         vor(n)    = zero
+                        vgrad(:, n) = zero
                     enddo
                     !$omp end do
                     !$omp end parallel
@@ -322,6 +323,7 @@ subroutine grid2par(vel, vor, vgrad, add)
                 do n = 1, n_parcels
                     vel(:, n) = zero
                     vor(n)    = zero
+                    vgrad(:, n) = zero
                 enddo
                 !$omp end do
                 !$omp end parallel
@@ -331,10 +333,7 @@ subroutine grid2par(vel, vor, vgrad, add)
             !$omp do private(n, p, l, points, weight, is, js, weights)
             do n = 1, n_parcels
 
-                vgrad(:, n) = zero
-
                 points = get_ellipse_points(parcels%position(:, n), &
-                                            parcels%volume(n),      &
                                             parcels%B(:, n))
 
                 ! we have 2 points per ellipse

src/2d/parcels/parcel_merge.f90

@@ -9,7 +9,7 @@ module parcel_merge
                                , n_parcels              &
                                , parcel_replace         &
                                , get_delx
-    use parcel_ellipse, only : get_B22, get_ab
+    use parcel_ellipse, only : get_ab
     use options, only : parcel, verbose
     use parcel_bc
     use timer, only : start_timer, stop_timer
@@ -80,7 +80,7 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm)
             integer                                    :: m, ic, is, l, n
             integer                                    :: loca(n_parcels)
             double precision                           :: x0(n_merge)
-            double precision                           :: posm(2, n_merge), delx, vmerge, dely, B22, mu
+            double precision                           :: posm(2, n_merge), delx, vmerge, dely, mu
             double precision                           :: buoym(n_merge), vortm(n_merge)
 #ifndef ENABLE_DRY_MODE
             double precision                           :: hum(n_merge)
@@ -186,15 +186,13 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm)
 
                     vmerge = one / vm(l)
 
-                    B22 = get_B22(parcels%B(1, ic), parcels%B(2, ic), parcels%volume(ic))
-
                     delx = get_delx(parcels%position(1, ic), posm(1, l))
                     dely = parcels%position(2, ic) - posm(2, l)
 
                     mu = parcels%volume(ic) * vmerge
                     B11m(l) = mu * (four * delx ** 2 + parcels%B(1, ic))
                     B12m(l) = mu * (four * delx * dely + parcels%B(2, ic))
-                    B22m(l) = mu * (four * dely ** 2 + B22)
+                    B22m(l) = mu * (four * dely ** 2 + parcels%B(3, ic))
 
                     parcels%volume(ic)  = vm(l)
                     parcels%position(1, ic) = posm(1, l)
@@ -216,14 +214,12 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm)
                 delx = get_delx(parcels%position(1, is), posm(1, n))
                 dely = parcels%position(2, is) - posm(2, n)
 
-                B22 = get_B22(parcels%B(1, is), parcels%B(2, is), parcels%volume(is))
-
                 ! volume fraction A_{is} / A
                 mu = vmerge * parcels%volume(is)
 
                 B11m(n) = B11m(n) + mu * (four * delx ** 2   + parcels%B(1, is))
                 B12m(n) = B12m(n) + mu * (four * delx * dely + parcels%B(2, is))
-                B22m(n) = B22m(n) + mu * (four * dely ** 2   + B22)
+                B22m(n) = B22m(n) + mu * (four * dely ** 2   + parcels%B(3, is))
             enddo
 
         end subroutine do_group_merge
@@ -266,6 +262,7 @@ subroutine geometric_merge(parcels, isma, iclo, n_merge)
 
                     parcels%B(1, ic) = B11(l) * factor
                     parcels%B(2, ic) = B12(l) * factor
+                    parcels%B(3, ic) = B22(l) * factor
 
                     call apply_periodic_bc(parcels%position(:, ic))
                 endif

src/2d/parcels/parcel_netcdf.f90

@@ -21,7 +21,7 @@ module parcel_netcdf
     integer            :: ncid
     integer            :: npar_dim_id, vol_id, buo_id,  &
                           x_pos_id, z_pos_id, vor_id,   &
-                          b11_id, b12_id,               &
+                          b11_id, b12_id, b22_id,       &
                           t_axis_id, t_dim_id
     double precision   :: restart_time
 
@@ -133,6 +133,15 @@ subroutine create_netcdf_parcel_file(basename, overwrite, l_restart)
                                        dimids=dimids,                           &
                                        varid=b12_id)
 
+            call define_netcdf_dataset(ncid=ncid,                               &
+                                       name='B22',                              &
+                                       long_name='B22 element of shape matrix', &
+                                       std_name='',                             &
+                                       unit='m^2',                              &
+                                       dtype=NF90_DOUBLE,                       &
+                                       dimids=dimids,                           &
+                                       varid=b22_id)
+
             call define_netcdf_dataset(ncid=ncid,                               &
                                        name='volume',                           &
                                        long_name='parcel volume',               &
@@ -205,6 +214,7 @@ subroutine write_netcdf_parcels(t)
 
             call write_netcdf_dataset(ncid, b11_id, parcels%B(1, 1:n_parcels), start, cnt)
             call write_netcdf_dataset(ncid, b12_id, parcels%B(2, 1:n_parcels), start, cnt)
+            call write_netcdf_dataset(ncid, b22_id, parcels%B(3, 1:n_parcels), start, cnt)
 
             call write_netcdf_dataset(ncid, vol_id, parcels%volume(1:n_parcels), start, cnt)
 
@@ -263,6 +273,13 @@ subroutine read_netcdf_parcels(fname)
                 stop
             endif
 
+            if (has_dataset(ncid, 'B22')) then
+                call read_netcdf_dataset(ncid, 'B22', parcels%B(3, 1:n_parcels), start, cnt)
+            else
+                print *, "The parcel shape component B22 must be present! Exiting."
+                stop
+            endif
+
             if (has_dataset(ncid, 'x_position')) then
                 call read_netcdf_dataset(ncid, 'x_position', &
                                          parcels%position(1, 1:n_parcels), start, cnt)

src/2d/parcels/parcel_split.f90

@@ -45,7 +45,7 @@ subroutine split_ellipses(parcels, threshold)
                 B11 = parcels%B(1, n)
                 B12 = parcels%B(2, n)
                 V = parcels%volume(n)
-                B22 = get_B22(B11, B12, V)
+                B22 = parcels%B(3, n)
 
                 a2 = get_eigenvalue(B11, B12, B22)
 
@@ -64,6 +64,7 @@ subroutine split_ellipses(parcels, threshold)
 
                 parcels%B(1, n) = B11 - f34 * a2 * evec(1) ** 2
                 parcels%B(2, n) = B12 - f34 * a2 * (evec(1) * evec(2))
+                parcels%B(3, n) = B22 - f34 * a2 * evec(2) ** 2
 
                 h = f14 * dsqrt(three * a2)
                 parcels%volume(n) = f12 * V