[1592] | 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 |
---|