From 9705ebf530e6ff740d86bb7c009dad2df8df432c Mon Sep 17 00:00:00 2001
From: Gordey Goyman <goyman@localhost.localdomain>
Date: Tue, 17 Sep 2024 13:20:37 +0300
Subject: [PATCH] cartesian metric on cube to test sat hordiff
---
.../sbp_diff2_variable_21_mod.f90 | 143 ++-----------
src/metric/metric_2d_cartesian_mod.f90 | 193 ++++++++++++++++++
src/metric/metric_factory_mod.f90 | 63 +++++-
.../test_hordiff/test_hordiff_scalar_mod.f90 | 10 +-
4 files changed, 274 insertions(+), 135 deletions(-)
create mode 100644 src/metric/metric_2d_cartesian_mod.f90
diff --git a/src/differential_operators/sbp_operators/sbp_diff2_variable_21_mod.f90 b/src/differential_operators/sbp_operators/sbp_diff2_variable_21_mod.f90
index 13c24361..d539f7f1 100644
--- a/src/differential_operators/sbp_operators/sbp_diff2_variable_21_mod.f90
+++ b/src/differential_operators/sbp_operators/sbp_diff2_variable_21_mod.f90
@@ -13,7 +13,7 @@ real(wp), parameter :: h11 = 1.0_wp / 2.0_wp, h22 = 1.0_wp, h33 = 1.0_wp
real(wp), parameter :: s11 = 3.0_wp / 2.0_wp, s12 = -2.0_wp, s13 = 1.0_wp / 2.0_wp !boundary operator + BS, B = diag(-1, 0, ..., 0, 1)
-real(wp), parameter :: q = 5.5_wp
+real(wp), parameter :: q = 0*5.5_wp
type, public, extends(sbp_diff2_variable_t) :: sbp_diff2_variable_21_t
contains
@@ -146,7 +146,6 @@ subroutine add_SAT_correction(this, f_out, f, df_norm, gI, b, mesh, scale_coefs,
integer(kind=4) :: k, i, j, is, ie, js, je
real(wp) :: l_flux, r_flux, df, dfg, df1, df2
real(wp) :: scale_coef_n, scale_coef_s, scale_coef_w, scale_coef_e, scale_coef
- real(wp) :: alpha_12, alpha_11, alpha_22, alpha1, alpha2
is = 1
ie = mesh%nx
@@ -165,55 +164,27 @@ subroutine add_SAT_correction(this, f_out, f, df_norm, gI, b, mesh, scale_coefs,
if (mesh%is<=is .and. mesh%ie>=is) then
do k = mesh%ks, mesh%ke
do j = mesh%js, mesh%je
- alpha_12 = mesh%J(is,j,k)*mesh%Qi(2,is,j,k)
- alpha1 = mesh%J(is,j,k)*mesh%Qi(1,is,j,k)
- alpha2 = mesh%J(is-1,j,k)*mesh%Qi(1,is-1,j,k)
l_flux = df_norm%p(is-1,j,k) !filling in halo zone
r_flux = df_norm%p(is,j,k)
df = (f%p(is,j,k) - f%p(is-1,j,k))
dfg = (f%p(is,j,k) - gI%p(is-1,j,k))
- f_out%p(is,j,k) = f_out%p(is,j,k) - b%p(is,j,k)*(alpha1 + alpha2)*0.25_wp*q/h11*df &
- + 0.5_wp*b%p(is,j,k)*(alpha1*s11)/h11*dfg &
- - 0.5_wp*(l_flux + r_flux)/h11
- if (j == js) then
- df1 = dfg
- df2 = (f%p(is,2,k) - gI%p(is-1,2,k))
- f_out%p(is,j,k) = f_out%p(is,j,k) - 0.5_wp*b%p(is,j,k)*alpha_12*(-df1 - 0.5_wp*df2)/h11
- else if (j == js+1) then
- df1 = (f%p(is,1,k) - gI%p(is-1,1,k))
- df2 = (f%p(is,3,k) - gI%p(is-1,3,k))
- f_out%p(is,j,k) = f_out%p(is,j,k) - 0.5_wp*b%p(is,j,k)*alpha_12*(df1 - 0.5_wp*df2)/h11
- else if (j == je) then
- df1 = dfg
- df2 = (f%p(is,je-1,k) - gI%p(is-1,je-1,k))
- f_out%p(is,j,k) = f_out%p(is,j,k) - 0.5_wp*b%p(is,j,k)*alpha_12*(df1 + 0.5_wp*df2)/h11
- else if (j == je-1) then
- df1 = (f%p(is,je,k) - gI%p(is-1,je,k))
- df2 = (f%p(is,je-2,k) - gI%p(is-1,je-2,k))
- f_out%p(is,j,k) = f_out%p(is,j,k) - 0.5_wp*b%p(is,j,k)*alpha_12*(-df1 + 0.5_wp*df2)/h11
- else
- df1 = (f%p(is,j+1,k) - gI%p(is-1,j+1,k))
- df2 = (f%p(is,j-1,k) - gI%p(is-1,j-1,k))
- f_out%p(is,j,k) = f_out%p(is,j,k) - 0.5_wp*b%p(is,j,k)*alpha_12*(0.5_wp*(df2 - df1))/h11
- end if
+ f_out%p(is,j,k) = f_out%p(is,j,k) + 0.5_wp*b%p(is,j,k)*(s11)/h11*dfg - 0.5_wp*(l_flux + r_flux)/h11
end do
end do
end if
if (mesh%is<=is+1 .and. mesh%ie>=is+1) then
do k = mesh%ks, mesh%ke
do j = mesh%js, mesh%je
- alpha_11 = mesh%J(is,j,k)*mesh%Qi(1,is,j,k)
dfg = (f%p(is,j,k) - gI%p(is-1,j,k))
- f_out%p(is+1,j,k) = f_out%p(is+1,j,k) + b%p(is,j,k)*0.5_wp*(alpha_11*s12)/h22*dfg
+ f_out%p(is+1,j,k) = f_out%p(is+1,j,k) + b%p(is,j,k)*0.5_wp*(s12)/h22*dfg
end do
end do
end if
if (mesh%is<=is+2 .and. mesh%ie>=is+2) then
do k = mesh%ks, mesh%ke
do j = mesh%js, mesh%je
- alpha_11 = mesh%J(is,j,k)*mesh%Qi(1,is,j,k)
dfg = (f%p(is,j,k) - gI%p(is-1,j,k))
- f_out%p(is+2,j,k) = f_out%p(is+2,j,k) + b%p(is,j,k)*alpha_11*0.5_wp*s13/h33*dfg
+ f_out%p(is+2,j,k) = f_out%p(is+2,j,k) + b%p(is,j,k)*0.5_wp*s13/h33*dfg
end do
end do
end if
@@ -221,56 +192,27 @@ subroutine add_SAT_correction(this, f_out, f, df_norm, gI, b, mesh, scale_coefs,
if (mesh%is<=ie .and. mesh%ie>=ie) then
do k = mesh%ks, mesh%ke
do j = mesh%js, mesh%je
- alpha_12 = mesh%J(ie,j,k)*mesh%Qi(2,ie,j,k)
- alpha1 = mesh%J(ie,j,k)*mesh%Qi(1,ie,j,k)
- alpha2 = mesh%J(ie+1,j,k)*mesh%Qi(1,ie+1,j,k)
l_flux = df_norm%p(ie,j,k)
r_flux = df_norm%p(ie+1,j,k)
df = (f%p(ie,j,k) - f%p(ie+1,j,k))
dfg = (f%p(ie,j,k) - gI%p(ie+1,j,k))
- f_out%p(ie,j,k) = f_out%p(ie,j,k) - b%p(ie,j,k)*(alpha1 + alpha2)*0.25_wp*q/h11*df &
- + b%p(ie,j,k)*0.5_wp*(alpha1*s11)/h11*dfg &
- - 0.5_wp*(l_flux + r_flux)/h11
-
- if (j == js) then
- df1 = dfg
- df2 = (f%p(ie,2,k) - gI%p(ie+1,2,k))
- f_out%p(ie,j,k) = f_out%p(ie,j,k) + 0.5_wp*b%p(ie,j,k)*alpha_12*(-df1 - 0.5_wp*df2)/h11
- else if (j == js+1) then
- df1 = (f%p(ie,1,k) - gI%p(ie+1,1,k))
- df2 = (f%p(ie,3,k) - gI%p(ie+1,3,k))
- f_out%p(ie,j,k) = f_out%p(ie,j,k) + 0.5_wp*b%p(ie,j,k)*alpha_12*(df1 - 0.5_wp*df2)/h11
- else if (j == je) then
- df1 = dfg
- df2 = (f%p(ie,je-1,k) - gI%p(ie+1,je-1,k))
- f_out%p(ie,j,k) = f_out%p(ie,j,k) + 0.5_wp*b%p(ie,j,k)*alpha_12*(df1 + 0.5_wp*df2)/h11
- else if (j == je-1) then
- df1 = (f%p(ie,je,k) - gI%p(ie+1,je,k))
- df2 = (f%p(ie,je-2,k) - gI%p(ie+1,je-2,k))
- f_out%p(ie,j,k) = f_out%p(ie,j,k) + 0.5_wp*b%p(ie,j,k)*alpha_12*(-df1 + 0.5_wp*df2)/h11
- else
- df1 = (f%p(ie,j+1,k) - gI%p(ie+1,j+1,k))
- df2 = (f%p(ie,j-1,k) - gI%p(ie+1,j-1,k))
- f_out%p(ie,j,k) = f_out%p(ie,j,k) + 0.5_wp*b%p(ie,j,k)*alpha_12*(0.5_wp*(df2 - df1))/h11
- end if
+ f_out%p(ie,j,k) = f_out%p(ie,j,k) + b%p(ie,j,k)*0.5_wp*(s11)/h11*dfg - 0.5_wp*(l_flux + r_flux)/h11
end do
end do
end if
if (mesh%is<=ie-1 .and. mesh%ie>=ie-1) then
do k = mesh%ks, mesh%ke
do j = mesh%js, mesh%je
- alpha_11 = mesh%J(ie,j,k)*mesh%Qi(1,ie,j,k)
dfg = (f%p(ie,j,k) - gI%p(ie+1,j,k))
- f_out%p(ie-1,j,k) = f_out%p(ie-1,j,k) + b%p(ie,j,k)*0.5_wp*(alpha_11*s12)/h22*dfg
+ f_out%p(ie-1,j,k) = f_out%p(ie-1,j,k) + b%p(ie,j,k)*0.5_wp*(s12)/h22*dfg
end do
end do
end if
if (mesh%is<=ie-2 .and. mesh%ie>=ie-2) then
do k = mesh%ks, mesh%ke
do j = mesh%js, mesh%je
- alpha_11 = mesh%J(ie,j,k)*mesh%Qi(1,ie,j,k)
dfg = (f%p(ie,j,k) - gI%p(ie+1,j,k))
- f_out%p(ie-2,j,k) = f_out%p(ie-2,j,k) + b%p(ie,j,k)*alpha_11*0.5_wp*s13/h33*dfg
+ f_out%p(ie-2,j,k) = f_out%p(ie-2,j,k) + b%p(ie,j,k)*0.5_wp*s13/h33*dfg
end do
end do
end if
@@ -279,57 +221,27 @@ subroutine add_SAT_correction(this, f_out, f, df_norm, gI, b, mesh, scale_coefs,
if (mesh%js<=js .and. mesh%je>=js) then
do k = mesh%ks, mesh%ke
do i = mesh%is, mesh%ie
- alpha_12 = mesh%J(i,js,k)*mesh%Qi(2,i,js,k)
- alpha1 = mesh%J(i,js,k)*mesh%Qi(3,i,js,k)
- alpha2 = mesh%J(i,js-1,k)*mesh%Qi(3,i,js-1,k)
-
l_flux = df_norm%p(i,js-1,k)!filling in halo zone
r_flux = df_norm%p(i,js,k)
df = (f%p(i,js,k) - f%p(i,js-1,k))
dfg = (f%p(i,js,k) - gI%p(i,js-1,k))
- f_out%p(i,js,k) = f_out%p(i,js,k) - b%p(i,js,k)*(alpha1 + alpha2)*0.25_wp*q/h11*df &
- + b%p(i,js,k)*0.5_wp*(alpha1*s11)/h11*dfg &
- - 0.5_wp*(l_flux + r_flux)/h11
-
- if (i == is) then
- df1 = dfg
- df2 = (f%p(2,js,k) - gI%p(2,js-1,k))
- f_out%p(i,js,k) = f_out%p(i,js,k) - 0.5_wp*b%p(i,js,k)*alpha_12*(-df1 - 0.5_wp*df2)/h11
- else if (i == is+1) then
- df1 = (f%p(1,js,k) - gI%p(1,js-1,k))
- df2 = (f%p(3,js,k) - gI%p(3,js-1,k))
- f_out%p(i,js,k) = f_out%p(i,js,k) - 0.5_wp*b%p(i,js,k)*alpha_12*(df1 - 0.5_wp*df2)/h11
- else if (i == ie) then
- df1 = dfg
- df2 = (f%p(ie-1,js,k) - gI%p(ie-1,js-1,k))
- f_out%p(i,js,k) = f_out%p(i,js,k) - 0.5_wp*b%p(i,js,k)*alpha_12*(df1 + 0.5_wp*df2)/h11
- else if (i == ie-1) then
- df1 = (f%p(ie,js,k) - gI%p(ie,js-1,k))
- df2 = (f%p(ie-2,js,k) - gI%p(ie-2,js-1,k))
- f_out%p(i,js,k) = f_out%p(i,js,k) - 0.5_wp*b%p(i,js,k)*alpha_12*(-df1 + 0.5_wp*df2)/h11
- else
- df1 = (f%p(i+1,js,k) - gI%p(i+1,js-1,k))
- df2 = (f%p(i-1,js,k) - gI%p(i-1,js-1,k))
- f_out%p(i,js,k) = f_out%p(i,js,k) - 0.5_wp*b%p(i,js,k)*alpha_12*(0.5_wp*(df2 - df1))/h11
- end if
+ f_out%p(i,js,k) = f_out%p(i,js,k) + b%p(i,js,k)*0.5_wp*(s11)/h11*dfg - 0.5_wp*(l_flux + r_flux)/h11
end do
end do
end if
if (mesh%js<=js+1 .and. mesh%je>=js+1) then
do k = mesh%ks, mesh%ke
do i = mesh%is, mesh%ie
- alpha_22 = mesh%J(i,js,k)*mesh%Qi(3,i,js,k)
dfg = (f%p(i,js,k) - gI%p(i,js-1,k))
- f_out%p(i,js+1,k) = f_out%p(i,js+1,k) + b%p(i,js,k)*0.5_wp*(alpha_22*s12)/h22*dfg
+ f_out%p(i,js+1,k) = f_out%p(i,js+1,k) + b%p(i,js,k)*0.5_wp*(s12)/h22*dfg
end do
end do
end if
if (mesh%js<=js+2 .and. mesh%je>=js+2) then
do k = mesh%ks, mesh%ke
do i = mesh%is, mesh%ie
- alpha_22 = mesh%J(i,js,k)*mesh%Qi(3,i,js,k)
dfg = (f%p(i,js,k) - gI%p(i,js-1,k))
- f_out%p(i,js+2,k) = f_out%p(i,js+2,k) + b%p(i,js,k)*alpha_22*0.5_wp*s13/h33*dfg
+ f_out%p(i,js+2,k) = f_out%p(i,js+2,k) + b%p(i,js,k)*0.5_wp*s13/h33*dfg
end do
end do
end if
@@ -337,56 +249,27 @@ subroutine add_SAT_correction(this, f_out, f, df_norm, gI, b, mesh, scale_coefs,
if (mesh%js<=je .and. mesh%je>=je) then
do k = mesh%ks, mesh%ke
do i = mesh%is, mesh%ie
- alpha_12 = mesh%J(i,je,k)*mesh%Qi(2,i,je,k)
- alpha1 = mesh%J(i,je,k)*mesh%Qi(3,i,je,k)
- alpha2 = mesh%J(i,je+1,k)*mesh%Qi(3,i,je+1,k)
l_flux = df_norm%p(i,je,k)
r_flux = df_norm%p(i,je+1,k)
df = (f%p(i,je,k) - f%p(i,je+1,k))
dfg = (f%p(i,je,k) - gI%p(i,je+1,k))
- f_out%p(i,je,k) = f_out%p(i,je,k) - b%p(i,je,k)*(alpha1 + alpha2)*0.25_wp*q/h11*df &
- + b%p(i,je,k)*0.5_wp*(alpha1*s11)/h11*dfg &
- - 0.5_wp*(l_flux + r_flux)/h11
-
- if (i == is) then
- df1 = dfg
- df2 = (f%p(2,je,k) - gI%p(2,je+1,k))
- f_out%p(i,je,k) = f_out%p(i,je,k) + 0.5_wp*b%p(i,je,k)*alpha_12*(-df1 - 0.5_wp*df2)/h11
- else if (i == is+1) then
- df1 = (f%p(1,je,k) - gI%p(1,je+1,k))
- df2 = (f%p(3,je,k) - gI%p(3,je+1,k))
- f_out%p(i,je,k) = f_out%p(i,je,k) + 0.5_wp*b%p(i,je,k)*alpha_12*(df1 - 0.5_wp*df2)/h11
- else if (i == ie) then
- df1 = dfg
- df2 = (f%p(ie-1,je,k) - gI%p(ie-1,je+1,k))
- f_out%p(i,je,k) = f_out%p(i,je,k) + 0.5_wp*b%p(i,je,k)*alpha_12*(df1 + 0.5_wp*df2)/h11
- else if (i == ie-1) then
- df1 = (f%p(ie,je,k) - gI%p(ie,je+1,k))
- df2 = (f%p(ie-2,je,k) - gI%p(ie-2,je+1,k))
- f_out%p(i,je,k) = f_out%p(i,je,k) + 0.5_wp*b%p(i,je,k)*alpha_12*(-df1 + 0.5_wp*df2)/h11
- else
- df1 = (f%p(i+1,je,k) - gI%p(i+1,je+1,k))
- df2 = (f%p(i-1,je,k) - gI%p(i-1,je+1,k))
- f_out%p(i,je,k) = f_out%p(i,je,k) + 0.5_wp*b%p(i,je,k)*alpha_12*(0.5_wp*(df2 - df1))/h11
- end if
+ f_out%p(i,je,k) = f_out%p(i,je,k) + b%p(i,je,k)*0.5_wp*(s11)/h11*dfg - 0.5_wp*(l_flux + r_flux)/h11
end do
end do
end if
if (mesh%js<=je-1 .and. mesh%je>=je-1) then
do k = mesh%ks, mesh%ke
do i = mesh%is, mesh%ie
- alpha_22 = mesh%J(i,je,k)*mesh%Qi(3,i,je,k)
dfg = (f%p(i,je,k) - gI%p(i,je+1,k))
- f_out%p(i,je-1,k) = f_out%p(i,je-1,k) + b%p(i,je,k)*0.5_wp*(alpha_22*s12)/h22*dfg
+ f_out%p(i,je-1,k) = f_out%p(i,je-1,k) + b%p(i,je,k)*0.5_wp*(s12)/h22*dfg
end do
end do
end if
if (mesh%js<=je-2 .and. mesh%je>=je-2) then
do k = mesh%ks, mesh%ke
do i = mesh%is, mesh%ie
- alpha_22 = mesh%J(i,je,k)*mesh%Qi(3,i,je,k)
dfg = (f%p(i,je,k) - gI%p(i,je+1,k))
- f_out%p(i,je-2,k) = f_out%p(i,je-2,k) + b%p(i,je,k)*alpha_22*0.5_wp*s13/h33*dfg
+ f_out%p(i,je-2,k) = f_out%p(i,je-2,k) + b%p(i,je,k)*0.5_wp*s13/h33*dfg
end do
end do
end if
diff --git a/src/metric/metric_2d_cartesian_mod.f90 b/src/metric/metric_2d_cartesian_mod.f90
new file mode 100644
index 00000000..a52f3a17
--- /dev/null
+++ b/src/metric/metric_2d_cartesian_mod.f90
@@ -0,0 +1,193 @@
+module metric_2d_cartesian_mod
+
+ use metric_2d_mod, only : metric_2d_t
+ use topology_mod, only : topology_t
+ use parcomm_mod, only : parcomm_global
+ use mesh_mod, only : mesh_t
+
+ implicit none
+ ! planar cartesian metric.
+ ! alpha = x, beta = y, z = 0
+ type, extends(metric_2d_t) :: metric_2d_cartesian_t
+ class(topology_t), allocatable :: topology
+ contains
+ procedure :: calc_r => calculate_cartesian_r
+
+ procedure :: calc_a1 => calculate_cartesian_a1
+ procedure :: calc_a2 => calculate_cartesian_a2
+
+ procedure :: calc_b1 => calculate_cartesian_b1
+ procedure :: calc_b2 => calculate_cartesian_b2
+
+ procedure :: calc_Q => calculate_cartesian_Q
+ procedure :: calc_Qi => calculate_cartesian_Qi
+
+ procedure :: calc_J => calculate_cartesian_J
+
+ procedure :: calc_G => calculate_cartesian_G
+
+ procedure :: transform_cart2native => transform_cartesian_to_native
+
+ procedure :: set_curvilinear_mesh
+
+ end type metric_2d_cartesian_t
+
+ contains
+
+ pure function calculate_cartesian_r(this, panel_ind, alpha, beta) result(r)
+ class(metric_2d_cartesian_t), intent(in) :: this
+ integer(kind=4), intent(in) :: panel_ind
+ real(kind=8), intent(in) :: alpha, beta
+ real(kind=8) :: r(3)
+
+ r = [alpha, beta, 0.0_8]
+
+ end function calculate_cartesian_r
+
+ pure function calculate_cartesian_a1(this, panel_ind, alpha, beta) result(a1)
+ class(metric_2d_cartesian_t), intent(in) :: this
+ integer(kind=4), intent(in) :: panel_ind
+ real(kind=8), intent(in) :: alpha, beta
+ real(kind=8) :: a1(3)
+
+ a1 = [1.0_8, 0.0_8, 0.0_8]
+
+ end function calculate_cartesian_a1
+
+ pure function calculate_cartesian_a2(this, panel_ind, alpha, beta) result(a2)
+ class(metric_2d_cartesian_t), intent(in) :: this
+ integer(kind=4), intent(in) :: panel_ind
+ real(kind=8), intent(in) :: alpha, beta
+ real(kind=8) :: a2(3)
+
+ a2 = [0.0_8, 1.0_8, 0.0_8]
+
+ end function calculate_cartesian_a2
+
+ pure function calculate_cartesian_b1(this, panel_ind, alpha, beta) result(b1)
+ class(metric_2d_cartesian_t), intent(in) :: this
+ integer(kind=4), intent(in) :: panel_ind
+ real(kind=8), intent(in) :: alpha, beta
+ real(kind=8) :: b1(3)
+
+ b1 = [1.0_8, 0.0_8, 0.0_8]
+
+ end function calculate_cartesian_b1
+
+ pure function calculate_cartesian_b2(this, panel_ind, alpha, beta) result(b2)
+ class(metric_2d_cartesian_t), intent(in) :: this
+ integer(kind=4), intent(in) :: panel_ind
+ real(kind=8), intent(in) :: alpha, beta
+ real(kind=8) :: b2(3)
+
+ b2 = [0.0_8, 1.0_8, 0.0_8]
+
+ end function calculate_cartesian_b2
+
+ pure function calculate_cartesian_Q(this, panel_ind, alpha, beta) result(Q)
+ !Calculate metric tensor
+ !Q = |a1*a1 a1*a2| = |Q(1) Q(2)|
+ ! |a1*a2 a2*a2| |Q(2) Q(3)|
+ class(metric_2d_cartesian_t), intent(in) :: this
+ integer(kind=4), intent(in) :: panel_ind
+ real(kind=8), intent(in) :: alpha, beta
+ real(kind=8) :: Q(3)
+
+ Q = [1.0_8, 0.0_8, 1.0_8]
+
+ end function calculate_cartesian_Q
+
+ pure function calculate_cartesian_Qi(this,panel_ind, alpha, beta) result(Qi)
+ !Calculate inverse metric tensor
+ !QI = inv |a1*a1 a1*a2| = |QI(1) QI(2)|
+ ! |a1*a2 a2*a2| |QI(2) QI(3)|
+ class(metric_2d_cartesian_t), intent(in) :: this
+ integer(kind=4), intent(in) :: panel_ind
+ real(kind=8), intent(in) :: alpha,beta
+ real(kind=8) :: Qi(3)
+
+ Qi = [1.0_8, 0.0_8, 1.0_8]
+
+ end function calculate_cartesian_Qi
+
+ pure function calculate_cartesian_G(this, panel_ind, alpha, beta) result(G)
+ !Calculate christoffel symbols
+ class(metric_2d_cartesian_t), intent(in) :: this
+ integer(kind=4), intent(in) :: panel_ind
+ real(kind=8), intent(in) :: alpha,beta
+ real(kind=8) :: G(2, 2, 2)
+
+ G(1:2, 1:2, 1:2) = 0.0_8
+
+ end function calculate_cartesian_G
+
+ pure function calculate_cartesian_J(this, panel_ind, alpha, beta) result(J)
+ !Compute sqrt of metric tensor determinant
+ class(metric_2d_cartesian_t), intent(in) :: this
+ integer(kind=4), intent(in) :: panel_ind
+ real(kind=8), intent(in) :: alpha, beta
+ real(kind=8) :: J
+
+ J = 1.0_8
+
+ end function calculate_cartesian_J
+
+ subroutine transform_cartesian_to_native(this, panel_ind, alpha, beta, r)
+ class(metric_2d_cartesian_t), intent(in) :: this
+ integer(kind=4), intent(out) :: panel_ind
+ real(kind=8), intent(out) :: alpha, beta
+ real(kind=8), intent(in) :: r(3)
+
+ alpha = r(1)
+ beta = r(2)
+ panel_ind = 1
+
+ end subroutine transform_cartesian_to_native
+
+ subroutine set_curvilinear_mesh(this, mesh)
+
+ use const_mod, only : pi
+
+ class(metric_2d_cartesian_t), intent(in) :: this
+ type(mesh_t), intent(inout) :: mesh
+
+ integer(kind=4) :: t, i, j, k, ks, ke, panel_ind, halo_width
+ real(kind=8) :: alpha, beta, ta, tb, sigm, r(3), a1(3), a2(3), b1(3), b2(3)
+
+ do t = mesh%ts, mesh%te
+ halo_width = mesh%tile(t)%halo_width
+ ks = mesh%tile(t)%ks
+ ke = mesh%tile(t)%ke
+ panel_ind = mesh%tile(t)%panel_ind
+ do j = mesh%tile(t)%js-halo_width, mesh%tile(t)%je+halo_width
+ beta = mesh%tile(t)%get_beta(j)
+ do i = mesh%tile(t)%is-halo_width, mesh%tile(t)%ie+halo_width
+ alpha = mesh%tile(t)%get_alpha(i)
+
+ mesh%tile(t)%rx(i,j,ks) = alpha
+ mesh%tile(t)%ry(i,j,ks) = beta
+ mesh%tile(t)%rz(i,j,ks) = 0.0_8
+
+ mesh%tile(t)%h(i,j,ks) = 0.0_8
+
+ mesh%tile(t)%a1(1:4,i,j,ks) = [1.0_8, 0.0_8, 0.0_8, 0.0_8]
+ mesh%tile(t)%a2(1:4,i,j,ks) = [0.0_8, 1.0_8, 0.0_8, 0.0_8]
+ mesh%tile(t)%a3(1:4,i,j,ks) = [0.0_8, 0.0_8, 0.0_8, 1.0_8]
+
+ mesh%tile(t)%b1(1:3,i,j,ks) = [1.0_8, 0.0_8, 0.0_8]
+ mesh%tile(t)%b2(1:3,i,j,ks) = [0.0_8, 1.0_8, 0.0_8]
+
+ mesh%tile(t)%b3(1:4,i,j,ks) = [0.0_8, 0.0_8, 0.0_8, 1.0_8]
+
+ mesh%tile(t)%Q(1:3,i,j,ks) = [1.0_8, 0.0_8, 1.0_8]
+
+ mesh%tile(t)%Qi(1:3,i,j,ks) = [1.0_8, 0.0_8, 1.0_8]
+
+ mesh%tile(t)%J(i,j,ks) = 1.0_8
+ end do
+ end do
+ end do
+
+ end subroutine set_curvilinear_mesh
+
+ end module metric_2d_cartesian_mod
diff --git a/src/metric/metric_factory_mod.f90 b/src/metric/metric_factory_mod.f90
index a02424c5..c8556b8f 100644
--- a/src/metric/metric_factory_mod.f90
+++ b/src/metric/metric_factory_mod.f90
@@ -112,7 +112,7 @@ subroutine create_shallow_atmosphere_metric(metric, topology, config)
allocate(metric_shallow_atm)
- call config%get(metric_2d_type,"metric_2d_type",default=METRIC_2D_DEFAULT)
+ call config%get(metric_2d_type,"metric_2d_type",default="cartesian")
call create_metric_2d(metric_shallow_atm%metric_2d, topology, metric_2d_type, config)
metric_shallow_atm%scale = metric_shallow_atm%metric_2d%scale
@@ -240,6 +240,9 @@ subroutine create_metric_2d(metric, topology, metric_type, config)
case("ecs","equiangular_cubed_sphere")
call create_ecs_2d_metric(metric, topology, &
scale, omega, rotation_matrix, rotation_axis)
+ case("cartesian")
+ call create_cartesian_2d_metric(metric, topology, &
+ scale, omega, rotation_matrix, rotation_axis)
case default
call parcomm_global%abort("Unknown metric_type var " // metric_type // &
" in metric_factory_mod")
@@ -345,5 +348,63 @@ subroutine create_ecs_2d_metric(metric, topology, sphere_r, sphere_omega, rotati
end subroutine create_ecs_2d_metric
+subroutine create_cartesian_2d_metric(metric, topology, sphere_r, sphere_omega, rotation_matrix, rotation_axis)
+
+ use topology_mod, only : topology_t
+ use metric_2d_mod, only : metric_2d_t
+ use metric_2d_cartesian_mod, only : metric_2d_cartesian_t
+ use const_mod, only : pi
+
+ class(topology_t), intent(in) :: topology
+ class(metric_2d_t), allocatable, intent(out) :: metric
+ real(kind=8), optional, intent(in) :: sphere_r, sphere_omega
+ real(kind=8), optional, intent(in) :: rotation_matrix(3,3), rotation_axis(3)
+
+ type(metric_2d_cartesian_t), allocatable :: cartesian_metric
+ real(kind=8) :: local_rotation_matrix(3, 3), local_rotation_axis(3)
+ real(kind=8) :: local_sphere_r, local_sphere_omega
+
+ local_rotation_matrix = reshape([1,0,0, 0,1,0, 0,0,1], [3,3])
+ if (present(rotation_matrix)) local_rotation_matrix = rotation_matrix
+
+ local_rotation_axis = [0, 0, 1]
+ if (present(rotation_axis)) local_rotation_axis = rotation_axis
+
+ local_sphere_omega = 1.0_8
+ if (present(sphere_omega)) local_sphere_omega = sphere_omega
+
+ local_sphere_r = 1.0_8
+ if (present(sphere_r)) local_sphere_r = sphere_r
+
+ allocate(cartesian_metric)
+
+ cartesian_metric%scale = local_sphere_r
+ cartesian_metric%omega = local_sphere_omega
+
+ cartesian_metric%rotation_matrix = local_rotation_matrix
+ cartesian_metric%rotation_axis = local_rotation_axis
+
+ cartesian_metric%topology = topology
+ cartesian_metric%alpha0 = 0.0_8
+ cartesian_metric%beta0 = 0.0_8
+ cartesian_metric%alpha1 = 1.0_8
+ cartesian_metric%beta1 = 1.0_8
+
+ ! select type(topology)
+ ! class is (torus_topology_t)
+ ! cartesian_metric%topology = topology
+ ! cartesian_metric%alpha0 = 0.0_8
+ ! cartesian_metric%beta0 = 0.0_8
+ ! cartesian_metric%alpha1 = 1.0_8 / topology%panels_x_num
+ ! cartesian_metric%beta1 = 1.0_8 / topology%panels_x_num
+ ! class default
+ ! call parcomm_global%abort("Wrong topology class in create_cartesian_2d_metric!")
+ ! end select
+
+ call move_alloc(cartesian_metric, metric)
+
+end subroutine create_cartesian_2d_metric
+
+
end module metric_factory_mod
diff --git a/src/test/test_hordiff/test_hordiff_scalar_mod.f90 b/src/test/test_hordiff/test_hordiff_scalar_mod.f90
index eb56eaa4..c6049144 100644
--- a/src/test/test_hordiff/test_hordiff_scalar_mod.f90
+++ b/src/test/test_hordiff/test_hordiff_scalar_mod.f90
@@ -48,26 +48,28 @@ subroutine hordiff_scalar_test()
call create_grid_field(h, halo_width, 0, domain%mesh_p)
call create_grid_field(h_tend, halo_width, 0, domain%mesh_p)
- call create_hordiff_operator(diff_op, "hordiff_scalar_Ah_sbp_63_narrow", diff_coeff, domain)
+ call create_hordiff_operator(diff_op, "hordiff_scalar_Ah_sbp_sat_21_narrow", diff_coeff, domain)
call set_scalar_test_field(h, field, domain%mesh_p, 0)
+ ! call h%assign(1.0_8, domain%mesh_p)
+
if(trim(staggering)=="Ah") then
call create_halo_procedure(Ah_sync, domain, 1, "Ah_scalar_sync")
call Ah_sync%get_halo_scalar(h, domain, 1)
end if
- call create_latlon_outputer(outputer, 2*N+1, 4*N, "o", domain)
+ ! call create_latlon_outputer(outputer, 2*N+1, 4*N, "o", domain)
call create_master_paneled_outputer(outputer_mp, "p", domain)
- call create_quadrature(quadrature, "SBP_Ah63_quadrature", domain%mesh_p)
+ call create_quadrature(quadrature, "SBP_Ah21_quadrature", domain%mesh_p)
do it = 1, 600
if (domain%parcomm%myid==0) print*, it
- call outputer%write(h, domain, 'h_diff.dat', it)
+ ! call outputer%write(h, domain, 'h_diff.dat', it)
call outputer_mp%write(h, domain, 'h_diff_mp.dat', it)
--
GitLab