| 1 | |
|---|
| 2 | ! |
|---|
| 3 | ! This program may be freely redistributed under the |
|---|
| 4 | ! condition that the copyright notices (including this |
|---|
| 5 | ! entire header) are not removed, and no compensation |
|---|
| 6 | ! is received through use of the software. Private, |
|---|
| 7 | ! research, and institutional use is free. You may |
|---|
| 8 | ! distribute modified versions of this code UNDER THE |
|---|
| 9 | ! CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE |
|---|
| 10 | ! TO IT IN THE SAME FILE REMAIN UNDER COPYRIGHT OF THE |
|---|
| 11 | ! ORIGINAL AUTHOR, BOTH SOURCE AND OBJECT CODE ARE |
|---|
| 12 | ! MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR |
|---|
| 13 | ! NOTICE IS GIVEN OF THE MODIFICATIONS. Distribution |
|---|
| 14 | ! of this code as part of a commercial system is |
|---|
| 15 | ! permissible ONLY BY DIRECT ARRANGEMENT WITH THE |
|---|
| 16 | ! AUTHOR. (If you are not directly supplying this |
|---|
| 17 | ! code to a customer, and you are instead telling them |
|---|
| 18 | ! how they can obtain it for free, then you are not |
|---|
| 19 | ! required to make any arrangement with me.) |
|---|
| 20 | ! |
|---|
| 21 | ! Disclaimer: Neither I nor: Columbia University, the |
|---|
| 22 | ! National Aeronautics and Space Administration, nor |
|---|
| 23 | ! the Massachusetts Institute of Technology warrant |
|---|
| 24 | ! or certify this code in any way whatsoever. This |
|---|
| 25 | ! code is provided "as-is" to be used at your own risk. |
|---|
| 26 | ! |
|---|
| 27 | ! |
|---|
| 28 | |
|---|
| 29 | ! |
|---|
| 30 | ! FFSL1D.f90: upwind-biased flux-reconstruction scheme. |
|---|
| 31 | ! |
|---|
| 32 | ! Darren Engwirda |
|---|
| 33 | ! 31-Mar-2019 |
|---|
| 34 | ! de2363 [at] columbia [dot] edu |
|---|
| 35 | ! |
|---|
| 36 | ! |
|---|
| 37 | |
|---|
| 38 | subroutine ffsl1d(npos,nvar,ndof,spac,tDEL, & |
|---|
| 39 | & mask,uvel,qbar,qedg,bclo, & |
|---|
| 40 | & bchi,work,opts) |
|---|
| 41 | |
|---|
| 42 | ! |
|---|
| 43 | ! NPOS no. edges over grid. |
|---|
| 44 | ! NVAR no. state variables. |
|---|
| 45 | ! NDOF no. degrees-of-freedom per grid-cell. |
|---|
| 46 | ! SPAC grid-cell spacing array. LENGTH(SPAC) == +1 if |
|---|
| 47 | ! spacing is uniform . |
|---|
| 48 | ! TDEL time-step . |
|---|
| 49 | ! MASK logical grid-cell masking array. |
|---|
| 50 | ! UVEL edge-centred velocity vectors. UVEL has SIZE = |
|---|
| 51 | ! NPOS-by-1 . |
|---|
| 52 | ! QBAR cell-centred integral moments. QBAR has SIZE = |
|---|
| 53 | ! NDOF-by-NVAR-by-NPOS-1 . |
|---|
| 54 | ! QEDG edge-centred upwind flux eval. QEDG has SIZE = |
|---|
| 55 | ! NVAR-by-NPOS . |
|---|
| 56 | ! BCLO boundary condition at lower endpoint . |
|---|
| 57 | ! BCHI boundary condition at upper endpoint . |
|---|
| 58 | ! WORK method work-space. See RCON-WORK for details . |
|---|
| 59 | ! OPTS method parameters. See RCON-OPTS for details . |
|---|
| 60 | ! |
|---|
| 61 | |
|---|
| 62 | implicit none |
|---|
| 63 | |
|---|
| 64 | !------------------------------------------- arguments ! |
|---|
| 65 | integer, intent(in) :: npos,nvar,ndof |
|---|
| 66 | class(rmap_work), intent(inout):: work |
|---|
| 67 | class(rmap_opts), intent(inout):: opts |
|---|
| 68 | real*8 , intent(in) :: spac(:) |
|---|
| 69 | real*8 , intent(in) :: tDEL |
|---|
| 70 | logical, intent(in) :: mask(:) |
|---|
| 71 | real*8 , intent(in) :: qbar(:,:,:) |
|---|
| 72 | real*8 , intent(in) :: uvel(:) |
|---|
| 73 | real*8 , intent(out) :: qedg(:,:) |
|---|
| 74 | class(rcon_ends), intent(in) :: bclo(:) |
|---|
| 75 | class(rcon_ends), intent(in) :: bchi(:) |
|---|
| 76 | |
|---|
| 77 | !------------------------------------------- variables ! |
|---|
| 78 | integer :: head,tail,nprt |
|---|
| 79 | |
|---|
| 80 | head = +0 ; tail = +0 ; qedg = 0.d+0 |
|---|
| 81 | |
|---|
| 82 | do while (.true.) |
|---|
| 83 | |
|---|
| 84 | !--------------------------------- 1. find active part ! |
|---|
| 85 | |
|---|
| 86 | do head = tail+1, npos-1 |
|---|
| 87 | if (mask(head) .eqv..true.) exit |
|---|
| 88 | end do |
|---|
| 89 | |
|---|
| 90 | do tail = head+1, npos-1 |
|---|
| 91 | if (mask(tail).neqv..true.) exit |
|---|
| 92 | end do |
|---|
| 93 | tail = tail - 1 |
|---|
| 94 | |
|---|
| 95 | if (head.ge.npos) exit |
|---|
| 96 | |
|---|
| 97 | !--------------------------------- 2. rcon active part ! |
|---|
| 98 | |
|---|
| 99 | nprt = tail - head + 1 |
|---|
| 100 | |
|---|
| 101 | if (size(spac).ne.+1) then |
|---|
| 102 | |
|---|
| 103 | call rcon1d(nprt+1,nvar,ndof , & |
|---|
| 104 | & spac( head:tail), & |
|---|
| 105 | & qbar(:,:,head:tail), & |
|---|
| 106 | & bclo,bchi,work%cell_func, & |
|---|
| 107 | & work,opts ) |
|---|
| 108 | |
|---|
| 109 | else |
|---|
| 110 | |
|---|
| 111 | call rcon1d(nprt+1,nvar,ndof , & |
|---|
| 112 | & spac,qbar(:,:,head:tail), & |
|---|
| 113 | & bclo,bchi,work%cell_func, & |
|---|
| 114 | & work,opts ) |
|---|
| 115 | |
|---|
| 116 | end if |
|---|
| 117 | |
|---|
| 118 | !--------------------------------- 3. int. active part ! |
|---|
| 119 | |
|---|
| 120 | select case(opts%cell_meth) |
|---|
| 121 | case(pcm_method) !! 1st-order scheme |
|---|
| 122 | |
|---|
| 123 | if (size(spac).ne.+1) then |
|---|
| 124 | |
|---|
| 125 | call flux1d(nprt+1,nvar,1, & |
|---|
| 126 | & spac( head:tail+0) , & |
|---|
| 127 | & tDEL, & |
|---|
| 128 | & uvel( head:tail+1) , & |
|---|
| 129 | & work%cell_func, & |
|---|
| 130 | & qedg(:,head:tail+1) ) |
|---|
| 131 | |
|---|
| 132 | else |
|---|
| 133 | |
|---|
| 134 | call flux1d(nprt+1,nvar,1, & |
|---|
| 135 | & spac,tDEL , & |
|---|
| 136 | & uvel( head:tail+1) , & |
|---|
| 137 | & work%cell_func, & |
|---|
| 138 | & qedg(:,head:tail+1) ) |
|---|
| 139 | |
|---|
| 140 | end if |
|---|
| 141 | |
|---|
| 142 | case(plm_method) !! 2nd-order scheme |
|---|
| 143 | |
|---|
| 144 | if (size(spac).ne.+1) then |
|---|
| 145 | |
|---|
| 146 | call flux1d(nprt+1,nvar,2, & |
|---|
| 147 | & spac( head:tail+0) , & |
|---|
| 148 | & tDEL, & |
|---|
| 149 | & uvel( head:tail+1) , & |
|---|
| 150 | & work%cell_func, & |
|---|
| 151 | & qedg(:,head:tail+1) ) |
|---|
| 152 | |
|---|
| 153 | else |
|---|
| 154 | |
|---|
| 155 | call flux1d(nprt+1,nvar,2, & |
|---|
| 156 | & spac,tDEL , & |
|---|
| 157 | & uvel( head:tail+1) , & |
|---|
| 158 | & work%cell_func, & |
|---|
| 159 | & qedg(:,head:tail+1) ) |
|---|
| 160 | |
|---|
| 161 | end if |
|---|
| 162 | |
|---|
| 163 | case(ppm_method) !! 3rd-order scheme |
|---|
| 164 | |
|---|
| 165 | if (size(spac).ne.+1) then |
|---|
| 166 | |
|---|
| 167 | call flux1d(nprt+1,nvar,3, & |
|---|
| 168 | & spac( head:tail+0) , & |
|---|
| 169 | & tDEL, & |
|---|
| 170 | & uvel( head:tail+1) , & |
|---|
| 171 | & work%cell_func, & |
|---|
| 172 | & qedg(:,head:tail+1) ) |
|---|
| 173 | |
|---|
| 174 | else |
|---|
| 175 | |
|---|
| 176 | call flux1d(nprt+1,nvar,3, & |
|---|
| 177 | & spac,tDEL , & |
|---|
| 178 | & uvel( head:tail+1) , & |
|---|
| 179 | & work%cell_func, & |
|---|
| 180 | & qedg(:,head:tail+1) ) |
|---|
| 181 | |
|---|
| 182 | end if |
|---|
| 183 | |
|---|
| 184 | case(pqm_method) !! 5th-order scheme |
|---|
| 185 | |
|---|
| 186 | if (size(spac).ne.+1) then |
|---|
| 187 | |
|---|
| 188 | call flux1d(nprt+1,nvar,5, & |
|---|
| 189 | & spac( head:tail+0) , & |
|---|
| 190 | & tDEL, & |
|---|
| 191 | & uvel( head:tail+1) , & |
|---|
| 192 | & work%cell_func, & |
|---|
| 193 | & qedg(:,head:tail+1) ) |
|---|
| 194 | |
|---|
| 195 | else |
|---|
| 196 | |
|---|
| 197 | call flux1d(nprt+1,nvar,5, & |
|---|
| 198 | & spac,tDEL , & |
|---|
| 199 | & uvel( head:tail+1) , & |
|---|
| 200 | & work%cell_func, & |
|---|
| 201 | & qedg(:,head:tail+1) ) |
|---|
| 202 | |
|---|
| 203 | end if |
|---|
| 204 | |
|---|
| 205 | end select |
|---|
| 206 | |
|---|
| 207 | end do |
|---|
| 208 | |
|---|
| 209 | return |
|---|
| 210 | |
|---|
| 211 | end subroutine |
|---|
| 212 | |
|---|
| 213 | ! FLUX1D: a degree-k, upwind-type flux reconstruction. ! |
|---|
| 214 | |
|---|
| 215 | pure subroutine flux1d(npos,nvar,mdof,SPAC, & |
|---|
| 216 | & tDEL,uvel,QHAT,qedg) |
|---|
| 217 | |
|---|
| 218 | ! |
|---|
| 219 | ! NPOS no. edges over grid. |
|---|
| 220 | ! NVAR no. state variables. |
|---|
| 221 | ! MDOF no. degrees-of-freedom per QHAT. |
|---|
| 222 | ! SPAC grid spacing vector. SIZE(SPAC)==+1 if uniform . |
|---|
| 223 | ! TDEL time-step . |
|---|
| 224 | ! UVEL edge-centred velocity vectors. UVEL has SIZE = |
|---|
| 225 | ! NPOS-by-1 . |
|---|
| 226 | ! QHAT cell-centred polynomial recon. QHAT has SIZE = |
|---|
| 227 | ! NDOF-by-NVAR-by-NPOS-1 . |
|---|
| 228 | ! QEDG edge-centred upwind flux eval. QEDG has SIZE = |
|---|
| 229 | ! NVAR-by-NPOS . |
|---|
| 230 | ! |
|---|
| 231 | |
|---|
| 232 | implicit none |
|---|
| 233 | |
|---|
| 234 | !------------------------------------------- arguments ! |
|---|
| 235 | integer, intent(in) :: npos,nvar,mdof |
|---|
| 236 | real*8 , intent(in) :: SPAC(:) |
|---|
| 237 | real*8 , intent(in) :: tDEL |
|---|
| 238 | real*8 , intent(in) :: uvel(:) |
|---|
| 239 | real*8 , intent(in) :: QHAT(:,:,:) |
|---|
| 240 | real*8 , intent(out) :: qedg(:,:) |
|---|
| 241 | |
|---|
| 242 | !------------------------------------------- variables ! |
|---|
| 243 | integer :: ipos,ivar |
|---|
| 244 | real*8 :: uCFL,xhat,ss11,ss22,flux |
|---|
| 245 | real*8 :: vv11(1:5) |
|---|
| 246 | real*8 :: vv22(1:5) |
|---|
| 247 | real*8 :: ivec(1:5) |
|---|
| 248 | |
|---|
| 249 | !----------- single-cell, lagrangian-type upwind rcon. ! |
|---|
| 250 | |
|---|
| 251 | do ipos = +2 , npos - 1 |
|---|
| 252 | |
|---|
| 253 | if (uvel(ipos) .gt. +0.d0) then |
|---|
| 254 | |
|---|
| 255 | !----------- integrate profile over upwind cell IPOS-1 ! |
|---|
| 256 | |
|---|
| 257 | if (size(SPAC).ne.+1) then |
|---|
| 258 | xhat = .5d0 * SPAC(ipos-1) |
|---|
| 259 | uCFL = uvel(ipos) & |
|---|
| 260 | & * tDEL / SPAC(ipos-1) |
|---|
| 261 | else |
|---|
| 262 | xhat = .5d0 * SPAC( +1) |
|---|
| 263 | uCFL = uvel(ipos) & |
|---|
| 264 | & * tDEL / SPAC( +1) |
|---|
| 265 | end if |
|---|
| 266 | |
|---|
| 267 | ss11 = +1.d0 - 2.d0 * uCFL |
|---|
| 268 | ss22 = +1.d0 |
|---|
| 269 | |
|---|
| 270 | call bfun1d(-1,mdof,ss11,vv11) |
|---|
| 271 | call bfun1d(-1,mdof,ss22,vv22) |
|---|
| 272 | |
|---|
| 273 | ivec = vv22 - vv11 |
|---|
| 274 | |
|---|
| 275 | do ivar = +1, nvar |
|---|
| 276 | |
|---|
| 277 | flux = dot_product ( & |
|---|
| 278 | & ivec(1:mdof), & |
|---|
| 279 | & QHAT(1:mdof,ivar,ipos-1) ) |
|---|
| 280 | |
|---|
| 281 | flux = flux * xhat |
|---|
| 282 | |
|---|
| 283 | qedg(ivar,ipos) = flux |
|---|
| 284 | |
|---|
| 285 | end do |
|---|
| 286 | |
|---|
| 287 | else & |
|---|
| 288 | & if (uvel(ipos) .lt. -0.d0) then |
|---|
| 289 | |
|---|
| 290 | !----------- integrate profile over upwind cell IPOS+0 ! |
|---|
| 291 | |
|---|
| 292 | if (size(SPAC).ne.+1) then |
|---|
| 293 | xhat = .5d0 * SPAC(ipos-0) |
|---|
| 294 | uCFL = uvel(ipos) & |
|---|
| 295 | & * tDEL / SPAC(ipos-0) |
|---|
| 296 | else |
|---|
| 297 | xhat = .5d0 * SPAC( +1) |
|---|
| 298 | uCFL = uvel(ipos) & |
|---|
| 299 | & * tDEL / SPAC( +1) |
|---|
| 300 | end if |
|---|
| 301 | |
|---|
| 302 | ss11 = -1.d0 - 2.d0 * uCFL |
|---|
| 303 | ss22 = -1.d0 |
|---|
| 304 | |
|---|
| 305 | call bfun1d(-1,mdof,ss11,vv11) |
|---|
| 306 | call bfun1d(-1,mdof,ss22,vv22) |
|---|
| 307 | |
|---|
| 308 | ivec = vv22 - vv11 |
|---|
| 309 | |
|---|
| 310 | do ivar = +1, nvar |
|---|
| 311 | |
|---|
| 312 | flux = dot_product ( & |
|---|
| 313 | & ivec(1:mdof), & |
|---|
| 314 | & QHAT(1:mdof,ivar,ipos-0) ) |
|---|
| 315 | |
|---|
| 316 | flux = flux * xhat |
|---|
| 317 | |
|---|
| 318 | qedg(ivar,ipos) = flux |
|---|
| 319 | |
|---|
| 320 | end do |
|---|
| 321 | |
|---|
| 322 | end if |
|---|
| 323 | |
|---|
| 324 | end do |
|---|
| 325 | |
|---|
| 326 | return |
|---|
| 327 | |
|---|
| 328 | end subroutine |
|---|
| 329 | |
|---|
| 330 | |
|---|
| 331 | |
|---|
