-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathutils.f90
156 lines (130 loc) · 4.7 KB
/
utils.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
!------------------------------------------------------------------------------!
! NDSPMHD: A Smoothed Particle (Magneto)Hydrodynamics code for (astrophysical) !
! fluid dynamics simulations in 1, 2 and 3 spatial dimensions. !
! !
! (c) 2002-2015 Daniel Price !
! !
! http://users.monash.edu.au/~dprice/ndspmhd !
! daniel.price@monash.edu -or- dprice@cantab.net (forwards to current address) !
! !
! NDSPMHD comes with ABSOLUTELY NO WARRANTY. !
! This is free software; and you are welcome to redistribute !
! it under the terms of the GNU General Public License !
! (see LICENSE file for details) and the provision that !
! this notice remains intact. If you modify this file, please !
! note section 2a) of the GPLv2 states that: !
! !
! a) You must cause the modified files to carry prominent notices !
! stating that you changed the files and the date of any change. !
! !
! ChangeLog: !
!------------------------------------------------------------------------------!
!-------------------------------------------------------------------
! This file contains some simple utilities used in various places
! throughout the code
!------------------------------------------------------------------
module utils
implicit none
public :: minmaxave,cross_product3D,curl3D_epsijk,det,delta_fn
contains
!-------------------------------------------------------------------
! simple routine to take min, max and average of a quantity
!-------------------------------------------------------------------
subroutine minmaxave(x,xmin,xmax,xav,npts)
integer :: i
integer, intent(in) :: npts
real, intent(in), dimension(npts) :: x
real, intent(out) :: xmin,xmax,xav
xav = 0.
xmin = huge(xmin)
xmax = -xmin
do i=1,npts
xav = xav + x(i)
xmin = min(xmin,x(i))
xmax = max(xmax,x(i))
enddo
xav = xav/real(npts)
return
end subroutine minmaxave
!-------------------------------
!+
! cross product of two vectors
!+
!-------------------------------
pure subroutine cross_product3D(veca,vecb,vecc)
real, dimension(3), intent(in) :: veca,vecb
real, dimension(3), intent(out) :: vecc
vecc(1) = veca(2)*vecb(3) - veca(3)*vecb(2)
vecc(2) = veca(3)*vecb(1) - veca(1)*vecb(3)
vecc(3) = veca(1)*vecb(2) - veca(2)*vecb(1)
end subroutine cross_product3D
pure subroutine curl3D_epsijk(gradAvec,curlA)
real, dimension(3,3), intent(in) :: gradAvec
real, dimension(3), intent(out) :: curlA
curlA(1) = gradAvec(2,3) - gradAvec(3,2)
curlA(2) = gradAvec(3,1) - gradAvec(1,3)
curlA(3) = gradAvec(1,2) - gradAvec(2,1)
end subroutine curl3D_epsijk
!-------------------------------
!+
! Kronecker delta
!+
!-------------------------------
pure integer function delta_fn(i,j)
integer, intent(in) :: i,j
if (i==j) then
delta_fn = 1
else
delta_fn = 0
endif
end function delta_fn
!----------------------------------------------------------------
!+
! Internal subroutine that inverts a 3x3 matrix
!+
!----------------------------------------------------------------
subroutine matrixinvert3D(A,Ainv,ierr)
real, intent(in), dimension(3,3) :: A
real, intent(out), dimension(3,3) :: Ainv
integer, intent(out) :: ierr
real, dimension(3) :: x0,x1,x2,result
real :: det, ddet
ierr = 0
x0 = A(1,:)
x1 = A(2,:)
x2 = A(3,:)
call cross_product3D(x1,x2,result)
det = dot_product(x0,result)
if (abs(det).gt.tiny(det)) then
ddet = 1./det
else
print*,' matrix = ',A(:,:)
print*,' determinant = ',det
ddet = 0.
Ainv = 0.
ierr = 1
return
endif
Ainv(:,1) = result(:)*ddet
call cross_product3D(x2,x0,result)
Ainv(:,2) = result(:)*ddet
call cross_product3D(x0,x1,result)
Ainv(:,3) = result(:)*ddet
return
end subroutine matrixinvert3D
!-----------------------------------
!+
! get determinant of 3 x 3 matrix
!+
!-----------------------------------
pure real function det(A)
real, intent(in), dimension(3,3) :: A
real, dimension(3) :: x0,x1,x2,result
x0 = A(1,:)
x1 = A(2,:)
x2 = A(3,:)
call cross_product3D(x1,x2,result)
det = dot_product(x0,result)
return
end function det
end module utils