From abbfb5339858296d5897e4aaf744d222e9d3bc93 Mon Sep 17 00:00:00 2001 From: Knut Date: Tue, 5 Jul 2022 15:47:35 +0200 Subject: [PATCH 1/3] diffusion1: test_tridiagonal --- tests/CMakeLists.txt | 9 ++++++-- tests/test_tridiagonal.F90 | 45 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+), 2 deletions(-) create mode 100644 tests/test_tridiagonal.F90 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1084beac9..997cc1569 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -25,7 +25,12 @@ add_executable(test_time EXCLUDE_FROM_ALL ) target_link_libraries(test_time PRIVATE util) +add_executable(test_tridiagonal EXCLUDE_FROM_ALL + test_tridiagonal.F90 + ) +target_link_libraries(test_tridiagonal PRIVATE util) + add_custom_target(gotm_tests) -add_dependencies(gotm_tests test_airsea test_albedo test_eqstate test_time) +add_dependencies(gotm_tests test_airsea test_albedo test_eqstate test_time test_tridiagonal) -set_property(TARGET gotm_tests test_airsea test_albedo test_bulk test_eqstate test_time PROPERTY FOLDER tests) +set_property(TARGET gotm_tests test_airsea test_albedo test_bulk test_eqstate test_time test_tridiagonal PROPERTY FOLDER tests) diff --git a/tests/test_tridiagonal.F90 b/tests/test_tridiagonal.F90 new file mode 100644 index 000000000..46c0af607 --- /dev/null +++ b/tests/test_tridiagonal.F90 @@ -0,0 +1,45 @@ +#include "cppdefs.h" +! + program test_tridiagonal +! +! !DESCRIPTION: +! program to test the tridiagonal inter action routines +! +! To execute: +! make test_tridiagonal +! +! !USES: + use mtridiagonal, only: init_tridiagonal, tridiagonal, clean_tridiagonal + use mtridiagonal, only: au, bu, cu, du + implicit none +! +! !LOCAL VARIABLES + integer, parameter :: rk = kind(_ONE_) + integer, parameter :: Nmax = 4 + integer :: N + REALTYPE :: res(0:Nmax) + +!EOP +!----------------------------------------------------------------------- +!BOC + + LEVEL1 "tridiag( 1 , 2 , 1 ) X = 1" + do N = Nmax,0,-1 + LEVEL0 LINE + LEVEL1 "N = ", N + call init_tridiagonal(N) + au = 1._rk + bu = 2._rk + cu = 1._rk + du = 1._rk + call tridiagonal(N, 0, N, res(0:N)) + LEVEL1 "X = ", real(res(0:N)) + call clean_tridiagonal() + end do + + end program test_tridiagonal +!EOC + +!----------------------------------------------------------------------- +! Copyright by the GOTM-team under the GNU Public License - www.gnu.org +!----------------------------------------------------------------------- From 50a6784a73da6190f34ced975ecc978d2911149d Mon Sep 17 00:00:00 2001 From: Knut Date: Tue, 5 Jul 2022 16:13:09 +0200 Subject: [PATCH 2/3] diffusion1: bugfix - tridiagonal now also works for one element --- src/util/tridiagonal.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/util/tridiagonal.F90 b/src/util/tridiagonal.F90 index a8323db66..742095754 100644 --- a/src/util/tridiagonal.F90 +++ b/src/util/tridiagonal.F90 @@ -117,6 +117,7 @@ subroutine tridiagonal(N,fi,lt,value) !EOP ! ! !LOCAL VARIABLES: + REALTYPE :: tmp integer :: i ! !----------------------------------------------------------------------- @@ -124,13 +125,12 @@ subroutine tridiagonal(N,fi,lt,value) ru(lt)=au(lt)/bu(lt) qu(lt)=du(lt)/bu(lt) - do i=lt-1,fi+1,-1 - ru(i)=au(i)/(bu(i)-cu(i)*ru(i+1)) - qu(i)=(du(i)-cu(i)*qu(i+1))/(bu(i)-cu(i)*ru(i+1)) + do i=lt-1,fi,-1 + tmp = _ONE_ / (bu(i)-cu(i)*ru(i+1)) + ru(i)=au(i) * tmp + qu(i)=(du(i)-cu(i)*qu(i+1)) * tmp end do - qu(fi)=(du(fi)-cu(fi)*qu(fi+1))/(bu(fi)-cu(fi)*ru(fi+1)) - value(fi)=qu(fi) do i=fi+1,lt value(i)=qu(i)-ru(i)*value(i-1) From 7d49c3d822155430ff8af360e5a36c4f21d18c90 Mon Sep 17 00:00:00 2001 From: Knut Date: Thu, 7 Jul 2022 15:17:16 +0200 Subject: [PATCH 3/3] diffusion1: rewrite of diffusion now also works for one layer --- src/util/CMakeLists.txt | 1 + src/util/diff_center.F90 | 8 +- src/util/diff_face.F90 | 20 +++- src/util/diffusion.F90 | 212 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 239 insertions(+), 2 deletions(-) create mode 100644 src/util/diffusion.F90 diff --git a/src/util/CMakeLists.txt b/src/util/CMakeLists.txt index 7d81cb632..52aad51ce 100644 --- a/src/util/CMakeLists.txt +++ b/src/util/CMakeLists.txt @@ -9,6 +9,7 @@ configure_file(gotm_version.F90.in gotm_version.F90) add_library(util adv_center.F90 convert_fluxes.F90 + diffusion.F90 diff_center.F90 diff_face.F90 eqstate.F90 diff --git a/src/util/diff_center.F90 b/src/util/diff_center.F90 index 3d4059da1..6fe64014c 100644 --- a/src/util/diff_center.F90 +++ b/src/util/diff_center.F90 @@ -114,10 +114,15 @@ subroutine diff_center(N,dt,cnpar,posconc,h,Bcup,Bcdw, & ! !LOCAL VARIABLES: integer :: i REALTYPE :: a,c,l + REALTYPE :: Af(0:N) ! !----------------------------------------------------------------------- !BOC -! +#if 1 + Af = _ONE_ + call diffusion(N, dt, cnpar, posconc, h, h, Af, Bcup, Bcdw, & + Yup, Ydw, nuY, Lsour, Qsour, Taur, Yobs, Y) +#else ! set up matrix do i=2,N-1 c = 2.0d0*dt*nuY(i) /(h(i)+h(i+1))/h(i) @@ -191,6 +196,7 @@ subroutine diff_center(N,dt,cnpar,posconc,h,Bcup,Bcdw, & ! solve linear system call tridiagonal(N,1,N,Y) +#endif return end subroutine diff_center diff --git a/src/util/diff_face.F90 b/src/util/diff_face.F90 index 4e2fb3893..4ec5b7546 100644 --- a/src/util/diff_face.F90 +++ b/src/util/diff_face.F90 @@ -64,10 +64,27 @@ subroutine diff_face(N,dt,cnpar,h,Bcup,Bcdw,Yup,Ydw,nuY,Lsour,Qsour,Y) ! !LOCAL VARIABLES: integer :: i REALTYPE :: a,c,l +#ifdef _NEW_DIFF_FACE_ + REALTYPE :: hf(0:N-1) + REALTYPE :: nuYc(0:N-1) + REALTYPE :: Taur(0:N-1) + REALTYPE :: Yobs(0:N-1) + REALTYPE :: Ac(0:N-1) +#endif ! !----------------------------------------------------------------------- !BOC -! + + if (N .eq. 1) return + +#ifdef _NEW_DIFF_FACE_ + hf (1:N-1) = _HALF_ * ( h (1:N-1) + h (2:N ) ) + nuYc(1:N-2) = _HALF_ * ( nuY(1:N-2) + nuY(2:N-1) ) + Taur = 1.0d15 + Ac = _ONE_ + call diffusion(N-1, dt, cnpar, 0, hf, hf, Ac, Bcup, Bcdw, & + Yup, Ydw, nuYc, Lsour, Qsour, Taur, Yobs, Y(0:N-1)) +#else ! set up matrix do i=2,N-2 c = dt*( nuY(i+1) + nuY(i ) ) / ( h(i)+h(i+1) ) / h(i+1) @@ -124,6 +141,7 @@ subroutine diff_face(N,dt,cnpar,h,Bcup,Bcdw,Yup,Ydw,nuY,Lsour,Qsour,Y) ! solve linear system call tridiagonal(N,1,N-1,Y) +#endif return end subroutine diff_face diff --git a/src/util/diffusion.F90 b/src/util/diffusion.F90 new file mode 100644 index 000000000..a6d0a6869 --- /dev/null +++ b/src/util/diffusion.F90 @@ -0,0 +1,212 @@ +#include"cppdefs.h" +!----------------------------------------------------------------------- +!BOP +! +! !ROUTINE: Diffusion schemes --- grid centers\label{sec:diffusionMean} +! +! !INTERFACE: + subroutine diffusion(N, dt, cnpar, posconc, h, Vc, Af, Bcup, Bcdw, & + Yup, Ydw, nuY, Lsour, Qsour, Taur, Yobs, Y) +! +! !DESCRIPTION: +! This subroutine solves the one-dimensional diffusion equation +! including source terms, +! \begin{equation} +! \label{YdiffCenter} +! \partder{Y}{t} +! = \partder{}{z} \left( \nu_Y \partder{Y}{z} \right) +! - \frac{1}{\tau_R}(Y-Y_{obs}) +! + Y L_{\text{sour}} + Q_{\text{sour}} +! \comma +! \end{equation} +! for al variables defined at the centers of the grid cells, and +! a diffusion coefficient $\nu_Y$ defined at the faces. +! Relaxation with time scale $\tau_R$ towards observed values +! $Y_{\text{obs}}$ is possible. $L_{\text{sour}}$ specifies a +! linear source term, and $Q_{\text{sour}}$ a constant source term. +! Central differences are used to discretize the problem +! as discussed in \sect{SectionNumericsMean}. The diffusion term, +! the linear source term, and the linear part arising from the +! relaxation term are treated +! with an implicit method, whereas the constant source term is treated +! fully explicit. +! +! The input parameters {\tt Bcup} and {\tt Bcdw} specify the type +! of the upper and lower boundary conditions, which can be either +! Dirichlet or Neumann-type. {\tt Bcup} and {\tt Bcdw} must have integer +! values corresponding to the parameters {\tt Dirichlet} and {\tt Neumann} +! defined in the module {\tt util}, see \sect{sec:utils}. +! {\tt Yup} and {\tt Ydw} are the values of the boundary conditions at +! the surface and the bottom. Depending on the values of {\tt Bcup} and +! {\tt Bcdw}, they represent either fluxes or prescribed values. +! The integer {\tt posconc} indicates if a quantity is +! non-negative by definition ({\tt posconc}=1, such as for concentrations) +! or not ({\tt posconc}=0). For {\tt posconc}=1 and negative +! boundary fluxes, the source term linearisation according to +! \cite{Patankar80} is applied. +! +! Note that fluxes \emph{entering} a boundary cell are counted positive +! by convention. The lower and upper position for prescribing these fluxes +! are located at the lowest und uppermost grid faces with index "0" and +! index "N", respectively. If values are prescribed, they are located at +! the centers with index "1" and index "N", respectively. +! +! !USES: + use util, only : Dirichlet, Neumann + use mtridiagonal + + IMPLICIT NONE +! +! !INPUT PARAMETERS: + +! number of vertical layers + integer, intent(in) :: N + +! time step (s) + REALTYPE, intent(in) :: dt + +! "implicitness" parameter + REALTYPE, intent(in) :: cnpar + +! 1: non-negative concentration, 0: else + integer, intent(in) :: posconc + +! layer thickness (m) + REALTYPE, intent(in) :: h(0:N) + +! cell volume + REALTYPE, intent(in) :: Vc(0:N) + +! interface area + REALTYPE, intent(in) :: Af(0:N) + +! type of upper BC + integer, intent(in) :: Bcup + +! type of lower BC + integer, intent(in) :: Bcdw + +! value of upper BC + REALTYPE, intent(in) :: Yup + +! value of lower BC + REALTYPE, intent(in) :: Ydw + +! diffusivity of Y + REALTYPE, intent(in) :: nuY(0:N) + +! linear source term +! (treated implicitly) + REALTYPE, intent(in) :: Lsour(0:N) + +! constant source term +! (treated explicitly) + REALTYPE, intent(in) :: Qsour(0:N) + +! relaxation time (s) + REALTYPE, intent(in) :: Taur(0:N) + +! observed value of Y + REALTYPE, intent(in) :: Yobs(0:N) +! +! !INPUT/OUTPUT PARAMETERS: + REALTYPE, intent(inout) :: Y(0:N) +! +! !REVISION HISTORY: +! Original author(s): Lars Umlauf +! +!EOP +! +! !LOCAL VARIABLES: + integer :: i + REALTYPE :: dV(0:N), dY(0:N) + REALTYPE :: fac, wup, wdw +! +!----------------------------------------------------------------------- +!BOC +! +! set up matrix + dV(0) = _ZERO_ + dY(0) = _ZERO_ + dV(N) = _ZERO_ + dY(N) = _ZERO_ + do concurrent ( i = 1 : N-1 ) + dV(i) = dt * Af(i) * nuY(i) / ( _HALF_ * ( h(i) + h(i+1) ) ) + dY (i) = Y(i+1) - Y(i) + end do + + do concurrent ( i = 1 : N ) + cu(i) = - cnpar * dV(i ) + au(i) = - cnpar * dV(i-1) + bu(i) = Vc(i) * ( _ONE_ - dt*Lsour(i) ) & + + cnpar * ( dV(i-1) + dV(i) ) + du(i) = Vc(i) * ( Y(i) + dt*Qsour(i) ) & + + ( _ONE_ - cnpar ) * ( dV(i)*dY(i) - dV(i-1)*dY(i-1) ) + end do + +! set up upper boundary condition + wup = _ZERO_ + select case(Bcup) + case(Neumann) + if (posconc.eq.1 .and. Yup.lt._ZERO_) then ! Patankar (1980) trick + fac = -Y(N) / ( dt * Af(N) * Yup ) + if ( fac.gt.TINY(_ONE_) ) bu(N) = bu(N) + _ONE_/fac + else + du(N) = du(N) + dt * Af(N) * Yup + end if + case(Dirichlet) + au(N) = _ZERO_ + bu(N) = Vc(N) + du(N) = Vc(N) * Yup + wup = _ONE_ + case default + FATAL 'invalid boundary condition type for upper boundary' + stop 'diffusion.F90' + end select + +! set up lower boundary condition + wdw = _ZERO_ + select case(Bcdw) + case(Neumann) + if (posconc.eq.1 .and. Ydw.lt._ZERO_) then ! Patankar (1980) trick + fac = -Y(1) / ( dt * Af(1) * Ydw ) + if ( fac.gt.TINY(_ONE_) ) bu(1) = bu(1) + _ONE_/fac + else + du(1) = du(1) + dt * Af(1) * Ydw + end if + case(Dirichlet) + cu(1) = _ZERO_ + bu(1) = Vc(1) + du(1) = Vc(1) * Ydw + wdw = _ONE_ + case default + FATAL 'invalid boundary condition type for lower boundary' + stop 'diffusion.F90' + end select + + if (N .eq. 1) then + if (Bcup.eq.Dirichlet .or. Bcdw.eq.Dirichlet) then +! overwrite all earlier assignments + bu(1) = Vc(1) + du(1) = Vc(1) * ( wup*Yup + wdw*Ydw ) / ( wup + wdw ) + end if + end if + +! relaxation to observed value + if (minval(Taur).lt.1.d10) then + do concurrent ( i = 1 : N ) + bu(i) = bu(i) + dt / Taur(i) * Vc(i) + du(i) = du(i) + dt / Taur(i) * Vc(i) * Yobs(i) + end do + end if + +! solve linear system + call tridiagonal(N, 1, N, Y) + + return + end subroutine diffusion +!EOC + +!----------------------------------------------------------------------- +! Copyright by the GOTM-team under the GNU Public License - www.gnu.org +!-----------------------------------------------------------------------