-
Notifications
You must be signed in to change notification settings - Fork 1
/
jacobians.f90
executable file
·118 lines (81 loc) · 2.12 KB
/
jacobians.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
module jacobians
use global
use write_pack
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function matrix_action(x) result(y)
implicit none
real(dp), intent(in) :: x(:)
real(dp), allocatable :: y(:)
real(dp), allocatable :: x1(:), x2(:), x3(:)
real(dp) :: jj
integer :: n, ii
n = size(x)
allocate(x1(n), x2(n), x3(n), y(n), stat=alloc_err)
x1 = 0.0_dp
x2 = 0.0_dp
x3 = 0.0_dp
y = 0.0_dp
x1(1) = 0.0_dp
x1(2:n) = x(1:n-1)
x2(1) = real((n-1) / 2, kind=dp)
jj = 1.0_dp
do ii = 2,n
if (ii <= (n+1)/2) then
x2(ii) = x2(ii-1) - 1.0_dp
else
x2(ii) = jj
jj = jj + 1
end if
end do
x2 = x2*x
x3 = x(2:n)
y = x1 + x2 + x3
end function matrix_action
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function jac_approx_flow_map(x) result(jact)
use global
use time_integrators
implicit none
real(dp), intent(in) :: x(:)
real(dp), allocatable :: jact(:)
real(dp), allocatable :: x0_delta(:), xT_delta(:), GT_delta(:)
real(dp) :: normx, dot_prod, eps
real(dp) :: dnrm2, ddot
integer :: n, ii
integer, parameter :: incx = 1, incy = 1
n = size(x)
allocate(jact(n), stat=alloc_err)
allocate(x0_delta(n), xT_Delta(n), GT_Delta(n), stat=alloc_err)
call check_alloc_err(alloc_err)
jact = 0.0_dp
x0_delta = 0.0_dp
xT_delta = 0.0_dp
GT_delta = 0.0_dp
normx = dnrm2(n,x,incx)
dot_prod = ddot(n,x,incx,x0,incy)
eps = sqrt(epsilon(1.0_dp)) / (normx**2.0_dp) * &
&max(dot_prod, normx)
x0_delta = x0 + eps * x
! Upack fields
call unpackx(x0_delta)
! Calculate phi and ux from uy
do ii = 1,Nx
if (kx(ii) /= 0.0_dp) then
tmp_uy = uy(:,ii)
ux(:,ii) = -CI*d1y(tmp_uy)/kx(ii)
else if (kx(ii) == 0.0_dp) then
ux(:,ii) = cmplx(0.0_dp, 0.0_dp, kind=C_DOUBLE_COMPLEX) ! Zero mean flow!
end if
phi(:,ii) = -kx(ii)**2.0_dp*uy(:,ii) + d2y(uy(:,ii))
end do
! Integrate out to some time
call imex_rk
! Form xT_delta from the fields
call packx(xT_delta)
! Compute flow map
GT_delta = (xT_delta - x0_delta) / t_final
jact = (GT_delta - GT) / eps
end function jac_approx_flow_map
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module jacobians