/*@@ @file apply_dissipation.F90 @date Feb 23 2013 @author Kentaro TAKAMI @desc Original one is "apply_dissipation.F77" which is written by Erik Schnetter. This is rewritten by Fortran 90 style, and "DO" loops are removed. @enddesc @@*/ #include "cctk.h" SUBROUTINE apply_dissipation( var_in, rhs, ni, nj, nk, dx, order, epsdis ) IMPLICIT NONE CCTK_INT, INTENT(IN) :: ni, nj, nk CCTK_REAL, INTENT(INOUT) :: rhs ( 1:ni, 1:nj, 1:nk ) CCTK_REAL, INTENT(IN) :: var_in( 1:ni, 1:nj, 1:nk ) CCTK_REAL, INTENT(IN) :: dx( 1:3 ) CCTK_INT, INTENT(IN) :: order CCTK_REAL, INTENT(IN) :: epsdis( 1:ni, 1:nj, 1:nk ) CCTK_REAL :: idx( 1:3 ) CCTK_INT :: ii, jj, kk !!* In order to avoid random results which is "maybe" caused by !!* variations in the starting address and alignment of the global stack, !!* we have to copy "var_in" to "var". CCTK_REAL :: var( 1:ni, 1:nj, 1:nk ) var = var_in idx( 1:3 ) = 1.0d0 / dx( 1:3 ) IF_ORDER: IF ( order == 1 ) THEN ii = ( ni - 1 ) jj = ( nj - 1 ) kk = ( nk - 1 ) rhs( 2:ii, 2:jj, 2:kk ) = & & rhs( 2:ii, 2:jj, 2:kk ) & & + epsdis( 2:ii, 2:jj, 2:kk ) & & * ( ( var( 1:(ii-1), 2:jj, 2:kk ) & & - 2.0d0 * var( 2:ii, 2:jj, 2:kk ) & & + var( 3:ni, 2:jj, 2:kk ) & & ) * idx(1) & & + ( var( 2:ii, 1:(jj-1), 2:kk ) & & - 2.0d0 * var( 2:ii, 2:jj, 2:kk ) & & + var( 2:ii, 3:nj, 2:kk ) & & ) * idx(2) & & + ( var( 2:ii, 2:jj, 1:(kk-1) ) & & - 2.0d0 * var( 2:ii, 2:jj, 2:kk ) & & + var( 2:ii, 2:jj, 3:nk ) & & ) * idx(3) & & ) ELSE IF ( order == 3 ) THEN IF_ORDER ii = ( ni - 2 ) jj = ( nj - 2 ) kk = ( nk - 2 ) rhs( 3:ii, 3:jj, 3:kk ) = & & rhs( 3:ii, 3:jj, 3:kk ) & & - epsdis( 3:ii, 3:jj, 3:kk ) / 16.0d0 & & * ( ( var( 1:(ii-2), 3:jj, 3:kk ) - 4.0d0 * var( 2:(ii-1), 3:jj, 3:kk ) & & + 6.0d0 * var( 3:ii, 3:jj, 3:kk ) & & - 4.0d0 * var( 4:(ii+1), 3:jj, 3:kk ) + var( 5:ni, 3:jj, 3:kk ) & & ) * idx(1) & & + ( var( 3:ii, 1:(jj-2), 3:kk ) - 4.0d0 * var( 3:ii, 2:(jj-1), 3:kk ) & & + 6.0d0 * var( 3:ii, 3:jj, 3:kk ) & & - 4.0d0 * var( 3:ii, 4:(jj+1), 3:kk ) + var( 3:ii, 5:nj, 3:kk ) & & ) * idx(2) & & + ( var( 3:ii, 3:jj, 1:(kk-2) ) - 4.0d0 * var( 3:ii, 3:jj, 2:(kk-1) ) & & + 6.0d0 * var( 3:ii, 3:jj, 3:kk ) & & - 4.0d0 * var( 3:ii, 3:jj, 4:(kk+1) ) + var( 3:ii, 3:jj, 5:nk ) & & ) * idx(3) & & ) ELSE IF ( order == 5 ) THEN IF_ORDER ii = ( ni - 3 ) jj = ( nj - 3 ) kk = ( nk - 3 ) rhs( 4:ii, 4:jj, 4:kk ) = & & rhs( 4:ii, 4:jj, 4:kk ) & & + epsdis( 4:ii, 4:jj, 4:kk ) / 64.0d0 & & * ( ( var( 1:(ii-3), 4:jj, 4:kk ) - 6.0d0 * var( 2:(ii-2), 4:jj, 4:kk ) & & + 15.0d0 * var( 3:(ii-1), 4:jj, 4:kk ) & & - 20.0d0 * var( 4:ii, 4:jj, 4:kk ) & & + 15.0d0 * var( 5:(ii+1), 4:jj, 4:kk ) & & - 6.0d0 * var( 6:(ii+2), 4:jj, 4:kk ) + var( 7:ni, 4:jj, 4:kk ) & & ) * idx(1) & & + ( var( 4:ii, 1:(jj-3), 4:kk ) - 6.0d0 * var( 4:ii, 2:(jj-2), 4:kk ) & & + 15.0d0 * var( 4:ii, 3:(jj-1), 4:kk ) & & - 20.0d0 * var( 4:ii, 4:jj, 4:kk ) & & + 15.0d0 * var( 4:ii, 5:(jj+1), 4:kk ) & & - 6.0d0 * var( 4:ii, 6:(jj+2), 4:kk ) + var( 4:ii, 7:nj, 4:kk ) & & ) * idx(2) & & + ( var( 4:ii, 4:jj, 1:(kk-3) ) - 6.0d0 * var( 4:ii, 4:jj, 2:(kk-2) ) & & + 15.0d0 * var( 4:ii, 4:jj, 3:(kk-1) ) & & - 20.0d0 * var( 4:ii, 4:jj, 4:kk ) & & + 15.0d0 * var( 4:ii, 4:jj, 5:(kk+1) ) & & - 6.0d0 * var( 4:ii, 4:jj, 6:(kk+2) ) + var( 4:ii, 4:jj, 7:nk ) & & ) * idx(3) & & ) ELSE IF ( order == 7 ) THEN IF_ORDER ii = ( ni - 4 ) jj = ( nj - 4 ) kk = ( nk - 4 ) rhs( 5:ii, 5:jj, 5:kk ) = & & rhs( 5:ii, 5:jj, 5:kk ) & & - epsdis( 5:ii, 5:jj, 5:kk ) / 256.0d0 & & * ( ( var( 1:(ii-4), 5:jj, 5:kk ) - 8.0d0 * var( 2:(ii-3), 5:jj, 5:kk ) & & + 28.0d0 * var( 3:(ii-2), 5:jj, 5:kk ) - 56.0d0 * var( 4:(ii-1), 5:jj, 5:kk ) & & + 70.0d0 * var( 5:ii, 5:jj, 5:kk ) & & - 56.0d0 * var( 6:(ii+1), 5:jj, 5:kk ) + 28.0d0 * var( 7:(ii+2), 5:jj, 5:kk ) & & - 8.0d0 * var( 8:(ii+3), 5:jj, 5:kk ) + var( 9:ni, 5:jj, 5:kk ) & & ) * idx(1) & & + ( var( 5:ii, 1:(jj-4), 5:kk ) - 8.0d0 * var( 5:ii, 2:(jj-3), 5:kk ) & & + 28.0d0 * var( 5:ii, 3:(jj-2), 5:kk ) - 56.0d0 * var( 5:ii, 4:(jj-1), 5:kk ) & & + 70.0d0 * var( 5:ii, 5:jj, 5:kk ) & & - 56.0d0 * var( 5:ii, 6:(jj+1), 5:kk ) + 28.0d0 * var( 5:ii, 7:(jj+2), 5:kk ) & & - 8.0d0 * var( 5:ii, 8:(jj+3), 5:kk ) + var( 5:ii, 9:nj, 5:kk ) & & ) * idx(2) & & + ( var( 5:ii, 5:jj, 1:(kk-4) ) - 8.0d0 * var( 5:ii, 5:jj, 2:(kk-3) ) & & + 28.0d0 * var( 5:ii, 5:jj, 3:(kk-2) ) - 56.0d0 * var( 5:ii, 5:jj, 4:(kk-1) ) & & + 70.0d0 * var( 5:ii, 5:jj, 5:kk ) & & - 56.0d0 * var( 5:ii, 5:jj, 6:(kk+1) ) + 28.0d0 * var( 5:ii, 5:jj, 7:(kk+2) ) & & - 8.0d0 * var( 5:ii, 5:jj, 8:(kk+3) ) + var( 5:ii, 5:jj, 9:nk ) & & ) * idx(3) & & ) ELSE IF ( order == 9 ) THEN IF_ORDER ii = ( ni - 5 ) jj = ( nj - 5 ) kk = ( nk - 5 ) rhs( 6:ii, 6:jj, 6:kk ) = & & rhs( 6:ii, 6:jj, 6:kk ) & & + epsdis( 6:ii, 6:jj, 6:kk ) / 1024.0d0 & & * ( ( var( 1:(ii-5), 6:jj, 6:kk ) - 10.0d0 * var( 2:(ii-4), 6:jj, 6:kk ) & & + 45.0d0 * var( 3:(ii-3), 6:jj, 6:kk ) - 120.0d0 * var( 4:(ii-2), 6:jj, 6:kk ) & & + 210.0d0 * var( 5:(ii-1), 6:jj, 6:kk ) & & - 252.0d0 * var( 6:ii, 6:jj, 6:kk ) & & + 210.0d0 * var( 7:(ii+1), 6:jj, 6:kk ) & & - 120.0d0 * var( 8:(ii+2), 6:jj, 6:kk ) + 45.0d0 * var( 9:(ii+3), 6:jj, 6:kk ) & & - 10.0d0 * var( 10:(ii+4), 6:jj, 6:kk ) + var( 11:ni, 6:jj, 6:kk ) & & ) * idx(1) & & + ( var( 6:ii, 1:(jj-5), 6:kk ) - 10.0d0 * var( 6:ii, 2:(jj-4), 6:kk ) & & + 45.0d0 * var( 6:ii, 3:(jj-3), 6:kk ) - 120.0d0 * var( 6:ii, 4:(jj-2), 6:kk ) & & + 210.0d0 * var( 6:ii, 5:(jj-1), 6:kk ) & & - 252.0d0 * var( 6:ii, 6:jj, 6:kk ) & & + 210.0d0 * var( 6:ii, 7:(jj+1), 6:kk ) & & - 120.0d0 * var( 6:ii, 8:(jj+2), 6:kk ) + 45.0d0 * var( 6:ii, 9:(jj+3), 6:kk ) & & - 10.0d0 * var( 6:ii, 10:(jj+4), 6:kk ) + var( 6:ii, 11:nj, 6:kk ) & & ) * idx(2) & & + ( var( 6:ii, 6:jj, 1:(kk-5) ) - 10.0d0 * var( 6:ii, 6:jj, 2:(kk-4) ) & & + 45.0d0 * var( 6:ii, 6:jj, 3:(kk-3) ) - 120.0d0 * var( 6:ii, 6:jj, 4:(kk-2) ) & & + 210.0d0 * var( 6:ii, 6:jj, 5:(kk-1) ) & & - 252.0d0 * var( 6:ii, 6:jj, 6:kk ) & & + 210.0d0 * var( 6:ii, 6:jj, 7:(kk+1) ) & & - 120.0d0 * var( 6:ii, 6:jj, 8:(kk+2) ) + 45.0d0 * var( 6:ii, 6:jj, 9:(kk+3) ) & & - 10.0d0 * var( 6:ii, 6:jj, 10:(kk+4) ) + var( 6:ii, 6:jj, 11:nk ) & & ) * idx(3) & & ) ELSE IF_ORDER CALL CCTK_WARN( 0, "internal error" ) END IF IF_ORDER RETURN END SUBROUTINE apply_dissipation