Context Navigation

slopes_m.F90 @ 5326

Last change on this file since 5326 was 2471, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes r2434:2457 into testing branch
File size: 6.3 KB

Line
1	module slopes_m
2
3	! Author: Lionel GUEZ
4
5	implicit none
6
7	interface slopes
8	! This generic function computes second order slopes with Van
9	! Leer slope-limiting, given cell averages. Reference: Dukowicz,
10	! 1987, SIAM Journal on Scientific and Statistical Computing, 8,
11	! 305.
12
13	! The only difference between the specific functions is the rank
14	! of the first argument and the equal rank of the result.
15
16	! real, intent(in), rank >= 1:: f ! (n, ...) cell averages, n must be >= 1
17	! real, intent(in):: x(:) ! (n + 1) cell edges
18	! real slopes, same shape as f ! (n, ...)
19
20	module procedure slopes1, slopes2, slopes3, slopes4
21	end interface
22
23	private
24	public slopes
25
26	contains
27
28	pure function slopes1(f, x)
29
30	real, intent(in):: f(:)
31	real, intent(in):: x(:)
32	real slopes1(size(f))
33
34	! Local:
35	integer n, i
36	real xc(size(f)) ! (n) cell centers
37	real h
38
39	!------------------------------------------------------
40
41	n = size(f)
42	forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2.
43	slopes1(1) = 0.
44	slopes1(n) = 0.
45
46	do i = 2, n - 1
47	if (f(i) >= max(f(i - 1), f(i + 1)) .or. f(i) &
48	<= min(f(i - 1), f(i + 1))) then
49	! Local extremum
50	slopes1(i) = 0.
51	else
52	! (f(i - 1), f(i), f(i + 1)) strictly monotonous
53
54	! Second order slope:
55	slopes1(i) = (f(i + 1) - f(i - 1)) / (xc(i + 1) - xc(i - 1))
56
57	! Slope limitation:
58	h = abs(x(i + 1) - xc(i))
59	slopes1(i) = sign(min(abs(slopes1(i)), abs(f(i + 1) - f(i)) / h, &
60	abs(f(i) - f(i - 1)) / h), slopes1(i))
61	end if
62	end do
63
64	end function slopes1
65
66	!*************************************************************
67
68	pure function slopes2(f, x)
69
70	real, intent(in):: f(:, :)
71	real, intent(in):: x(:)
72	real slopes2(size(f, 1), size(f, 2))
73
74	! Local:
75	integer n, i, j
76	real xc(size(f, 1)) ! (n) cell centers
77	real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1)
78
79	!------------------------------------------------------
80
81	n = size(f, 1)
82	forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2.
83
84	forall (i = 2:n - 1)
85	h(i) = abs(x(i + 1) - xc(i))
86	delta_xc(i) = xc(i + 1) - xc(i - 1)
87	end forall
88
89	do j = 1, size(f, 2)
90	slopes2(1, j) = 0.
91
92	do i = 2, n - 1
93	if (f(i, j) >= max(f(i - 1, j), f(i + 1, j)) .or. &
94	f(i, j) <= min(f(i - 1, j), f(i + 1, j))) then
95	! Local extremum
96	slopes2(i, j) = 0.
97	else
98	! (f(i - 1, j), f(i, j), f(i + 1, j))
99	! strictly monotonous
100
101	! Second order slope:
102	slopes2(i, j) = (f(i + 1, j) - f(i - 1, j)) / delta_xc(i)
103
104	! Slope limitation:
105	slopes2(i, j) = sign(min(abs(slopes2(i, j)), &
106	abs(f(i + 1, j) - f(i, j)) / h(i), &
107	abs(f(i, j) - f(i - 1, j)) / h(i)), slopes2(i, j))
108	end if
109	end do
110
111	slopes2(n, j) = 0.
112	end do
113
114	end function slopes2
115
116	!*************************************************************
117
118	pure function slopes3(f, x)
119
120	real, intent(in):: f(:, :, :)
121	real, intent(in):: x(:)
122	real slopes3(size(f, 1), size(f, 2), size(f, 3))
123
124	! Local:
125	integer n, i, j, k
126	real xc(size(f, 1)) ! (n) cell centers
127	real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1)
128
129	!------------------------------------------------------
130
131	n = size(f, 1)
132	forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2.
133
134	forall (i = 2:n - 1)
135	h(i) = abs(x(i + 1) - xc(i))
136	delta_xc(i) = xc(i + 1) - xc(i - 1)
137	end forall
138
139	do k = 1, size(f, 3)
140	do j = 1, size(f, 2)
141	slopes3(1, j, k) = 0.
142
143	do i = 2, n - 1
144	if (f(i, j, k) >= max(f(i - 1, j, k), f(i + 1, j, k)) .or. &
145	f(i, j, k) <= min(f(i - 1, j, k), f(i + 1, j, k))) then
146	! Local extremum
147	slopes3(i, j, k) = 0.
148	else
149	! (f(i - 1, j, k), f(i, j, k), f(i + 1, j, k))
150	! strictly monotonous
151
152	! Second order slope:
153	slopes3(i, j, k) = (f(i + 1, j, k) - f(i - 1, j, k)) &
154	/ delta_xc(i)
155
156	! Slope limitation:
157	slopes3(i, j, k) = sign(min(abs(slopes3(i, j, k)), &
158	abs(f(i + 1, j, k) - f(i, j, k)) / h(i), &
159	abs(f(i, j, k) - f(i - 1, j, k)) / h(i)), slopes3(i, j, k))
160	end if
161	end do
162
163	slopes3(n, j, k) = 0.
164	end do
165	end do
166
167	end function slopes3
168
169	!*************************************************************
170
171	pure function slopes4(f, x)
172
173	real, intent(in):: f(:, :, :, :)
174	real, intent(in):: x(:)
175	real slopes4(size(f, 1), size(f, 2), size(f, 3), size(f, 4))
176
177	! Local:
178	integer n, i, j, k, l
179	real xc(size(f, 1)) ! (n) cell centers
180	real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1)
181
182	!------------------------------------------------------
183
184	n = size(f, 1)
185	forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2.
186
187	forall (i = 2:n - 1)
188	h(i) = abs(x(i + 1) - xc(i))
189	delta_xc(i) = xc(i + 1) - xc(i - 1)
190	end forall
191
192	do l = 1, size(f, 4)
193	do k = 1, size(f, 3)
194	do j = 1, size(f, 2)
195	slopes4(1, j, k, l) = 0.
196
197	do i = 2, n - 1
198	if (f(i, j, k, l) >= max(f(i - 1, j, k, l), f(i + 1, j, k, l)) &
199	.or. f(i, j, k, l) &
200	<= min(f(i - 1, j, k, l), f(i + 1, j, k, l))) then
201	! Local extremum
202	slopes4(i, j, k, l) = 0.
203	else
204	! (f(i - 1, j, k, l), f(i, j, k, l), f(i + 1, j, k, l))
205	! strictly monotonous
206
207	! Second order slope:
208	slopes4(i, j, k, l) = (f(i + 1, j, k, l) &
209	- f(i - 1, j, k, l)) / delta_xc(i)
210
211	! Slope limitation:
212	slopes4(i, j, k, l) = sign(min(abs(slopes4(i, j, k, l)), &
213	abs(f(i + 1, j, k, l) - f(i, j, k, l)) / h(i), &
214	abs(f(i, j, k, l) - f(i - 1, j, k, l)) / h(i)), &
215	slopes4(i, j, k, l))
216	end if
217	end do
218
219	slopes4(n, j, k, l) = 0.
220	end do
221	end do
222	end do
223
224	end function slopes4
225
226	end module slopes_m

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format