From d2ca2551dda911547deee09ca55e841e9daeb0b8 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 13 Dec 2015 08:48:20 -0200 Subject: [PATCH 1/8] INTERPOLATIONS: Make the extrema points monotonic in TVD. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 8d5eba4..b8d7821 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -484,12 +484,12 @@ module interpolations ! ! calculate the left- and right-side interface interpolations ! - do i = 1, n + do i = 2, n - 1 ! calculate left and right indices ! - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 ! calculate left and right side derivatives ! @@ -505,11 +505,14 @@ module interpolations fl(i ) = f(i) + df fr(im1) = f(i) - df - end do ! i = 1, n + end do ! i = 2, n - 1 ! update the interpolation of the first and last points ! - fl(1) = f(1) + i = n - 1 + fl(1) = 0.5d+00 * (f(1) + f(2)) + fr(i) = 0.5d+00 * (f(i) + f(n)) + fl(n) = f(n) fr(n) = f(n) !------------------------------------------------------------------------------- From a876581745a9aa7b42a2cf7b7d02624a769c2980 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 13 Dec 2015 08:53:15 -0200 Subject: [PATCH 2/8] INTERPOLATIONS: Make the extrema points monotonic in WENO3. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index b8d7821..b0530b5 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -566,12 +566,12 @@ module interpolations ! iterate along the vector ! - do i = 1, n + do i = 2, n - 1 ! prepare neighbour indices ! - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 ! calculate the left and right derivatives ! @@ -628,12 +628,15 @@ module interpolations ! fr(im1) = (wp * fp + wm * fm) / ww - end do ! i = 1, n + end do ! i = 2, n - 1 ! update the interpolation of the first and last points ! - fl(1) = f (1) - fr(n) = fl(n) + i = n - 1 + fl(1) = 0.5d+00 * (f(1) + f(2)) + fr(i) = 0.5d+00 * (f(i) + f(n)) + fl(n) = f(n) + fr(n) = f(n) !------------------------------------------------------------------------------- ! From dbe5fe71a0bf941a3e7eb8cfeb34fd94907e45a5 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 13 Dec 2015 08:59:21 -0200 Subject: [PATCH 3/8] INTERPOLATIONS: Make the extrema points monotonic in LIMO3. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index b0530b5..45bc64d 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -692,12 +692,12 @@ module interpolations ! iterate over positions and interpolate states ! - do i = 1, n + do i = 2, n - 1 ! prepare neighbour indices ! - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 ! prepare left and right differences ! @@ -773,12 +773,15 @@ module interpolations end if - end do ! i = 1, n + end do ! i = 2, n - 1 ! update the interpolation of the first and last points ! - fl(1) = f (1) - fr(n) = fl(n) + i = n - 1 + fl(1) = 0.5d+00 * (f(1) + f(2)) + fr(i) = 0.5d+00 * (f(i) + f(n)) + fl(n) = f(n) + fr(n) = f(n) !------------------------------------------------------------------------------- ! From f96c0d0297db4ad99f33d292d113778e275a80ff Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 13 Dec 2015 09:14:34 -0200 Subject: [PATCH 4/8] INTERPOLATIONS: Make the extrema points monotonic in WENO5Z. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 45bc64d..02d87d6 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -872,13 +872,13 @@ module interpolations ! iterate along the vector ! - do i = 1, n + do i = 2, n - 1 ! prepare neighbour indices ! - im1 = max(1, i - 1) im2 = max(1, i - 2) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 ip2 = min(n, i + 2) ! calculate βₖ (eqs. 9-11 in [1]) @@ -937,12 +937,15 @@ module interpolations ! fr(im1) = (wl * ql + wr * qr) + wc * qc - end do ! i = 1, n + end do ! i = 2, n - 1 ! update the interpolation of the first and last points ! - fl(1) = fr(1) - fr(n) = fl(n) + i = n - 1 + fl(1) = 0.5d+00 * (f(1) + f(2)) + fr(i) = 0.5d+00 * (f(i) + f(n)) + fl(n) = f(n) + fr(n) = f(n) !------------------------------------------------------------------------------- ! From 5b5efea431f6e4d4ec961f46170256dbc9fbd966 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 13 Dec 2015 10:56:28 -0200 Subject: [PATCH 5/8] INTERPOLATIONS: Make LIMO3 reconstruction symmetric. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 02d87d6..da8f3d9 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -715,11 +715,7 @@ module interpolations ! calculate values at i + ½ ! - if (dfr == 0.0d+00) then - - fl(i) = f(i) - - else + if (abs(dfr) > eps) then ! calculate the slope ratio (eq. 2.8 in [1]) ! @@ -741,15 +737,15 @@ module interpolations ! fl(i) = f(i) + dfr * (xl * f1 + xi * f2) + else + + fl(i) = f(i) + end if ! calculate values at i - ½ ! - if (dfl == 0.0d+00) then - - fr(im1) = f(i) - - else + if (abs(dfl) > eps) then ! calculate the slope ratio (eq. 2.8 in [1]) ! @@ -771,6 +767,10 @@ module interpolations ! fr(im1) = f(i) - dfl * (xl * f1 + xi * f2) + else + + fr(im1) = f(i) + end if end do ! i = 2, n - 1 From 28324447c4c8914553cf64c59ed490cc4fe8f083 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 13 Dec 2015 10:59:31 -0200 Subject: [PATCH 6/8] INTERPOLATIONS: Parameter eps should be at least 1.0d-12 for LIMO3. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index da8f3d9..78cc91a 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -184,6 +184,7 @@ module interpolations if (verbose .and. ng < 2) & call print_warning("interpolations:initialize_interpolation" & , "Increase the number of ghost cells (at least 2).") + eps = max(1.0d-12, eps) case ("weno5z", "weno5-z", "WENO5Z", "WENO5-Z") name_rec = "5th order WENO-Z (Borges et al. 2008)" reconstruct_states => reconstruct_weno5z From c22a5437b8e2c94b911dce6c9bfb537957635987 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 13 Dec 2015 11:10:31 -0200 Subject: [PATCH 7/8] INTERPOLATIONS: Make the extrema points monotonic in MP5. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 78cc91a..ed01497 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -2502,11 +2502,11 @@ module interpolations ! obtain the face values using high order interpolation ! - do i = 1, n + do i = 2, n - 1 im2 = max(1, i - 2) - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 ip2 = min(n, i + 2) fl(i) = (4.7d+01 * f(i ) + (2.7d+01 * f(ip1) - 1.3d+01 * f(im1)) & @@ -2516,14 +2516,14 @@ module interpolations - (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) & / 6.0d+01 - end do ! i = 1, n + end do ! i = 2, n - 1 ! apply monotonicity preserving limiting ! - do i = 1, n + do i = 2, n - 1 - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 if (dfm(i) * dfp(i) >= 0.0d+00) then sigma = kappa @@ -2593,12 +2593,15 @@ module interpolations ! fr(im1) = fr(i) - end do + end do ! n = 2, n - 1 ! update the interpolation of the first and last points ! - fl(1) = fr(1) - fr(n) = fl(n) + i = n - 1 + fl(1) = 0.5d+00 * (f(1) + f(2)) + fr(i) = 0.5d+00 * (f(i) + f(n)) + fl(n) = f(n) + fr(n) = f(n) !------------------------------------------------------------------------------- ! From 2f26b79e7977070d8353b321d8f926e9990fcd5e Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 13 Dec 2015 11:34:03 -0200 Subject: [PATCH 8/8] INTERPOLATIONS: Make the extrema points monotonic in CRMP5. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 131 +++++++++++++++++++++++------------------ 1 file changed, 75 insertions(+), 56 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index ed01497..f75fd03 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -2678,10 +2678,10 @@ module interpolations ! prepare the tridiagonal system coefficients for the interior ! - do i = 1, n + do i = ng, n - ng + 1 - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 a(i,1) = 3.0d-01 b(i,1) = 6.0d-01 @@ -2694,15 +2694,36 @@ module interpolations r(i,1) = (f(im1) + 1.9d+01 * f(i ) + 1.0d+01 * f(ip1)) / 3.0d+01 r(i,2) = (f(ip1) + 1.9d+01 * f(i ) + 1.0d+01 * f(im1)) / 3.0d+01 - end do ! i = 1, n + end do ! i = ng, n - ng + 1 ! interpolate ghost zones using explicit method (left-side reconstruction) ! - do i = 1, ng + do i = 2, ng im2 = max(1, i - 2) - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 + ip2 = i + 2 + + a(i,1) = 0.0d+00 + b(i,1) = 1.0d+00 + c(i,1) = 0.0d+00 + + r(i,1) = (4.7d+01 * f(i ) + (2.7d+01 * f(ip1) - 1.3d+01 * f(im1)) & + - (3.0d+00 * f(ip2) - 2.0d+00 * f(im2))) & + / 6.0d+01 + + end do ! i = 2, ng + a(1,1) = 0.0d+00 + b(1,1) = 1.0d+00 + c(1,1) = 0.0d+00 + r(1,1) = 0.5d+00 * (f(1) + f(2)) + + do i = n - ng, n - 1 + + im2 = i - 2 + im1 = i - 1 + ip1 = i + 1 ip2 = min(n, i + 2) a(i,1) = 0.0d+00 @@ -2713,32 +2734,40 @@ module interpolations - (3.0d+00 * f(ip2) - 2.0d+00 * f(im2))) & / 6.0d+01 - end do ! i = 1, ng - - do i = n - ng, n - - im2 = max(1, i - 2) - im1 = max(1, i - 1) - ip1 = min(n, i + 1) - ip2 = min(n, i + 2) - - a(i,1) = 0.0d+00 - b(i,1) = 1.0d+00 - c(i,1) = 0.0d+00 - - r(i,1) = (4.7d+01 * f(i ) + (2.7d+01 * f(ip1) - 1.3d+01 * f(im1)) & - - (3.0d+00 * f(ip2) - 2.0d+00 * f(im2))) & - / 6.0d+01 - - end do ! i = n - ng, n + end do ! i = n - ng, n - 1 + a(n,1) = 0.0d+00 + b(n,1) = 1.0d+00 + c(n,1) = 0.0d+00 + r(n,1) = f(n) ! interpolate ghost zones using explicit method (right-side reconstruction) ! - do i = 1, ng + 1 + do i = 2, ng + 1 im2 = max(1, i - 2) - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 + ip2 = i + 2 + + a(i,2) = 0.0d+00 + b(i,2) = 1.0d+00 + c(i,2) = 0.0d+00 + + r(i,2) = (4.7d+01 * f(i ) + (2.7d+01 * f(im1) - 1.3d+01 * f(ip1)) & + - (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) & + / 6.0d+01 + + end do ! i = 2, ng + 1 + a(1,2) = 0.0d+00 + b(1,2) = 1.0d+00 + c(1,2) = 0.0d+00 + r(1,2) = f(1) + + do i = n - ng + 1, n - 1 + + im2 = i - 2 + im1 = i - 1 + ip1 = i + 1 ip2 = min(n, i + 2) a(i,2) = 0.0d+00 @@ -2749,24 +2778,11 @@ module interpolations - (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) & / 6.0d+01 - end do ! i = 1, ng + 1 - - do i = n - ng + 1, n - - im2 = max(1, i - 2) - im1 = max(1, i - 1) - ip1 = min(n, i + 1) - ip2 = min(n, i + 2) - - a(i,2) = 0.0d+00 - b(i,2) = 1.0d+00 - c(i,2) = 0.0d+00 - - r(i,2) = (4.7d+01 * f(i ) + (2.7d+01 * f(im1) - 1.3d+01 * f(ip1)) & - - (3.0d+00 * f(im2) - 2.0d+00 * f(ip2))) & - / 6.0d+01 - - end do ! i = n - ng + 1, n + end do ! i = n - ng + 1, n - 1 + a(n,2) = 0.0d+00 + b(n,2) = 1.0d+00 + c(n,2) = 0.0d+00 + r(n,2) = 0.5d+00 * (f(n-1) + f(n)) ! solve the tridiagonal system of equations for the left-side interpolation ! @@ -2774,10 +2790,10 @@ module interpolations ! apply the monotonicity preserving limiting ! - do i = 1, n + do i = 2, n - 1 - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 if (dfm(i) * dfp(i) >= 0.0d+00) then sigma = kappa @@ -2814,7 +2830,7 @@ module interpolations end if - end do + end do ! i = 2, n - 1 ! solve the tridiagonal system of equations for the right-side interpolation ! @@ -2822,10 +2838,10 @@ module interpolations ! apply the monotonicity preserving limiting ! - do i = 1, n + do i = 2, n - 1 - im1 = max(1, i - 1) - ip1 = min(n, i + 1) + im1 = i - 1 + ip1 = i + 1 if (dfm(i) * dfp(i) >= 0.0d+00) then sigma = kappa @@ -2867,12 +2883,15 @@ module interpolations ! fr(im1) = fr(i) - end do + end do ! i = 2, n - 1 ! update the interpolation of the first and last points ! - fl(1) = fr(1) - fr(n) = fl(n) + i = n - 1 + fl(1) = 0.5d+00 * (f(1) + f(2)) + fr(i) = 0.5d+00 * (f(i) + f(n)) + fl(n) = f(n) + fr(n) = f(n) !------------------------------------------------------------------------------- !