diff options
author | diener <diener@f69c4107-0314-4c4f-9ad4-17e986b73f4a> | 2004-10-06 00:18:32 +0000 |
---|---|---|
committer | diener <diener@f69c4107-0314-4c4f-9ad4-17e986b73f4a> | 2004-10-06 00:18:32 +0000 |
commit | 9d639bc123f2c23f5db7c436ec63e21c07e8746b (patch) | |
tree | bab3917f0e11b0246725a7e0642f4f9224dd1979 | |
parent | 46f32e0ccbf4aa49a0de2425aa67338b53288cd2 (diff) |
First attempt at a summation by parts finite differencing thorn.
git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@2 f69c4107-0314-4c4f-9ad4-17e986b73f4a
-rw-r--r-- | README | 9 | ||||
-rw-r--r-- | doc/documentation.tex | 144 | ||||
-rw-r--r-- | interface.ccl | 31 | ||||
-rw-r--r-- | param.ccl | 8 | ||||
-rw-r--r-- | schedule.ccl | 2 | ||||
-rw-r--r-- | src/Derivatives_2_1.F90 | 96 | ||||
-rw-r--r-- | src/Derivatives_4_2.F90 | 160 | ||||
-rw-r--r-- | src/Derivatives_6_3.F90 | 239 | ||||
-rw-r--r-- | src/call_derivs.c | 105 | ||||
-rw-r--r-- | src/make.code.defn | 8 |
10 files changed, 802 insertions, 0 deletions
@@ -0,0 +1,9 @@ +CVS info : $Header$ + +Cactus Code Thorn SummationByParts +Thorn Author(s) : Peter Diener <diener@cct.lsu.edu> +Thorn Maintainer(s) : Peter Diener <diener@cct.lsu.edu> +-------------------------------------------------------------------------- + +Purpose of the thorn: + diff --git a/doc/documentation.tex b/doc/documentation.tex new file mode 100644 index 0000000..0b64bc1 --- /dev/null +++ b/doc/documentation.tex @@ -0,0 +1,144 @@ +% *======================================================================* +% Cactus Thorn template for ThornGuide documentation +% Author: Ian Kelley +% Date: Sun Jun 02, 2002 +% $Header$ +% +% Thorn documentation in the latex file doc/documentation.tex +% will be included in ThornGuides built with the Cactus make system. +% The scripts employed by the make system automatically include +% pages about variables, parameters and scheduling parsed from the +% relevant thorn CCL files. +% +% This template contains guidelines which help to assure that your +% documentation will be correctly added to ThornGuides. More +% information is available in the Cactus UsersGuide. +% +% Guidelines: +% - Do not change anything before the line +% % START CACTUS THORNGUIDE", +% except for filling in the title, author, date, etc. fields. +% - Each of these fields should only be on ONE line. +% - Author names should be separated with a \\ or a comma. +% - You can define your own macros, but they must appear after +% the START CACTUS THORNGUIDE line, and must not redefine standard +% latex commands. +% - To avoid name clashes with other thorns, 'labels', 'citations', +% 'references', and 'image' names should conform to the following +% convention: +% ARRANGEMENT_THORN_LABEL +% For example, an image wave.eps in the arrangement CactusWave and +% thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps +% - Graphics should only be included using the graphicx package. +% More specifically, with the "\includegraphics" command. Do +% not specify any graphic file extensions in your .tex file. This +% will allow us to create a PDF version of the ThornGuide +% via pdflatex. +% - References should be included with the latex "\bibitem" command. +% - Use \begin{abstract}...\end{abstract} instead of \abstract{...} +% - Do not use \appendix, instead include any appendices you need as +% standard sections. +% - For the benefit of our Perl scripts, and for future extensions, +% please use simple latex. +% +% *======================================================================* +% +% Example of including a graphic image: +% \begin{figure}[ht] +% \begin{center} +% \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure} +% \end{center} +% \caption{Illustration of this and that} +% \label{MyArrangement_MyThorn_MyLabel} +% \end{figure} +% +% Example of using a label: +% \label{MyArrangement_MyThorn_MyLabel} +% +% Example of a citation: +% \cite{MyArrangement_MyThorn_Author99} +% +% Example of including a reference +% \bibitem{MyArrangement_MyThorn_Author99} +% {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999), +% 1--16. {\tt http://www.nowhere.com/}} +% +% *======================================================================* + +% If you are using CVS use this line to give version information +% $Header$ + +\documentclass{article} + +% Use the Cactus ThornGuide style file +% (Automatically used from Cactus distribution, if you have a +% thorn without the Cactus Flesh download this from the Cactus +% homepage at www.cactuscode.org) +\usepackage{../../../../doc/latex/cactus} + +\begin{document} + +% The author of the documentation +\author{Peter Diener \textless diener@cct.lsu.edu\textgreater} + +% The title of the document (not necessarily the name of the Thorn) +\title{} + +% the date your document was last changed, if your document is in CVS, +% please use: +% \date{$ $Date$ $} +\date{September 29 2004} + +\maketitle + +% Do not delete next line +% START CACTUS THORNGUIDE + +% Add all definitions used in this documentation here +% \def\mydef etc + +% Add an abstract for this thorn's documentation +\begin{abstract} + +\end{abstract} + +% The following sections are suggestive only. +% Remove them or add your own. + +\section{Introduction} + +\section{Physical System} + +\section{Numerical Implementation} + +\section{Using This Thorn} + +\subsection{Obtaining This Thorn} + +\subsection{Basic Usage} + +\subsection{Special Behaviour} + +\subsection{Interaction With Other Thorns} + +\subsection{Examples} + +\subsection{Support and Feedback} + +\section{History} + +\subsection{Thorn Source Code} + +\subsection{Thorn Documentation} + +\subsection{Acknowledgements} + + +\begin{thebibliography}{9} + +\end{thebibliography} + +% Do not delete next line +% END CACTUS THORNGUIDE + +\end{document} diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..c306db4 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,31 @@ +# Interface definition for thorn SummationByParts +# $Header$ + +implements: SummationByParts +inherits: grid + +#SUBROUTINE Diff_gf ( CCTK_REAL IN ARRAY var, \ +# CCTK_INT IN ni, \ +# CCTK_INT IN nj, \ +# CCTK_INT IN nk, \ +# CCTK_INT IN dir, \ +# CCTK_REAL IN delta, \ +# CCTK_REAL OUT ARRAY dvar ) +#PROVIDES FUNCTION Diff_gf WITH deriv_gf LANGUAGE Fortran + +SUBROUTINE Diff_gf ( CCTK_POINTER IN cctkGH, \ + CCTK_INT IN dir, \ + CCTK_STRING IN var_name, \ + CCTK_STRING IN dvar_name ) +PROVIDES FUNCTION Diff_gf WITH DiffGf LANGUAGE C + +CCTK_INT FUNCTION GetDomainSpecification \ + (CCTK_INT IN size, \ + CCTK_REAL OUT ARRAY physical_min, \ + CCTK_REAL OUT ARRAY physical_max, \ + CCTK_REAL OUT ARRAY interior_min, \ + CCTK_REAL OUT ARRAY interior_max, \ + CCTK_REAL OUT ARRAY exterior_min, \ + CCTK_REAL OUT ARRAY exterior_max, \ + CCTK_REAL OUT ARRAY spacing) +USES FUNCTION GetDomainSpecification diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..ca42e42 --- /dev/null +++ b/param.ccl @@ -0,0 +1,8 @@ +# Parameter definitions for thorn SummationByParts +# $Header$ + +INT order "Order of accuracy" STEERABLE=always +{ + 2:6 :: "" +} 2 + diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..8ebd83a --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,2 @@ +# Schedule definitions for thorn SummationByParts +# $Header$ diff --git a/src/Derivatives_2_1.F90 b/src/Derivatives_2_1.F90 new file mode 100644 index 0000000..564e575 --- /dev/null +++ b/src/Derivatives_2_1.F90 @@ -0,0 +1,96 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine deriv_gf_2_1 ( var, ni, nj, nk, dir, bb, gsize, delta, dvar ) + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_REAL, parameter :: zero = 0.0 + integer, parameter :: wp = kind(zero) + CCTK_INT, intent(IN) :: ni, nj, nk + CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var + CCTK_INT, intent(IN) :: dir + CCTK_INT, intent(IN) :: bb(2) + CCTK_INT, intent(IN) :: gsize + CCTK_REAL, intent(IN) :: delta + CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar + + CCTK_REAL, dimension(1), save :: a + CCTK_REAL, dimension(3,2), save :: q + CCTK_REAL :: idel + + CCTK_INT :: il, ir, jl, jr, kl, kr + + logical, save :: first = .true. + + if ( first ) then + a(1) = 1.0_wp/2.0_wp + + q(1,1) = -1.0_wp; q(2,1) = 1.0_wp; q(3,1) = zero + q(1,2) = -1.0_wp/2.0_wp; q(2,2) = zero; q(3,2) = 1.0_wp/2.0_wp + first = .false. + end if + + idel = 1.0_wp / delta + + direction: select case (dir) + case (0) direction + if ( bb(1) == 0 ) then + il = 1 + gsize + else + dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) ) * idel + dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(3,2) * var(3,:,:) ) * idel + il = 3 + end if + if ( bb(2) == 0 ) then + ir = ni - gsize + else + dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + & + q(3,2) * var(ni-2,:,:) ) * idel + dvar(ni,:,:) = - ( q(1,1) * var(ni,:,:) + & + q(2,1) * var(ni-1,:,:) ) * idel + ir = ni - 2 + end if + dvar(il:ir,:,:) = a(1) * ( var(il+1:ir+1,:,:) - var(il-1:ir-1,:,:) ) * idel + case (1) direction + if ( bb(1) == 0 ) then + jl = 1 + gsize + else + dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) ) * idel + dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(3,2) * var(:,3,:) ) * idel + jl = 3 + end if + if ( bb(2) == 0 ) then + jr = nj - gsize + else + dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + & + q(3,2) * var(:,nj-2,:) ) * idel + dvar(:,nj,:) = - ( q(1,1) * var(:,nj,:) + & + q(2,1) * var(:,nj-1,:) ) * idel + jr = nj - 2 + end if + dvar(:,jl:jr,:) = a(1) * ( var(:,jl+1:jr+1,:) - var(:,jl-1:jr-1,:) ) * idel + case (2) direction + if ( bb(1) == 0 ) then + kl = 1 + gsize + else + dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) ) * idel + dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(3,2) * var(:,:,3) ) * idel + kl = 3 + end if + if ( bb(2) == 0 ) then + kr = nk - gsize + else + dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + & + q(3,2) * var(:,:,nk-2) ) * idel + dvar(:,:,nk) = - ( q(1,1) * var(:,:,nk) + & + q(2,1) * var(:,:,nk-1) ) * idel + kr = nk - 2 + end if + dvar(:,:,kl:kr) = a(1) * ( var(:,:,kl+1:kr+1) - var(:,:,kl-1:kr-1) ) * idel + end select direction +end subroutine deriv_gf_2_1 diff --git a/src/Derivatives_4_2.F90 b/src/Derivatives_4_2.F90 new file mode 100644 index 0000000..ec56613 --- /dev/null +++ b/src/Derivatives_4_2.F90 @@ -0,0 +1,160 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine deriv_gf_4_2 ( var, ni, nj, nk, dir, bb, gsize, delta, dvar ) + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_REAL, parameter :: zero = 0.0 + integer, parameter :: wp = kind(zero) + CCTK_INT, intent(IN) :: ni, nj, nk + CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var + CCTK_INT, intent(IN) :: dir + CCTK_INT, intent(IN) :: bb(2) + CCTK_INT, intent(IN) :: gsize + CCTK_REAL, intent(IN) :: delta + CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar + + CCTK_REAL, dimension(2), save :: a + CCTK_REAL, dimension(6,4), save :: q + CCTK_REAL :: idel + + CCTK_INT :: il, ir, jl, jr, kl, kr + + logical, save :: first = .true. + + if ( first ) then + a(1) = 2.0_wp/3.0_wp; a(2) = -1.0_wp/12.0_wp + + q(1,1) = -24.0_wp/17.0_wp; q(2,1) = 59.0_wp/34.0_wp + q(3,1) = -4.0_wp/17.0_wp; q(4,1) = -3.0_wp/34.0_wp + q(5,1) = zero; q(6,1) = zero + q(1,2) = -1.0_wp/2.0_wp; q(2,2) = zero + q(3,2) = 1.0_wp/2.0_wp; q(4,2) = zero + q(5,2) = zero; q(6,2) = zero + q(1,3) = 4.0_wp/43.0_wp; q(2,3) = -59.0_wp/86.0_wp + q(3,3) = zero; q(4,3) = 59.0_wp/86.0_wp + q(5,3) = -4.0_wp/43.0_wp; q(6,3) = zero + q(1,4) = 3.0_wp/98.0_wp; q(2,4) = zero + q(3,4) = -59.0_wp/98.0_wp; q(4,4) = zero + q(5,4) = 32.0_wp/49.0_wp; q(6,4) = -4.0_wp/49.0_wp + first = .false. + end if + + idel = 1.0_wp / delta + + direction: select case (dir) + case (0) direction + if ( bb(1) == 0 ) then + il = 1 + gsize + else + dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) + & + q(3,1) * var(3,:,:) + q(4,1) * var(4,:,:) ) * idel + dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(3,2) * var(3,:,:) ) * idel + dvar(3,:,:) = ( q(1,3) * var(1,:,:) + q(2,3) * var(2,:,:) + & + q(4,3) * var(4,:,:) + q(5,3) * var(5,:,:) ) * idel + dvar(4,:,:) = ( q(1,4) * var(1,:,:) + q(3,4) * var(3,:,:) + & + q(5,4) * var(5,:,:) + q(6,4) * var(6,:,:) ) * idel + il = 5 + end if + if ( bb(2) == 0 ) then + ir = ni - gsize + else + dvar(ni-3,:,:) = - ( q(1,4) * var(ni,:,:) + & + q(3,4) * var(ni-2,:,:) + & + q(5,4) * var(ni-4,:,:) + & + q(6,4) * var(ni-5,:,:) ) * idel + dvar(ni-2,:,:) = - ( q(1,3) * var(ni,:,:) + & + q(2,3) * var(ni-1,:,:) + & + q(4,3) * var(ni-3,:,:) + & + q(5,3) * var(ni-4,:,:) ) * idel + dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + & + q(3,2) * var(ni-2,:,:) ) * idel + dvar(ni,:,:) = - ( q(1,1) * var(ni,:,:) + & + q(2,1) * var(ni-1,:,:) + & + q(3,1) * var(ni-2,:,:) + & + q(4,1) * var(ni-3,:,:) ) * idel + ir = ni - 4 + end if + dvar(il:ir,:,:) = ( a(1) * ( var(il+1:ir+1,:,:) - & + var(il-1:ir-1,:,:) ) + & + a(2) * ( var(il+2:ir+2,:,:) - & + var(il-2:ir-2,:,:) ) ) * idel + case (1) direction + if ( bb(1) == 0 ) then + jl = 1 + gsize + else + dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + & + q(3,1) * var(:,3,:) + q(4,1) * var(:,4,:) ) * idel + dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(3,2) * var(:,3,:) ) * idel + dvar(:,3,:) = ( q(1,3) * var(:,1,:) + q(2,3) * var(:,2,:) + & + q(4,3) * var(:,4,:) + q(5,3) * var(:,5,:) ) * idel + dvar(:,4,:) = ( q(1,4) * var(:,1,:) + q(3,4) * var(:,3,:) + & + q(5,4) * var(:,5,:) + q(6,4) * var(:,6,:) ) * idel + jl = 5 + end if + if ( bb(2) == 0 ) then + jr = nj - gsize + else + dvar(:,nj-3,:) = - ( q(1,4) * var(:,nj,:) + & + q(3,4) * var(:,nj-2,:) + & + q(5,4) * var(:,nj-4,:) + & + q(6,4) * var(:,nj-5,:) ) * idel + dvar(:,nj-2,:) = - ( q(1,3) * var(:,nj,:) + & + q(2,3) * var(:,nj-1,:) + & + q(4,3) * var(:,nj-3,:) + & + q(5,3) * var(:,nj-4,:) ) * idel + dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + & + q(3,2) * var(:,nj-2,:) ) * idel + dvar(:,nj,:) = - ( q(1,1) * var(:,nj,:) + & + q(2,1) * var(:,nj-1,:) + & + q(3,1) * var(:,nj-2,:) + & + q(4,1) * var(:,nj-3,:) ) * idel + jr = nj - 4 + end if + dvar(:,jl:jr,:) = ( a(1) * ( var(:,jl+1:jr+1,:) - & + var(:,jl-1:jr-1,:) ) + & + a(2) * ( var(:,jl+2:jr+2,:) - & + var(:,jl-2:jr-2,:) ) ) * idel + case (2) direction + if ( bb(1) == 0 ) then + kl = 1 + gsize + else + dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + & + q(3,1) * var(:,:,3) + q(4,1) * var(:,:,4) ) * idel + dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(3,2) * var(:,:,3) ) * idel + dvar(:,:,3) = ( q(1,3) * var(:,:,1) + q(2,3) * var(:,:,2) + & + q(4,3) * var(:,:,4) + q(5,3) * var(:,:,5) ) * idel + dvar(:,:,4) = ( q(1,4) * var(:,:,1) + q(3,4) * var(:,:,3) + & + q(5,4) * var(:,:,5) + q(6,4) * var(:,:,6) ) * idel + kl = 5 + end if + if ( bb(2) == 0 ) then + kr = nk - gsize + else + dvar(:,:,nk-3) = - ( q(1,4) * var(:,:,nk) + & + q(3,4) * var(:,:,nk-2) + & + q(5,4) * var(:,:,nk-4) + & + q(6,4) * var(:,:,nk-5) ) * idel + dvar(:,:,nk-2) = - ( q(1,3) * var(:,:,nk) + & + q(2,3) * var(:,:,nk-1) + & + q(4,3) * var(:,:,nk-3) + & + q(5,3) * var(:,:,nk-4) ) * idel + dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + & + q(3,2) * var(:,:,nk-2) ) * idel + dvar(:,:,nk) = - ( q(1,1) * var(:,:,nk) + & + q(2,1) * var(:,:,nk-1) + & + q(3,1) * var(:,:,nk-2) + & + q(4,1) * var(:,:,nk-3) ) * idel + kr = nk - 4 + end if + dvar(:,:,kl:kr) = ( a(1) * ( var(:,:,kl+1:kr+1) - & + var(:,:,kl-1:kr-1) ) + & + a(2) * ( var(:,:,kl+2:kr+2) - & + var(:,:,kl-2:kr-2) ) ) * idel + end select direction +end subroutine deriv_gf_4_2 diff --git a/src/Derivatives_6_3.F90 b/src/Derivatives_6_3.F90 new file mode 100644 index 0000000..3da882d --- /dev/null +++ b/src/Derivatives_6_3.F90 @@ -0,0 +1,239 @@ +#include "cctk.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + +subroutine deriv_gf_6_3 ( var, ni, nj, nk, dir, bb, gsize, delta, dvar ) + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + CCTK_REAL, parameter :: zero = 0.0 + integer, parameter :: wp = kind(zero) + CCTK_INT, intent(IN) :: ni, nj, nk + CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var + CCTK_INT, intent(IN) :: dir + CCTK_INT, intent(IN) :: bb(2) + CCTK_INT, intent(IN) :: gsize + CCTK_REAL, intent(IN) :: delta + CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar + + CCTK_REAL, dimension(3), save :: a + CCTK_REAL, dimension(9,6), save :: q + CCTK_REAL :: idel + + CCTK_INT :: il, ir, jl, jr, kl, kr + + logical, save :: first = .true. + + if ( first ) then + a(1) = 3.0_wp/4.0_wp; a(2) = -3.0_wp/20.0_wp; a(3) = 1.0_wp/60.0_wp + + q(1,1) = -21600.0_wp/13649.0_wp; q(2,1) = 81763.0_wp/40947.0_wp + q(3,1) = 131.0_wp/27298.0_wp; q(4,1) = -9143.0_wp/13649.0_wp + q(5,1) = 20539.0_wp/81894.0_wp; q(6,1) = zero; + q(7,1) = zero; q(8,1) = zero; q(9,1) = zero + q(1,2) = -81763.0_wp/180195.0_wp; q(2,2) = zero + q(3,2) = 7357.0_wp/36039.0_wp; q(4,2) = 30637.0_wp/72078.0_wp + q(5,2) = -2328.0_wp/12013.0_wp; q(6,2) = 6611.0_wp/360390.0_wp + q(7,2) = zero; q(8,2) = zero; q(9,2) = zero + q(1,3) = -131.0_wp/54220.0_wp; q(2,3) = -7357.0_wp/16266.0_wp + q(3,3) = zero; q(4,3) = 645.0_wp/2711.0_wp + q(5,3) = 11237.0_wp/32532.0_wp; q(6,3) = -3487.0_wp/27110.0_wp + q(7,3) = zero; q(8,3) = zero; q(9,3) = zero + q(1,4) = 9143.0_wp/53590.0_wp; q(2,4) = -30637.0_wp/64308.0_wp + q(3,4) = -645.0_wp/5359.0_wp; q(4,4) = zero + q(5,4) = 13733.0_wp/32154.0_wp; q(6,4) = -67.0_wp/4660.0_wp + q(7,4) = 72.0_wp/5359.0_wp; q(8,4) = zero; q(9,4) = zero; + q(1,5) = -20539.0_wp/236310.0_wp; q(2,5) = 2328.0_wp/7877.0_wp + q(3,5) = -11237.0_wp/47262.0_wp; q(4,5) = -13733.0_wp/23631.0_wp + q(5,5) = zero; q(6,5) = 89387.0_wp/118155.0_wp + q(7,5) = -1296.0_wp/7877.0_wp; q(8,5) = 144.0_wp/7877.0_wp; q(9,5) = zero + q(1,6) = zero; q(2,6) = -6611.0_wp/262806.0_wp + q(3,6) = 3487.0_wp/43801.0_wp; q(4,6) = 1541.0_wp/87602.0_wp + q(5,6) = -89387.0_wp/131403.0_wp; q(6,6) = zero + q(7,6) = 32400.0_wp/43801.0_wp; q(8,6) = -6480.0_wp/43801.0_wp + q(9,6) = 720.0_wp/43801.0_wp + first = .false. + end if + + idel = 1.0_wp / delta + + direction: select case (dir) + case (0) direction + if ( bb(1) == 0 ) then + il = 1 + gsize + else + dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) + & + q(3,1) * var(3,:,:) + q(4,1) * var(4,:,:) + & + q(5,1) * var(5,:,:) ) * idel + dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(3,2) * var(3,:,:) + & + q(4,2) * var(4,:,:) + q(5,2) * var(5,:,:) + & + q(6,2) * var(6,:,:) ) * idel + dvar(3,:,:) = ( q(1,3) * var(1,:,:) + q(2,3) * var(2,:,:) + & + q(4,3) * var(4,:,:) + q(5,3) * var(5,:,:) + & + q(6,3) * var(6,:,:) ) * idel + dvar(4,:,:) = ( q(1,4) * var(1,:,:) + q(2,4) * var(2,:,:) + & + q(3,4) * var(3,:,:) + q(5,4) * var(5,:,:) + & + q(6,4) * var(6,:,:) + q(7,4) * var(7,:,:) ) * idel + dvar(5,:,:) = ( q(1,5) * var(1,:,:) + q(2,5) * var(2,:,:) + & + q(3,5) * var(3,:,:) + q(4,5) * var(4,:,:) + & + q(6,5) * var(6,:,:) + q(7,5) * var(7,:,:) + & + q(8,5) * var(8,:,:) ) * idel + dvar(6,:,:) = ( q(2,6) * var(2,:,:) + q(3,6) * var(3,:,:) + & + q(4,6) * var(4,:,:) + q(5,6) * var(5,:,:) + & + q(7,6) * var(7,:,:) + q(8,6) * var(8,:,:) + & + q(9,6) * var(9,:,:) ) * idel + il = 7 + end if + if ( bb(2) == 0 ) then + ir = ni - gsize + else + dvar(ni,:,:) = - ( q(1,1) * var(ni,:,:) + q(2,1) * var(ni-1,:,:) + & + q(3,1) * var(ni-2,:,:) + q(4,1) * var(ni-3,:,:) + & + q(5,1) * var(ni-4,:,:) ) * idel + dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + q(3,2) * var(ni-2,:,:) + & + q(4,2) * var(ni-3,:,:) + q(5,2) * var(ni-4,:,:) + & + q(6,2) * var(ni-5,:,:) ) * idel + dvar(ni-2,:,:) = - ( q(1,3) * var(ni,:,:) + q(2,3) * var(ni-1,:,:) + & + q(4,3) * var(ni-3,:,:) + q(5,3) * var(ni-4,:,:) + & + q(6,3) * var(ni-5,:,:) ) * idel + dvar(ni-3,:,:) = - ( q(1,4) * var(ni,:,:) + q(2,4) * var(ni-1,:,:) + & + q(3,4) * var(ni-2,:,:) + q(5,4) * var(ni-4,:,:) + & + q(6,4) * var(ni-5,:,:) + & + q(7,4) * var(ni-6,:,:) ) * idel + dvar(ni-4,:,:) = - ( q(1,5) * var(ni,:,:) + q(2,5) * var(ni-1,:,:) + & + q(3,5) * var(ni-2,:,:) + q(4,5) * var(ni-3,:,:) + & + q(6,5) * var(ni-5,:,:) + q(7,5) * var(ni-6,:,:) + & + q(8,5) * var(ni-7,:,:) ) * idel + dvar(ni-5,:,:) = - ( q(2,6) * var(ni-1,:,:) + q(3,6) * var(ni-2,:,:) + & + q(4,6) * var(ni-3,:,:) + q(5,6) * var(ni-4,:,:) + & + q(7,6) * var(ni-6,:,:) + q(8,6) * var(ni-7,:,:) + & + q(9,6) * var(ni-8,:,:) ) * idel + ir = ni - 6 + end if + dvar(il:ir,:,:) = ( a(1) * ( var(il+1:ir+1,:,:) - & + var(il-1:ir-1,:,:) ) + & + a(2) * ( var(il+2:ir+2,:,:) - & + var(il-2:ir-2,:,:) ) + & + a(3) * ( var(il+3:ir+3,:,:) - & + var(il-3:ir-3,:,:) ) ) * idel + case (1) direction + if ( bb(1) == 0 ) then + jl = 1 + gsize + else + dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + & + q(3,1) * var(:,3,:) + q(4,1) * var(:,4,:) + & + q(5,1) * var(:,5,:) ) * idel + dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(3,2) * var(:,3,:) + & + q(4,2) * var(:,4,:) + q(5,2) * var(:,5,:) + & + q(6,2) * var(:,6,:) ) * idel + dvar(:,3,:) = ( q(1,3) * var(:,1,:) + q(2,3) * var(:,2,:) + & + q(4,3) * var(:,4,:) + q(5,3) * var(:,5,:) + & + q(6,3) * var(:,6,:) ) * idel + dvar(:,4,:) = ( q(1,4) * var(:,1,:) + q(2,4) * var(:,2,:) + & + q(3,4) * var(:,3,:) + q(5,4) * var(:,5,:) + & + q(6,4) * var(:,6,:) + q(7,4) * var(:,7,:) ) * idel + dvar(:,5,:) = ( q(1,5) * var(:,1,:) + q(2,5) * var(:,2,:) + & + q(3,5) * var(:,3,:) + q(4,5) * var(:,4,:) + & + q(6,5) * var(:,6,:) + q(7,5) * var(:,7,:) + & + q(8,5) * var(:,8,:) ) * idel + dvar(:,6,:) = ( q(2,6) * var(:,2,:) + q(3,6) * var(:,3,:) + & + q(4,6) * var(:,4,:) + q(5,6) * var(:,5,:) + & + q(7,6) * var(:,7,:) + q(8,6) * var(:,8,:) + & + q(9,6) * var(:,9,:) ) * idel + jl = 7 + end if + if ( bb(2) == 0 ) then + jr = nj - gsize + else + dvar(:,nj,:) = - ( q(1,1) * var(:,nj,:) + q(2,1) * var(:,nj-1,:) + & + q(3,1) * var(:,nj-2,:) + q(4,1) * var(:,nj-3,:) + & + q(5,1) * var(:,nj-4,:) ) * idel + dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + q(3,2) * var(:,nj-2,:) + & + q(4,2) * var(:,nj-3,:) + q(5,2) * var(:,nj-4,:) + & + q(6,2) * var(:,nj-5,:) ) * idel + dvar(:,nj-2,:) = - ( q(1,3) * var(:,nj,:) + q(2,3) * var(:,nj-1,:) + & + q(4,3) * var(:,nj-3,:) + q(5,3) * var(:,nj-4,:) + & + q(6,3) * var(:,nj-5,:) ) * idel + dvar(:,nj-3,:) = - ( q(1,4) * var(:,nj,:) + q(2,4) * var(:,nj-1,:) + & + q(3,4) * var(:,nj-2,:) + q(5,4) * var(:,nj-4,:) + & + q(6,4) * var(:,nj-5,:) + & + q(7,4) * var(:,nj-6,:) ) * idel + dvar(:,nj-4,:) = - ( q(1,5) * var(:,nj,:) + q(2,5) * var(:,nj-1,:) + & + q(3,5) * var(:,nj-2,:) + q(4,5) * var(:,nj-3,:) + & + q(6,5) * var(:,nj-5,:) + q(7,5) * var(:,nj-6,:) + & + q(8,5) * var(:,nj-7,:) ) * idel + dvar(:,nj-5,:) = - ( q(2,6) * var(:,nj-1,:) + q(3,6) * var(:,nj-2,:) + & + q(4,6) * var(:,nj-3,:) + q(5,6) * var(:,nj-4,:) + & + q(7,6) * var(:,nj-6,:) + q(8,6) * var(:,nj-7,:) + & + q(9,6) * var(:,nj-8,:) ) * idel + jr = nj - 6 + end if + dvar(:,jl:jr,:) = ( a(1) * ( var(:,jl+1:jr+1,:) - & + var(:,jl-1:jr-1,:) ) + & + a(2) * ( var(:,jl+2:jr+2,:) - & + var(:,jl-2:jr-2,:) ) + & + a(3) * ( var(:,jl+3:jr+3,:) - & + var(:,jl-3:jr-3,:) ) ) * idel + case (2) direction + if ( bb(1) == 0 ) then + kl = 1 + gsize + else + dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + & + q(3,1) * var(:,:,3) + q(4,1) * var(:,:,4) + & + q(5,1) * var(:,:,5) ) * idel + dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(3,2) * var(:,:,3) + & + q(4,2) * var(:,:,4) + q(5,2) * var(:,:,5) + & + q(6,2) * var(:,:,6) ) * idel + dvar(:,:,3) = ( q(1,3) * var(:,:,1) + q(2,3) * var(:,:,2) + & + q(4,3) * var(:,:,4) + q(5,3) * var(:,:,5) + & + q(6,3) * var(:,:,6) ) * idel + dvar(:,:,4) = ( q(1,4) * var(:,:,1) + q(2,4) * var(:,:,2) + & + q(3,4) * var(:,:,3) + q(5,4) * var(:,:,5) + & + q(6,4) * var(:,:,6) + q(7,4) * var(:,:,7) ) * idel + dvar(:,:,5) = ( q(1,5) * var(:,:,1) + q(2,5) * var(:,:,2) + & + q(3,5) * var(:,:,3) + q(4,5) * var(:,:,4) + & + q(6,5) * var(:,:,6) + q(7,5) * var(:,:,7) + & + q(8,5) * var(:,:,8) ) * idel + dvar(:,:,6) = ( q(2,6) * var(:,:,2) + q(3,6) * var(:,:,3) + & + q(4,6) * var(:,:,4) + q(5,6) * var(:,:,5) + & + q(7,6) * var(:,:,7) + q(8,6) * var(:,:,8) + & + q(9,6) * var(:,:,9) ) * idel + kl = 7 + end if + if ( bb(2) == 0 ) then + kr = nk - gsize + else + dvar(:,:,nk) = - ( q(1,1) * var(:,:,nk) + q(2,1) * var(:,:,nk-1) + & + q(3,1) * var(:,:,nk-2) + q(4,1) * var(:,:,nk-3) + & + q(5,1) * var(:,:,nk-4) ) * idel + dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + q(3,2) * var(:,:,nk-2) + & + q(4,2) * var(:,:,nk-3) + q(5,2) * var(:,:,nk-4) + & + q(6,2) * var(:,:,nk-5) ) * idel + dvar(:,:,nk-2) = - ( q(1,3) * var(:,:,nk) + q(2,3) * var(:,:,nk-1) + & + q(4,3) * var(:,:,nk-3) + q(5,3) * var(:,:,nk-4) + & + q(6,3) * var(:,:,nk-5) ) * idel + dvar(:,:,nk-3) = - ( q(1,4) * var(:,:,nk) + q(2,4) * var(:,:,nk-1) + & + q(3,4) * var(:,:,nk-2) + q(5,4) * var(:,:,nk-4) + & + q(6,4) * var(:,:,nk-5) + & + q(7,4) * var(:,:,nk-6) ) * idel + dvar(:,:,nk-4) = - ( q(1,5) * var(:,:,nk) + q(2,5) * var(:,:,nk-1) + & + q(3,5) * var(:,:,nk-2) + q(4,5) * var(:,:,nk-3) + & + q(6,5) * var(:,:,nk-5) + q(7,5) * var(:,:,nk-6) + & + q(8,5) * var(:,:,nk-7) ) * idel + dvar(:,:,nk-5) = - ( q(2,6) * var(:,:,nk-1) + q(3,6) * var(:,:,nk-2) + & + q(4,6) * var(:,:,nk-3) + q(5,6) * var(:,:,nk-4) + & + q(7,6) * var(:,:,nk-6) + q(8,6) * var(:,:,nk-7) + & + q(9,6) * var(:,:,nk-8) ) * idel + kr = nk - 6 + end if + dvar(:,:,kl:kr) = ( a(1) * ( var(:,:,kl+1:kr+1) - & + var(:,:,kl-1:kr-1) ) + & + a(2) * ( var(:,:,kl+2:kr+2) - & + var(:,:,kl-2:kr-2) ) + & + a(3) * ( var(:,:,kl+3:kr+3) - & + var(:,:,kl-3:kr-3) ) ) * idel + end select direction +end subroutine deriv_gf_6_3 diff --git a/src/call_derivs.c b/src/call_derivs.c new file mode 100644 index 0000000..c696e4e --- /dev/null +++ b/src/call_derivs.c @@ -0,0 +1,105 @@ +#include "cctk.h" +#include "cctk_Parameters.h" + +#include <assert.h> + +void DiffGf ( const CCTK_POINTER cctkGH, const CCTK_INT dir, + const char *var_name, const char *dvar_name ) +{ + CCTK_REAL *var, *dvar; + CCTK_INT ni, nj, nk, gsize, ic; + CCTK_REAL delta; + CCTK_REAL phys_min[3], phys_max[3]; + CCTK_REAL int_min[3], int_max[3]; + CCTK_REAL ext_min[3], ext_max[3]; + CCTK_REAL spacing[3]; + CCTK_INT ierr; + CCTK_INT lsh[3], bbox[6], bb[2], nghostzones[3]; + void CCTK_FCALL CCTK_FNAME(deriv_gf_2_1)(CCTK_REAL *var, const CCTK_INT *ni, + const CCTK_INT *nj, const CCTK_INT *nk, + const CCTK_INT *dir, + const CCTK_INT *bb, + const CCTK_INT *gsize, + const CCTK_REAL *delta, + CCTK_REAL *dvar); + void CCTK_FCALL CCTK_FNAME(deriv_gf_4_2)(CCTK_REAL *var, const CCTK_INT *ni, + const CCTK_INT *nj, const CCTK_INT *nk, + const CCTK_INT *dir, + const CCTK_INT *bb, + const CCTK_INT *gsize, + const CCTK_REAL *delta, + CCTK_REAL *dvar); + void CCTK_FCALL CCTK_FNAME(deriv_gf_6_3)(CCTK_REAL *var, const CCTK_INT *ni, + const CCTK_INT *nj, const CCTK_INT *nk, + const CCTK_INT *dir, + const CCTK_INT *bb, + const CCTK_INT *gsize, + const CCTK_REAL *delta, + CCTK_REAL *dvar); + DECLARE_CCTK_PARAMETERS + + + if (CCTK_GrouplshVN(cctkGH, 3, lsh, var_name) < 0) + { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error getting local size of grid for variable '%s'",var_name); + } + if (CCTK_GroupbboxVN(cctkGH, 6, bbox, var_name) < 0) + { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error getting bounding box for variable '%s'",var_name); + } + if (CCTK_GroupnghostzonesVN(cctkGH, 3, nghostzones, var_name) < 0) + { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error getting number of ghostzones for variable '%s'",var_name); + } + + ni = lsh[0]; nj = lsh[1]; nk = lsh[2]; + + ierr = GetDomainSpecification ( 3, phys_min, phys_max, int_min, int_max, + ext_min, ext_max, spacing ); + + switch(dir) { + case 0: { + delta = spacing[0]; + bb[0] = bbox[0]; bb[1] = bbox[1]; + gsize = nghostzones[0]; + break; + } + case 1: { + delta = spacing[1]; + bb[0] = bbox[2]; bb[1] = bbox[3]; + gsize = nghostzones[1]; + break; + } + case 2: { + delta = spacing[2]; + bb[0] = bbox[4]; bb[1] = bbox[5]; + gsize = nghostzones[2]; + break; + } + default: + assert(0); + } + + var = (CCTK_REAL *)(CCTK_VarDataPtr(cctkGH,0,var_name)); + dvar = (CCTK_REAL *)(CCTK_VarDataPtr(cctkGH,0,dvar_name)); + + switch(order) { + case 2: { + CCTK_FNAME(deriv_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar); + break; + } + case 4: { + CCTK_FNAME(deriv_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar); + break; + } + case 6: { + CCTK_FNAME(deriv_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar); + break; + } + default: + assert(0); + } +} diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..eee5124 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,8 @@ +# Main make.code.defn file for thorn SummationByParts +# $Header$ + +# Source files in this directory +SRCS = Derivatives_2_1.F90 Derivatives_4_2.F90 Derivatives_6_3.F90 call_derivs.c + +# Subdirectories containing source files +SUBDIRS = |