Skip to content

Commit

Permalink
Revert "Revert to unscaled ellipses/ellipsoids"
Browse files Browse the repository at this point in the history
  • Loading branch information
@matt-frey
matt-frey authored Apr 19, 2021
1 parent bf5954e commit fe11c23
Show file tree
Hide file tree
Showing 11 changed files with 105 additions and 157 deletions.
16 changes: 8 additions & 8 deletions python-scripts/tools/h5_reader.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,10 +76,11 @@ def get_ellipses(self, step):
V = self.get_parcel_dataset(step, 'volume')
angle = self.get_parcel_dataset(step, 'orientation')

B22 = self._get_B22(B[0, :], B[1, :], V)
a2 = self._get_eigenvalue(B[0, :], B[1, :], B22)
B22 = self._get_B22(B[0, :], B[1, :])
lam = self._get_eigenvalue(B[0, :], B[1, :], B22)

b2 = (V / np.pi) ** 2 / a2
a2 = lam * V / np.pi
b2 = V / np.pi / lam

return [Ellipse(xy=position[:, i],
width=2 * np.sqrt(a2[i]),
Expand All @@ -92,13 +93,12 @@ def get_aspect_ratio(self, step):
B = self.get_parcel_dataset(step, 'B')
V = self.get_parcel_dataset(step, 'volume')

B22 = self._get_B22(B[0, :], B[1, :], V)
a2 = self._get_eigenvalue(B[0, :], B[1, :], B22)
return a2 / V * np.pi
B22 = self._get_B22(B[0, :], B[1, :])
return self._get_eigenvalue(B[0, :], B[1, :], B22)


def _get_B22(self, B11, B12, V):
return ((V / np.pi) ** 2 + B12 ** 2) / B11
def _get_B22(self, B11, B12):
return (1.0 + B12 ** 2) / B11


def _get_eigenvalue(self, B11, B12, B22):
Expand Down
12 changes: 3 additions & 9 deletions src/constants.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,20 +5,14 @@
module constants
implicit none

double precision, parameter :: zero = 0.0d0
double precision, parameter :: zero = 0.d0

double precision, parameter :: one = 1.0d0
double precision, parameter :: one = 1.d0

double precision, parameter :: two = 2.0d0

double precision, parameter :: three = 3.0d0

double precision, parameter :: four = 4.0d0
double precision, parameter :: two = 2.d0

double precision, parameter :: pi = acos(-one)

double precision, parameter :: one_over_pi = one / pi

! maximum number of allowed parcels
integer, parameter :: max_num_parcels = 1e6

Expand Down
11 changes: 4 additions & 7 deletions src/init.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@ module init
volume_f, &
vorticity_f, &
get_position
use ellipse, only : get_ab
use parcel_container, only : parcels, n_parcels
implicit none

Expand Down Expand Up @@ -39,13 +38,12 @@ subroutine init_parcels

call init_stretch

! initialize the volume of each parcel
parcels%volume(1:n_parcels, 1) = vcell / parcel_info%n_per_cell

call init_B_matrix

call init_velocity

! initialize the volume of each parcel
parcels%volume(1:n_parcels, 1) = vcell / parcel_info%n_per_cell

end subroutine init_parcels

Expand Down Expand Up @@ -108,9 +106,8 @@ end subroutine init_stretch

subroutine init_B_matrix
if (parcel_info%is_elliptic) then
! initialze circles
parcels%B(1:n_parcels, 1) = get_ab(parcels%volume(1:n_parcels, 1)) ! B11
parcels%B(1:n_parcels, 2) = 0.0 ! B12
parcels%B(1:n_parcels, 1) = 1.0 ! B11
parcels%B(1:n_parcels, 2) = 0.0 ! B12
else
deallocate(parcels%B)
endif
Expand Down
71 changes: 24 additions & 47 deletions src/parcels/ellipse.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,121 +2,98 @@
! This module contains all ellipse operations.
! =============================================================================
module ellipse
use constants, only : pi &
, one_over_pi &
, zero &
, two
use constants, only : pi
implicit none

contains

! the eigenvalue is the semi-major axis squared (a**2)
function get_eigenvalue(B11, B12, B22) result(a2)
function get_eigenvalue(B11, B12, B22) result(eval)
double precision, intent(in) :: B11
double precision, intent(in) :: B12
double precision, intent(in) :: B22
double precision :: a2
double precision :: eval

a2 = 0.5d0 * (B11 + B22) + sqrt(0.25d0 * (B11 - B22) ** 2 + B12 ** 2)
eval = 0.5 * (B11 + B22) + sqrt(0.25 * (B11 - B22) ** 2 + B12 ** 2)
end function get_eigenvalue

function get_eigenvector(a2, B12, B22) result(evec)
double precision, intent(in) :: a2
function get_eigenvector(eval, B12, B22) result(evec)
double precision, intent(in) :: eval
double precision, intent(in) :: B12
double precision, intent(in) :: B22
double precision :: evec(2)

evec(1) = a2 - B22
evec(1) = eval - B22
evec(2) = B12

if (abs(evec(1)) + abs(evec(2)) == zero) then
if (abs(evec(1)) + abs(evec(2)) == 0.0d0) then
evec = evec + epsilon(evec)
endif

evec = evec / norm2(evec)

end function get_eigenvector

function get_angle(B11, B12, volume) result(angle)
function get_angle(B11, B12) result(angle)
double precision, intent(in) :: B11
double precision, intent(in) :: B12
double precision, intent(in) :: volume
double precision :: B22
double precision :: a2
double precision :: eval
double precision :: evec(2)
double precision :: angle

B22 = get_B22(B11, B12, volume)
B22 = get_B22(B11, B12)

a2 = get_eigenvalue(B11, B12, B22)
eval = get_eigenvalue(B11, B12, B22)

evec = get_eigenvector(a2, B12, B22)
evec = get_eigenvector(eval, B12, B22)

angle = atan2(evec(2), evec(1))

end function get_angle

function get_angles(B, volume, n_parcels) result(angle)
function get_angles(B, n_parcels) result(angle)
double precision, intent(in) :: B(:, :)
double precision, intent(in) :: volume(:, :)
integer, intent(in) :: n_parcels
double precision :: angle(n_parcels)
integer :: n

do n = 1, n_parcels
angle(n) = get_angle(B(n, 1), B(n, 2), volume(n, 1))
angle(n) = get_angle(B(n, 1), B(n, 2))
enddo
end function get_angles


elemental function get_B22(B11, B12, volume) result(B22)
elemental function get_B22(B11, B12) result(B22)
double precision, intent(in) :: B11
double precision, intent(in) :: B12
double precision, intent(in) :: volume
double precision :: B22
double precision :: ab

ab = get_ab(volume)

B22 = (ab ** 2 + B12 ** 2) / B11
B22 = (1.0 + B12 ** 2) / B11

end function get_B22

elemental function get_ab(volume) result(ab)
double precision, intent(in) :: volume
double precision :: ab

ab = volume * one_over_pi
end function

elemental function get_aspect_ratio(a2, volume) result(lam)
double precision, intent(in) :: a2
double precision, intent(in) :: volume
double precision :: lam

lam = (a2 / volume) * pi
end function

function get_ellipse_points(position, volume, 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 :: c
double precision :: a2
double precision :: eval
double precision :: evec(2)
double precision :: dx(2)
double precision :: points(2, 2)

B22 = get_B22(B(1), B(2), volume)
B22 = get_B22(B(1), B(2))

a2 = get_eigenvalue(B(1), B(2), B22)
eval = get_eigenvalue(B(1), B(2), B22)

c = sqrt(abs(two * a2 - B(1) - B22))
! a * b = volume / pi with a and b semi-major and semi-minor axes
c = sqrt(volume / pi * abs(2.0 * eval - B(1) - B22))

evec = get_eigenvector(a2, B(2), B22)
evec = get_eigenvector(eval, B(2), B22)

dx = 0.5d0 * c * evec
dx = 0.5 * c * evec

points(1, :) = position - dx
points(2, :) = position + dx
Expand Down
2 changes: 1 addition & 1 deletion src/parcels/parcel_container.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ subroutine write_h5_parcels(iter)
if (allocated(parcels%B)) then
call write_h5_dataset_2d(name, "B", parcels%B(1:n_parcels, :))

angle = get_angles(parcels%B, parcels%volume, n_parcels)
angle = get_angles(parcels%B, n_parcels)
call write_h5_dataset_1d(name, "orientation", angle)
endif

Expand Down
Loading

0 comments on commit fe11c23

Please sign in to comment.