Commits July 2013

commits@lists.einsteintoolkit.org

7 participants
90 discussions

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 216)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:12 PM Modified: /trunk/ interface.ccl Log: GRHydro_InitData: correct GRHydro aliased function interfaces and usage From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +52 -51 =================================================================== --- trunk/interface.ccl 2013-07-06 18:12:15 UTC (rev 215) +++ trunk/interface.ccl 2013-07-06 18:12:54 UTC (rev 216) @@ -83,89 +83,90 @@ CCTK_REAL INOUT temperature, CCTK_REAL IN Y_e) -SUBROUTINE Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL OUT dens, \ - CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ - CCTK_REAL OUT sz, CCTK_REAL OUT tau, \ - CCTK_REAL OUT rho, CCTK_REAL OUT velx, \ - CCTK_REAL OUT vely, CCTK_REAL OUT velz, \ - CCTK_REAL OUT epsilon, CCTK_REAL OUT press, \ - CCTK_REAL OUT w_lorentz, CCTK_REAL IN uxx, \ - CCTK_REAL IN uxy, CCTK_REAL IN uxz, CCTK_REAL IN uyy, \ - CCTK_REAL IN uyz, CCTK_REAL IN uzz, CCTK_REAL IN det, \ - CCTK_REAL IN x, CCTK_REAL IN y, CCTK_REAL IN z, \ - CCTK_REAL IN r, CCTK_REAL IN rho_min, \ - CCTK_INT IN GRHydro_reflevel, CCTK_REAL IN GRHydro_C2P_failed) +SUBROUTINE Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL INOUT dens, \ + CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, \ + CCTK_REAL INOUT sz, CCTK_REAL OUT tau, \ + CCTK_REAL INOUT rho, CCTK_REAL INOUT velx, \ + CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ + CCTK_REAL INOUT epsilon, CCTK_REAL INOUT press, \ + CCTK_REAL INOUT w_lorentz, CCTK_REAL IN uxx, \ + CCTK_REAL IN uxy, CCTK_REAL IN uxz, \ + CCTK_REAL IN uyy, CCTK_REAL IN uyz, \ + CCTK_REAL IN uzz, CCTK_REAL IN det, \ + CCTK_REAL IN x, CCTK_REAL IN y, CCTK_REAL IN z, \ + CCTK_REAL IN r, CCTK_REAL IN rho_min, \ + CCTK_INT IN GRHydro_reflevel, CCTK_REAL OUT GRHydro_C2P_failed) -void FUNCTION Con2PrimGen(CCTK_INT INOUT handle, CCTK_REAL INOUT dens, \ +void FUNCTION Con2PrimGen(CCTK_INT IN handle, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ CCTK_REAL INOUT tau, CCTK_REAL INOUT rho, \ CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ CCTK_REAL INOUT w_lorentz, \ - CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ - CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ - CCTK_REAL INOUT det, \ - CCTK_REAL INOUT x, CCTK_REAL INOUT y, CCTK_REAL INOUT z, \ - CCTK_REAL INOUT r, \ + CCTK_REAL IN uxx, CCTK_REAL IN uxy, CCTK_REAL IN uxz, \ + CCTK_REAL IN uyy, CCTK_REAL IN uyz, CCTK_REAL IN uzz, \ + CCTK_REAL IN det, \ + CCTK_REAL IN x, CCTK_REAL IN y, CCTK_REAL IN z, \ + CCTK_REAL IN r, \ CCTK_INT OUT epsnegative, \ - CCTK_REAL INOUT GRHydro_rho_min, CCTK_REAL INOUT pmin, \ - CCTK_REAL INOUT epsmin, \ - CCTK_INT INOUT GRHydro_reflevel, \ + CCTK_REAL IN GRHydro_rho_min, CCTK_REAL IN pmin, \ + CCTK_REAL IN epsmin, \ + CCTK_INT IN GRHydro_reflevel, \ CCTK_REAL OUT retval) -void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ +void FUNCTION Con2PrimGenM(CCTK_INT IN handle, CCTK_INT IN keytemp, CCTK_REAL IN prec,CCTK_REAL IN gamma_eos, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ - CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ - CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ + CCTK_REAL INOUT tau, CCTK_REAL IN Bconsx, CCTK_REAL IN Bconsy, CCTK_REAL IN Bconsz, \ + CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ - CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ - CCTK_REAL INOUT Bvecsq, \ + CCTK_REAL OUT Bvecx, CCTK_REAL OUT Bvecy, CCTK_REAL OUT Bvecz, \ + CCTK_REAL OUT bvecsq, \ CCTK_REAL INOUT w_lorentz, \ - CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ - CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ - CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ - CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ - CCTK_REAL INOUT det, \ + CCTK_REAL IN gxx, CCTK_REAL IN gxy, CCTK_REAL IN gxz, \ + CCTK_REAL IN gyy, CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ + CCTK_REAL IN uxx, CCTK_REAL IN uxy, CCTK_REAL IN uxz, \ + CCTK_REAL IN uyy, CCTK_REAL IN uyz, CCTK_REAL IN uzz, \ + CCTK_REAL IN det, \ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) -void FUNCTION Con2PrimGenMee(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, \ - CCTK_REAL INOUT prec, CCTK_REAL INOUT gamma_eos, \ +void FUNCTION Con2PrimGenMee(CCTK_INT IN handle, CCTK_INT IN keytemp, \ + CCTK_REAL IN prec, CCTK_REAL IN gamma_eos, \ CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ CCTK_REAL INOUT tau, \ - CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ + CCTK_REAL IN Bconsx, CCTK_REAL IN Bconsy, CCTK_REAL IN Bconsz, \ CCTK_REAL INOUT entropycons, \ CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ - CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ - CCTK_REAL INOUT Bvecsq, \ + CCTK_REAL OUT Bvecx, CCTK_REAL OUT Bvecy, CCTK_REAL OUT Bvecz, \ + CCTK_REAL OUT Bvecsq, \ CCTK_REAL INOUT w_lorentz, \ - CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ - CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ - CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ - CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ - CCTK_REAL INOUT det, \ + CCTK_REAL IN gxx, CCTK_REAL IN gxy, CCTK_REAL IN gxz, \ + CCTK_REAL IN gyy, CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ + CCTK_REAL IN uxx, CCTK_REAL IN uxy, CCTK_REAL IN uxz, \ + CCTK_REAL IN uyy, CCTK_REAL IN uyz, CCTK_REAL IN uzz, \ + CCTK_REAL IN det, \ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) -void FUNCTION Con2PrimPolyM(CCTK_INT INOUT handle, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ +void FUNCTION Con2PrimPolyM(CCTK_INT IN handle, CCTK_REAL IN gamma_eos, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ - CCTK_REAL INOUT sc, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ - CCTK_REAL INOUT rho, \ + CCTK_REAL INOUT sc, CCTK_REAL IN Bconsx, CCTK_REAL IN Bconsy, CCTK_REAL IN Bconsz, \ + CCTK_REAL INOUT rho, \ CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ - CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ - CCTK_REAL INOUT Bvecsq, \ + CCTK_REAL OUT Bvecx, CCTK_REAL OUT Bvecy, CCTK_REAL OUT Bvecz, \ + CCTK_REAL OUT Bvecsq, \ CCTK_REAL INOUT w_lorentz, \ - CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ - CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ - CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ - CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ - CCTK_REAL INOUT det, \ + CCTK_REAL IN gxx, CCTK_REAL IN gxy, CCTK_REAL IN gxz, \ + CCTK_REAL IN gyy, CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ + CCTK_REAL IN uxx, CCTK_REAL IN uxy, CCTK_REAL IN uxz, \ + CCTK_REAL IN uyy, CCTK_REAL IN uyz, CCTK_REAL IN uzz, \ + CCTK_REAL IN det, \ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval)

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 215)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:12 PM Modified: /trunk/ schedule.ccl /trunk/src/ GRHydro_PoloidalMagFieldM.F90 Log: GRHydro_InitData: Reschedule to make things work with multipatch. Spit out a compile-time warning since poloidal ID B-field routine is not yet adapted to MP. From: Christian Reisswig <reisswig(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PoloidalMagFieldM.F90 Delta lines: +9 -1 =================================================================== --- trunk/src/GRHydro_PoloidalMagFieldM.F90 2013-05-27 22:32:37 UTC (rev 214) +++ trunk/src/GRHydro_PoloidalMagFieldM.F90 2013-07-06 18:12:15 UTC (rev 215) @@ -85,7 +85,8 @@ write(*,*)'GRHydro_InitData: Setting up initial poloidal magnetic field' - +! Initialize to zero + Bvec = 0.0d0 do i=2,nx-1 do j=2,ny-1 @@ -106,6 +107,7 @@ !! write(*,*)'delPcut = ', delPcut !! write(*,*)'maxP_Pcut = ', maxP_Pcut +#warning "This algorithm does only work on Cartesian grids!!" rho_dx = 0.5d0*(rho(i+1,j,k)-rho(i-1,j,k))/dx rho_dy = 0.5d0*(rho(i,j+1,k)-rho(i,j-1,k))/dy rho_dz = 0.5d0*(rho(i,j,k+1)-rho(i,j,k-1))/dz @@ -183,6 +185,12 @@ if (clean_divergence .ne. 0) then psidcrhs =0.0 endif + + !Bvec = 0 + !lBvec = 0 + !Bcons = 0 + !Avec = 0 + return end subroutine GRHydro_PoloidalMagFieldM Directory: /trunk/ ================== File [modified]: schedule.ccl Delta lines: +4 -2 =================================================================== --- trunk/schedule.ccl 2013-05-27 22:32:37 UTC (rev 214) +++ trunk/schedule.ccl 2013-07-06 18:12:15 UTC (rev 215) @@ -270,11 +270,13 @@ if (CCTK_EQUALS(initial_Bvec, "poloidalmagfield") || CCTK_EQUALS(initial_Avec, "poloidalmagfield")) { # SCHEDULE GRHydro_PoloidalMagFieldM AT CCTK_INITIAL AFTER IN HydroBase_Initial AFTER rnsid_init AFTER TOV_Initial_Data after CCCC_StarMapper_InitialData - SCHEDULE GRHydro_PoloidalMagFieldM AT CCTK_POSTINITIAL +# SCHEDULE GRHydro_PoloidalMagFieldM AT CCTK_POSTINITIAL + SCHEDULE GRHydro_PoloidalMagFieldM AT CCTK_INITIAL AFTER HydroBase_Initial BEFORE GRHydroTransformPrimToLocalBasis { LANGUAGE: Fortran } "Set up a poloidal magnetic field. It expects the other fluid variables already to be set, as for example in the TOV solution" - SCHEDULE group HydroBase_Boundaries IN HydroBase_Initial AFTER GRHydro_PoloidalMagFieldM +# SCHEDULE group HydroBase_Boundaries IN HydroBase_Initial AFTER GRHydro_PoloidalMagFieldM + SCHEDULE group HydroBase_Boundaries IN CCTK_INITIAL AFTER GRHydro_PoloidalMagFieldM BEFORE GRHydroTransformPrimToLocalBasis { } "Call boundary conditions after magnetic field initial data setup" }

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 560)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:11 PM Added: /trunk/src/ GRHydro_PPM.c, GRHydro_PPM.h, GRHydro_PPMReconstruct_drv_opt.F90 Log: GRHydro: reimplement PPM in C * add prototype of a re-implementation of PPM in C. So far we can do old PPM without and with temperature/Y_e reconstruction. * also add a version of the reconstruction driver that is optimized for the use of the new PPM routine. * no change to standard code version; one needs to manually replace source files to get the new stuff. * new version is factor 3 faster than old version (based on preliminary measurements). From: Christian Ott <cott(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [added]: GRHydro_PPM.c Delta lines: +437 -0 =================================================================== --- trunk/src/GRHydro_PPM.c (rev 0) +++ trunk/src/GRHydro_PPM.c 2013-07-06 18:11:31 UTC (rev 560) @@ -0,0 +1,437 @@ +#include <assert.h> +#include <math.h> +#include <stdio.h> +#include <stdlib.h> + +#if 1 +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#endif + +#include "GRHydro_PPM.h" + + +void SimplePPM_1d_c(int* restrict nx, + double* restrict dx, + double* restrict rho, + double* restrict velx, + double* restrict vely, + double* restrict velz, + double* restrict eps, + double* restrict press, + double* restrict rhominus, + double* restrict velxminus, + double* restrict velyminus, + double* restrict velzminus, + double* restrict epsminus, + double* restrict rhoplus, + double* restrict velxplus, + double* restrict velyplus, + double* restrict velzplus, + double* restrict epsplus) + +{ + double drho[*nx], dpress[*nx], deps[*nx], d2rho[*nx]; + double dmrho[*nx], dmpress[*nx], dmeps[*nx]; + double dvelx[*nx], dvely[*nx], dvelz[*nx]; + double dmvelx[*nx], dmvely[*nx], dmvelz[*nx]; + double tilde_flatten[*nx]; + + const double onesixth = 1.0/6.0; + + for(int i=1;i < *nx-1;i++) { + drho[i] = 0.5 * (rho[i+1]-rho[i-1]); + d2rho[i] = rho[i+1] - 2.0 * rho[i] + rho[i-1]; + dpress[i] = (press[i+1]-press[i-1]); + deps[i] = 0.5 * (eps[i+1]-eps[i-1]); + dvelx[i] = 0.5 * (velx[i+1]-velx[i-1]); + dvely[i] = 0.5 * (vely[i+1]-vely[i-1]); + dvelz[i] = 0.5 * (velz[i+1]-velz[i-1]); + } + + for(int i=1;i<*nx-1;i++) { + steep(rho,drho,dmrho,i); + steep(eps,deps,dmeps,i); + steep(velx,dvelx,dmvelx,i); + steep(vely,dvely,dmvely,i); + steep(velz,dvelz,dmvelz,i); + } + + for(int i=1;i<*nx-1;i++) { + rhoplus[i] = 0.5 * (rho[i] + rho[i+1]) + + (dmrho[i] - dmrho[i+1]) * onesixth; + epsplus[i] = 0.5 * (eps[i] + eps[i+1]) + + (dmeps[i] - dmeps[i+1]) * onesixth; + velxplus[i] = 0.5 * (velx[i] + velx[i+1]) + + (dmvelx[i] - dmvelx[i+1]) * onesixth; + velyplus[i] = 0.5 * (vely[i] + vely[i+1]) + + (dmvely[i] - dmvely[i+1]) * onesixth; + velzplus[i] = 0.5 * (velz[i] + velz[i+1]) + + (dmvelz[i] - dmvelz[i+1]) * onesixth; + } + + for(int i=1;i<*nx-2;i++) { + rhoplus[i] = approx_at_cell_interface(rho,i); + epsplus[i] = approx_at_cell_interface(eps,i); + velxplus[i] = approx_at_cell_interface(velx,i); + velyplus[i] = approx_at_cell_interface(vely,i); + velzplus[i] = approx_at_cell_interface(velz,i); + } + + for(int i=1;i<*nx-2;i++) { + rhominus[i+1] = rhoplus[i]; + epsminus[i+1] = epsplus[i]; + velxminus[i+1] = velxplus[i]; + velyminus[i+1] = velyplus[i]; + velzminus[i+1] = velzplus[i]; + } + + const double ppm_epsilon_shock = 0.01; + const double ppm_eta1 = 20.0; + const double ppm_eta2 = 0.05; + const double ppm_k0 = 0.2; + const double ppm_epsilon = 0.33; + const double ppm_omega1 = 0.75; + const double ppm_omega2 = 10.0; + const double ppm_small = 1.0e-7; + + // steepening + for(int i=2;i<*nx-2;i++) { + double etatilde = 0.0; + if ( d2rho[i+1]*d2rho[i-1] < 0.0 + && ( fabs(rho[i+1]-rho[i-1]) - ppm_epsilon_shock + * MIN(fabs(rho[i+1]), + fabs(rho[i-1])) > 0.0) ) + { + etatilde = (rho[i-2] - rho[i+2] + 4.0 * drho[i]) / (drho[i] * 12.0); + } + double eta = MAX(0.0, MIN(1.0, ppm_eta1 * (etatilde - ppm_eta2))); + if (ppm_k0 * fabs(drho[i]) * MIN(press[i-1],press[i+1]) + < fabs(dpress[i]) * MIN(rho[i-1], rho[i+1])) + { + eta = 0.0; + } + rhominus[i] = rhominus[i] * (1.0 - eta) + + (rho[i-1] + 0.5 * dmrho[i-1]) * eta; + rhoplus[i] = rhoplus[i] * (1.0 - eta) + + (rho[i+1] - 0.5 * dmrho[i+1]) * eta; + } + + // flattening + for(int i=2;i<*nx-2;i++) { + double dpress2 = press[i+2] - press[i-2]; + double dvel = velx[i+1] - velx[i-1]; + double w=0.0; + if ( (fabs(dpress[i]) > ppm_epsilon * MIN(press[i-1],press[i+1])) + && (dvel < 0.0) ) + { + w = 1.0; + } + if (fabs(dpress2) < ppm_small) + { + tilde_flatten[i] = 1.0; + } + else + { + tilde_flatten[i] = MAX(0.0, 1.0 - w * MAX(0.0, + ppm_omega2 * (dpress[i] / dpress2 - ppm_omega1))); + } + } + + for(int i=2;i<*nx-2;i++) { + double flatten = tilde_flatten[i]; + rhoplus[i] = flatten * rhoplus[i] + (1.0 - flatten) * rho[i]; + rhominus[i] = flatten * rhominus[i] + (1.0 - flatten) * rho[i]; + epsplus[i] = flatten * epsplus[i] + (1.0 - flatten) * eps[i]; + epsminus[i] = flatten * epsminus[i] + (1.0 - flatten) * eps[i]; + velxplus[i] = flatten * velxplus[i] + (1.0 - flatten) * velx[i]; + velxminus[i] = flatten * velxminus[i] + (1.0 - flatten) * velx[i]; + velyplus[i] = flatten * velyplus[i] + (1.0 - flatten) * vely[i]; + velyminus[i] = flatten * velyminus[i] + (1.0 - flatten) * vely[i]; + velzplus[i] = flatten * velzplus[i] + (1.0 - flatten) * velz[i]; + velzminus[i] = flatten * velzminus[i] + (1.0 - flatten) * velz[i]; + } // flattening + + for(int i=2;i<*nx-2;i++) { + monotonize(rhominus,rho,rhoplus,i); + monotonize(epsminus,eps,epsplus,i); + monotonize(velxminus,velx,velxplus,i); + monotonize(velyminus,vely,velyplus,i); + monotonize(velzminus,velz,velzplus,i); + } + + return; +} + +CCTK_FCALL void CCTK_FNAME(PPMc)(int* restrict nx, + double* restrict dx, + double* restrict rho, + double* restrict velx, + double* restrict vely, + double* restrict velz, + double* restrict eps, + double* restrict press, + double* restrict rhominus, + double* restrict velxminus, + double* restrict velyminus, + double* restrict velzminus, + double* restrict epsminus, + double* restrict rhoplus, + double* restrict velxplus, + double* restrict velyplus, + double* restrict velzplus, + double* restrict epsplus) { + +void SimplePPM_1d_c(nx, dx, + rho, + velx, + vely, + velz, + eps, + press, + rhominus, + velxminus, + velyminus, + velzminus, + epsminus, + rhoplus, + velxplus, + velyplus, + velzplus, + epsplus); + +} + + +void SimplePPM_1d_tyc(int* restrict nx, + double* restrict dx, + double* restrict rho, + double* restrict velx, + double* restrict vely, + double* restrict velz, + double* restrict eps, + double* restrict temp, + double* restrict ye, + double* restrict press, + double* restrict rhominus, + double* restrict velxminus, + double* restrict velyminus, + double* restrict velzminus, + double* restrict epsminus, + double * restrict tempminus, + double * restrict yeminus, + double* restrict rhoplus, + double* restrict velxplus, + double* restrict velyplus, + double* restrict velzplus, + double* restrict epsplus, + double* restrict tempplus, + double* restrict yeplus) + +{ + double drho[*nx], dpress[*nx], deps[*nx], d2rho[*nx]; + double dtemp[*nx], dye[*nx]; + double dmrho[*nx], dmpress[*nx], dmeps[*nx]; + double dmtemp[*nx], dmye[*nx]; + double dvelx[*nx], dvely[*nx], dvelz[*nx]; + double dmvelx[*nx], dmvely[*nx], dmvelz[*nx]; + double tilde_flatten[*nx]; + + const double onesixth = 1.0/6.0; + + for(int i=1;i < *nx-1;i++) { + drho[i] = 0.5 * (rho[i+1]-rho[i-1]); + d2rho[i] = rho[i+1] - 2.0 * rho[i] + rho[i-1]; + dpress[i] = (press[i+1]-press[i-1]); + deps[i] = 0.5 * (eps[i+1]-eps[i-1]); + dtemp[i] = 0.5 * (temp[i+1]-temp[i-1]); + dye[i] = 0.5 * (ye[i+1]-ye[i-1]); + dvelx[i] = 0.5 * (velx[i+1]-velx[i-1]); + dvely[i] = 0.5 * (vely[i+1]-vely[i-1]); + dvelz[i] = 0.5 * (velz[i+1]-velz[i-1]); + } + + for(int i=1;i<*nx-1;i++) { + steep(rho,drho,dmrho,i); + steep(eps,deps,dmeps,i); + steep(temp,dtemp,dmtemp,i); + steep(ye,dye,dmye,i); + steep(velx,dvelx,dmvelx,i); + steep(vely,dvely,dmvely,i); + steep(velz,dvelz,dmvelz,i); + } + + for(int i=1;i<*nx-1;i++) { + rhoplus[i] = 0.5 * (rho[i] + rho[i+1]) + + (dmrho[i] - dmrho[i+1]) * onesixth; + epsplus[i] = 0.5 * (eps[i] + eps[i+1]) + + (dmeps[i] - dmeps[i+1]) * onesixth; + tempplus[i] = 0.5 * (temp[i] + temp[i+1]) + + (dmtemp[i] - dmtemp[i+1]) * onesixth; + yeplus[i] = 0.5 * (ye[i] + ye[i+1]) + + (dmye[i] - dmye[i+1]) * onesixth; + velxplus[i] = 0.5 * (velx[i] + velx[i+1]) + + (dmvelx[i] - dmvelx[i+1]) * onesixth; + velyplus[i] = 0.5 * (vely[i] + vely[i+1]) + + (dmvely[i] - dmvely[i+1]) * onesixth; + velzplus[i] = 0.5 * (velz[i] + velz[i+1]) + + (dmvelz[i] - dmvelz[i+1]) * onesixth; + } + + for(int i=1;i<*nx-2;i++) { + rhoplus[i] = approx_at_cell_interface(rho,i); + epsplus[i] = approx_at_cell_interface(eps,i); + tempplus[i] = approx_at_cell_interface(temp,i); + yeplus[i] = approx_at_cell_interface(ye,i); + velxplus[i] = approx_at_cell_interface(velx,i); + velyplus[i] = approx_at_cell_interface(vely,i); + velzplus[i] = approx_at_cell_interface(velz,i); + } + + for(int i=1;i<*nx-2;i++) { + rhominus[i+1] = rhoplus[i]; + epsminus[i+1] = epsplus[i]; + tempminus[i+1] = tempplus[i]; + yeminus[i+1] = yeplus[i]; + velxminus[i+1] = velxplus[i]; + velyminus[i+1] = velyplus[i]; + velzminus[i+1] = velzplus[i]; + } + + const double ppm_epsilon_shock = 0.01; + const double ppm_eta1 = 20.0; + const double ppm_eta2 = 0.05; + const double ppm_k0 = 0.2; + const double ppm_epsilon = 0.33; + const double ppm_omega1 = 0.75; + const double ppm_omega2 = 10.0; + const double ppm_small = 1.0e-7; + + // steepening + for(int i=2;i<*nx-2;i++) { + double etatilde = 0.0; + if ( d2rho[i+1]*d2rho[i-1] < 0.0 + && ( fabs(rho[i+1]-rho[i-1]) - ppm_epsilon_shock + * MIN(fabs(rho[i+1]), + fabs(rho[i-1])) > 0.0) ) + { + etatilde = (rho[i-2] - rho[i+2] + 4.0 * drho[i]) / (drho[i] * 12.0); + } + double eta = MAX(0.0, MIN(1.0, ppm_eta1 * (etatilde - ppm_eta2))); + if (ppm_k0 * fabs(drho[i]) * MIN(press[i-1],press[i+1]) + < fabs(dpress[i]) * MIN(rho[i-1], rho[i+1])) + { + eta = 0.0; + } + rhominus[i] = rhominus[i] * (1.0 - eta) + + (rho[i-1] + 0.5 * dmrho[i-1]) * eta; + rhoplus[i] = rhoplus[i] * (1.0 - eta) + + (rho[i+1] - 0.5 * dmrho[i+1]) * eta; + } + + // flattening + for(int i=2;i<*nx-2;i++) { + double dpress2 = press[i+2] - press[i-2]; + double dvel = velx[i+1] - velx[i-1]; + double w=0.0; + if ( (fabs(dpress[i]) > ppm_epsilon * MIN(press[i-1],press[i+1])) + && (dvel < 0.0) ) + { + w = 1.0; + } + if (fabs(dpress2) < ppm_small) + { + tilde_flatten[i] = 1.0; + } + else + { + tilde_flatten[i] = MAX(0.0, 1.0 - w * MAX(0.0, + ppm_omega2 * (dpress[i] / dpress2 - ppm_omega1))); + } + } + + for(int i=2;i<*nx-2;i++) { + double flatten = tilde_flatten[i]; + rhoplus[i] = flatten * rhoplus[i] + (1.0 - flatten) * rho[i]; + rhominus[i] = flatten * rhominus[i] + (1.0 - flatten) * rho[i]; + epsplus[i] = flatten * epsplus[i] + (1.0 - flatten) * eps[i]; + epsminus[i] = flatten * epsminus[i] + (1.0 - flatten) * eps[i]; + tempplus[i] = flatten * tempplus[i] + (1.0 - flatten) * temp[i]; + tempminus[i] = flatten * tempminus[i] + (1.0 - flatten) * temp[i]; + yeplus[i] = flatten * yeplus[i] + (1.0 - flatten) * ye[i]; + yeminus[i] = flatten * yeminus[i] + (1.0 - flatten) * ye[i]; + velxplus[i] = flatten * velxplus[i] + (1.0 - flatten) * velx[i]; + velxminus[i] = flatten * velxminus[i] + (1.0 - flatten) * velx[i]; + velyplus[i] = flatten * velyplus[i] + (1.0 - flatten) * vely[i]; + velyminus[i] = flatten * velyminus[i] + (1.0 - flatten) * vely[i]; + velzplus[i] = flatten * velzplus[i] + (1.0 - flatten) * velz[i]; + velzminus[i] = flatten * velzminus[i] + (1.0 - flatten) * velz[i]; + } // flattening + + for(int i=2;i<*nx-2;i++) { + monotonize(rhominus,rho,rhoplus,i); + monotonize(epsminus,eps,epsplus,i); + monotonize(tempminus,temp,tempplus,i); + monotonize(yeminus,ye,yeplus,i); + monotonize(velxminus,velx,velxplus,i); + monotonize(velyminus,vely,velyplus,i); + monotonize(velzminus,velz,velzplus,i); + } + + return; +} + + + +CCTK_FCALL void CCTK_FNAME(PPMtyc)(int* restrict nx, + double* restrict dx, + double* restrict rho, + double* restrict velx, + double* restrict vely, + double* restrict velz, + double* restrict eps, + double* restrict temp, + double* restrict ye, + double* restrict press, + double* restrict rhominus, + double* restrict velxminus, + double* restrict velyminus, + double* restrict velzminus, + double* restrict epsminus, + double* restrict tempminus, + double* restrict yeminus, + double* restrict rhoplus, + double* restrict velxplus, + double* restrict velyplus, + double* restrict velzplus, + double* restrict epsplus, + double* restrict tempplus, + double* restrict yeplus) { + +void SimplePPM_1d_tyc(nx, dx, + rho, + velx, + vely, + velz, + eps, + temp, + ye, + press, + rhominus, + velxminus, + velyminus, + velzminus, + epsminus, + tempminus, + yeminus, + rhoplus, + velxplus, + velyplus, + velzplus, + epsplus, + tempplus, + yeplus); +} + File [added]: GRHydro_PPM.h Delta lines: +56 -0 =================================================================== --- trunk/src/GRHydro_PPM.h (rev 0) +++ trunk/src/GRHydro_PPM.h 2013-07-06 18:11:31 UTC (rev 560) @@ -0,0 +1,56 @@ +#include <assert.h> +#include <math.h> +#include <stdio.h> +#include <stdlib.h> + +#if 1 +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#endif + +#define MIN(a,b) (((a)<(b))?(a):(b)) +#define MAX(a,b) (((a)>(b))?(a):(b)) + +static inline void steep(double *x, double *dx, double* dmx, const int i) { + if ( (x[i+1] - x[i]) * (x[i]-x[i-1]) > 0.0 ) { + dmx[i] = copysign(1.0,dx[i]) * MIN(fabs(dx[i]), + MIN(2.0*fabs(x[i]-x[i-1]), + 2.0*fabs(x[i+1]-x[i]))); + } else { + dmx[i] = 0.0; + } +} + +static inline double approx_at_cell_interface(double* a, const int i) { + return 7.0/12.0*(a[i]+a[i+1]) - 1.0/12.0*(a[i-1]+a[i+2]); +} + +static inline void monotonize(double* restrict xminus, + double* restrict x, + double* restrict xplus, + const int i) { + + if ( !(xplus[i]==x[i] && x[i]==xminus[i]) + && ( (xplus[i]-x[i])*(x[i]-xminus[i]) <= 0.0 ) ) + { + xminus[i] = x[i]; + xplus[i] = x[i]; + } else if( 6.0 * (xplus[i]-xminus[i]) * + (x[i]-0.5*(xplus[i]+xminus[i])) > + (xplus[i]-xminus[i])*(xplus[i]-xminus[i]) ) + { + xminus[i] = 3.0*x[i]-2.0*xplus[i]; + } else if( 6.0 * (xplus[i]-xminus[i]) * + (x[i]-0.5*(xplus[i]+xminus[i])) < + -(xplus[i]-xminus[i])*(xplus[i]-xminus[i]) ) + { + xplus[i] = 3.0*x[i]-2.0*xminus[i]; + } + + return; +} + + + File [added]: GRHydro_PPMReconstruct_drv_opt.F90 Delta lines: +212 -0 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.F90 (rev 0) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.F90 2013-07-06 18:11:31 UTC (rev 560) @@ -0,0 +1,212 @@ + /*@@ + @file GRHydro_PPMReconstruct_drv.F90 + @date Tue Jul 19 13:22:03 EDT 2011 + @author Bruno C. Mundim, Joshua Faber, Christian D. Ott + @desc + Driver routine to perform the PPM reconstructions. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Functions.h" + +#include "SpaceMask.h" + +#define velx(i,j,k) vup(i,j,k,1) +#define vely(i,j,k) vup(i,j,k,2) +#define velz(i,j,k) vup(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) + + + /*@@ + @routine GRHydro_PPMReconstruct_drv + @date Tue Jul 19 13:24:34 EDT 2011 + @author Luca Baiotti, Ian Hawke, Bruno C. Mundim, Joshua Faber, Christian D. Ott + @desc + A driver routine to do PPM reconstructions. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) + + implicit none + + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET betaa, betab, betac + TARGET betax, betay, betaz + TARGET lvel, vel + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + integer :: nx, ny, nz, i, j, k + + logical, dimension(:,:,:), allocatable :: trivial_rp + + CCTK_INT :: type_bitsx, trivialx, not_trivialx, & + &type_bitsy, trivialy, not_trivialy, & + &type_bitsz, trivialz, not_trivialz + + CCTK_INT :: ierr + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + + logical :: apply_enhanced_ppm + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + beta1 => betaa + beta2 => betab + beta3 => betac + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + beta1 => betax + beta2 => betay + beta3 => betaz + vup => vel + end if + + allocate(trivial_rp(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) + if (ierr .ne. 0) then + call CCTK_WARN(0, "Allocation problems with trivial_rp") + end if + + call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX") + call SpaceMask_GetStateBits(trivialx, "Hydro_RiemannProblemX", & + &"trivial") + call SpaceMask_GetStateBits(not_trivialx, "Hydro_RiemannProblemX", & + &"not_trivial") + call SpaceMask_GetTypeBits(type_bitsy, "Hydro_RiemannProblemY") + call SpaceMask_GetStateBits(trivialy, "Hydro_RiemannProblemY", & + &"trivial") + call SpaceMask_GetStateBits(not_trivialy, "Hydro_RiemannProblemY", & + &"not_trivial") + call SpaceMask_GetTypeBits(type_bitsz, "Hydro_RiemannProblemZ") + call SpaceMask_GetStateBits(trivialz, "Hydro_RiemannProblemZ", & + &"trivial") + call SpaceMask_GetStateBits(not_trivialz, "Hydro_RiemannProblemZ", & + &"not_trivial") + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + ! if use_enhanced_ppm, allow old PPM on one level + if (GRHydro_oppm_reflevel .eq. (-1) .or. & + GRHydro_reflevel .ne. GRHydro_oppm_reflevel) then + apply_enhanced_ppm = use_enhanced_ppm .ne. 0 + else + apply_enhanced_ppm = .false. + end if + + + +!!$ PPM starts: + if (flux_direction == 1) then + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + + trivial_rp(:,j,k) = .false. + call PPMtyc(nx,CCTK_DELTA_SPACE(1),rho(:,j,k),velx(:,j,k),& + vely(:,j,k),velz(:,j,k),eps(:,j,k),temperature(:,j,k),Y_e(:,j,k),& + press(:,j,k),rhominus(:,j,k),& + velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),& + epsminus(:,j,k),tempminus(:,j,k), y_e_minus(:,j,k),& + rhoplus(:,j,k),& + velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& + epsplus(:,j,k),tempplus(:,j,k),y_e_plus(:,j,k)) + do i = 1, nx + if (trivial_rp(i,j,k)) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) + else + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else if (flux_direction == 2) then + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + trivial_rp(j,:,k) = .false. + call PPMtyc(ny,CCTK_DELTA_SPACE(2),rho(j,:,k),& + vely(j,:,k),velz(j,:,k),velx(j,:,k),& + eps(j,:,k),temperature(j,:,k),Y_e(j,:,k), & + press(j,:,k),rhominus(j,:,k),& + velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),& + epsminus(j,:,k),tempminus(j,:,k), y_e_minus(j,:,k),rhoplus(j,:,k),& + velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& + epsplus(j,:,k),tempplus(j,:,k),y_e_plus(j,:,k)) + do i = 1, ny + if (trivial_rp(j,i,k)) then + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) + else + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) + end if + end do + end do + end do + !$OMP END PARALLEL DO + + else if (flux_direction == 3) then + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + trivial_rp(j,k,:) = .false. + call PPMtyc(nz,CCTK_DELTA_SPACE(3),rho(j,k,:),& + velz(j,k,:),velx(j,k,:),vely(j,k,:),& + eps(j,k,:),temperature(j,k,:),Y_e(j,k,:), press(j,k,:),rhominus(j,k,:),& + velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),& + epsminus(j,k,:),tempminus(j,k,:),y_e_minus(j,k,:),rhoplus(j,k,:),& + velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& + epsplus(j,k,:),tempplus(j,k,:),y_e_plus(j,k,:)) + do i = 1, nz + if (trivial_rp(j,k,i)) then + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) + else + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if +!!$ PPM ends. + + deallocate(trivial_rp) + +end subroutine GRHydro_PPMReconstruct_drv +

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 559)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:11 PM Modified: /trunk/ param.ccl Log: GRHydro: make atmosphere threshold steerable on recovery From: Christian David Ott <cott(a)zwicky.(none)> File Changes: Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +2 -2 =================================================================== --- trunk/param.ccl 2013-07-06 18:11:25 UTC (rev 558) +++ trunk/param.ccl 2013-07-06 18:11:28 UTC (rev 559) @@ -482,12 +482,12 @@ 0: :: "" } 1.e-10 -REAL rho_abs_min "A minimum rho below which evolution is turned off (atmosphere). If negative, the relative minimum will be used instead." +REAL rho_abs_min "A minimum rho below which evolution is turned off (atmosphere). If negative, the relative minimum will be used instead." STEERABLE=recover { -1.0: :: "" } -1.0 -REAL rho_rel_min "A minimum relative rho (rhomin = centden * rho_rel_min) below which evolution is turned off (atmosphere). Only used if rho_abs_min < 0.0" +REAL rho_rel_min "A minimum relative rho (rhomin = centden * rho_rel_min) below which evolution is turned off (atmosphere). Only used if rho_abs_min < 0.0" STEERABLE=recover { 0: :: "" } 1.e-9

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 558)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:11 PM Modified: /trunk/ schedule.ccl Log: GRHydro: activate storage for Y_e, temperature, entropy when we evolve it From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/ ================== File [modified]: schedule.ccl Delta lines: +3 -0 =================================================================== --- trunk/schedule.ccl 2013-07-06 18:11:22 UTC (rev 557) +++ trunk/schedule.ccl 2013-07-06 18:11:25 UTC (rev 558) @@ -60,10 +60,13 @@ { STORAGE: Y_e_con[timelevels] STORAGE: Y_e_con_rhs, Y_e_con_flux, Y_e_plus, Y_e_minus + STORAGE: HydroBase::Y_e[timelevels] } if(CCTK_Equals(temperature_evolution_method,"GRHydro")) { STORAGE: tempplus, tempminus + STORAGE: HydroBase::temperature[timelevels] + STORAGE: HydroBase::entropy[timelevels] } if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) {

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 557)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:11 PM Added: /trunk/src/ GRHydro_Wrappers.F90 Modified: /trunk/ interface.ccl /trunk/src/ make.code.defn Log: GRHydro: provide wrapper routines that take detg as arguments then call internal routines taking sqrt(detg) as argument clean up and correct function prototypes in interface.ccl From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [added]: GRHydro_Wrappers.F90 Delta lines: +366 -0 =================================================================== --- trunk/src/GRHydro_Wrappers.F90 (rev 0) +++ trunk/src/GRHydro_Wrappers.F90 2013-07-06 18:11:22 UTC (rev 557) @@ -0,0 +1,366 @@ + /*@@ + @file GRHydro_Wrappers.F90 + @date Sat Jun 29 18:57:07 PDT 2013 + @author Roland Haas + @desc + Wrapper routines to provide a consistent interface to the aliased + functions even if the interface to the underlying implementation changes. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Functions.h" + +subroutine Con2PrimPolyWrapper(handle, dens, sx, sy, sz, tau, rho, & + velx, vely, velz, epsilon, press, w_lorentz, uxx, uxy, uxz, uyy, & + uyz, uzz, detg, x, y, z, r, GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed) + + use Con2Prim_fortran_interfaces + + implicit none + + DECLARE_CCTK_FUNCTIONS + + CCTK_REAL, INTENT(INOUT) :: dens, sx, sy, sz + CCTK_REAL, INTENT(OUT) :: tau + CCTK_REAL, INTENT(INOUT) :: rho + CCTK_REAL, INTENT(INOUT) :: velx, vely, velz, epsilon, press, w_lorentz + CCTK_REAL, INTENT(IN) :: uxx, uxy, uxz, uyy, uyz, uzz, detg + CCTK_REAL, INTENT(IN) :: x, y, z, r, GRHydro_rho_min + + CCTK_INT, INTENT(IN) :: handle, GRHydro_reflevel + CCTK_REAL, INTENT(OUT) :: GRHydro_C2P_failed + + CCTK_REAL :: sdetg + + sdetg = sqrt(detg) + + call Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & + velx, vely, velz, epsilon, press, w_lorentz, uxx, uxy, uxz, uyy, & + uyz, uzz, sdetg, x, y, z, r, GRHydro_rho_min, GRHydro_reflevel, & + GRHydro_C2P_failed) + +end subroutine + +subroutine Con2PrimGenWrapper(handle, dens, sx, sy, sz, tau, rho, velx, vely, & + velz, epsilon, pressure, w_lorentz, uxx, uxy, & + uxz, uyy, uyz, uzz, det, x, y, z, r, epsnegative, & + GRHydro_rho_min, pmin, epsmin, GRHydro_reflevel, & + retval) + + use Con2Prim_fortran_interfaces + + implicit none + + DECLARE_CCTK_FUNCTIONS + + CCTK_INT, INTENT(IN) :: handle + CCTK_REAL, INTENT(INOUT) :: dens, sx, sy, sz, tau + CCTK_REAL, INTENT(INOUT) :: rho, velx, vely, velz, epsilon, pressure + CCTK_REAL, INTENT(INOUT) :: w_lorentz + CCTK_REAL, INTENT(IN) :: uxx, uxy, uxz, uyy, uyz, uzz, det + CCTK_REAL, INTENT(IN) :: x, y, z, r + CCTK_INT, INTENT(OUT) :: epsnegative + CCTK_REAL, INTENT(IN) :: GRHydro_rho_min, pmin, epsmin + CCTK_INT, INTENT(IN) :: GRHydro_reflevel + CCTK_REAL, INTENT(OUT) :: retval + + CCTK_REAL :: sdetg + CCTK_INT, parameter :: cctk_iteration = -1, ii = -1, jj = -1, kk = -1 + logical :: log_epsnegative + + sdetg = sqrt(det) + + call Con2Prim_pt(cctK_iteration, ii, jj, kk, & + handle, dens, sx, sy, sz, tau, rho, velx, vely, & + velz, epsilon, pressure, w_lorentz, uxx, uxy, & + uxz, uyy, uyz, uzz, sdetg, x, y, z, r, log_epsnegative, & + GRHydro_rho_min, pmin, epsmin, GRHydro_reflevel, & + retval) + if(log_epsnegative) then + epsnegative = 1 + else + epsnegative = 0 + end if + +end subroutine + +subroutine Con2PrimGenMWrapper(handle, keytemp, prec, gamma_eos, dens, sx, sy, sz, & + tau, Bconsx, Bconsy, Bconsz, y_e, temp, rho, velx, & + vely, velz, epsilon, pressure, Bvecx, Bvecy, Bvecz, & + Bvecsq, w_lorentz, gxx, gxy, gxz, gyy, gyz, gzz, uxx, & + uxy, uxz, uyy, uyz, uzz, det, epsnegative, retval) + + use Con2PrimM_fortran_interfaces + + implicit none + + DECLARE_CCTK_FUNCTIONS + + CCTK_INT, INTENT(IN) :: handle, keytemp + CCTK_REAL, INTENT(IN) :: prec, gamma_eos + CCTK_REAL, INTENT(INOUT) :: dens, sx, sy, sz, tau + CCTK_REAL, INTENT(IN) :: Bconsx, Bconsy, Bconsz + CCTK_REAL, INTENT(INOUT) :: y_e, temp + CCTK_REAL, INTENT(INOUT) :: rho, velx, vely, velz, epsilon, pressure + CCTK_REAL, INTENT(OUT) :: Bvecx, Bvecy, Bvecz, Bvecsq + CCTK_REAL, INTENT(INOUT) :: w_lorentz + CCTK_REAL, INTENT(IN) :: gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, INTENT(IN) :: uxx, uxy, uxz, uyy, uyz, uzz + CCTK_REAL, INTENT(IN) :: det + CCTK_INT, INTENT(OUT) :: epsnegative + CCTK_REAL, INTENT(OUT) :: retval + + CCTK_REAL :: sdetg + + sdetg = sqrt(det) + + if(handle.eq.1 .or. handle.eq.2) then + call GRHydro_Con2PrimM_ptold(handle, gamma_eos, dens, sx, sy, sz, & + tau, Bconsx, Bconsy, Bconsz, rho, velx, & + vely, velz, epsilon, pressure, Bvecx, Bvecy, Bvecz, & + Bvecsq, w_lorentz, gxx, gxy, gxz, gyy, gyz, gzz, uxx, & + uxy, uxz, uyy, uyz, uzz, sdetg, epsnegative, retval) + else + call GRHydro_Con2PrimM_pt2(handle, keytemp, prec, prec, gamma_eos, dens, sx, sy, sz, & + tau, Bconsx, Bconsy, Bconsz, y_e, temp, rho, velx, & + vely, velz, epsilon, pressure, Bvecx, Bvecy, Bvecz, & + Bvecsq, w_lorentz, gxx, gxy, gxz, gyy, gyz, gzz, uxx, & + uxy, uxz, uyy, uyz, uzz, sdetg, epsnegative, retval) + end if + +end subroutine + +subroutine Con2PrimPolyMWrapper(handle, gamma_eos, dens, sx, sy, sz, sc,& + Bconsx, Bconsy, Bconsz, rho, velx, vely, velz,& + epsilon, pressure, Bvecx, Bvecy, Bvecz, Bvecsq,& + w_lorentz, gxx, gxy, gxz, gyy, gyz, gzz, uxx,& + uxy, uxz, uyy, uyz, uzz, det, epsnegative,& + retval) + + use Con2PrimM_fortran_interfaces + + implicit none + + DECLARE_CCTK_FUNCTIONS + + CCTK_INT, INTENT(IN) :: handle + CCTK_REAL, INTENT(IN) :: gamma_eos + CCTK_REAL, INTENT(INOUT) :: dens, sx, sy, sz, sc + CCTK_REAL, INTENT(INOUT) :: Bconsx, Bconsy, Bconsz + CCTK_REAL, INTENT(INOUT) :: rho, velx, vely, velz, epsilon, pressure + CCTK_REAL, INTENT(OUT) :: Bvecx, Bvecy, Bvecz + CCTK_REAL, INTENT(OUT) :: Bvecsq + CCTK_REAL, INTENT(INOUT) :: w_lorentz + CCTK_REAL, INTENT(IN) :: gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, INTENT(IN) :: uxx, uxy, uxz, uyy, uyz, uzz + CCTK_REAL, INTENT(IN) :: det + CCTK_INT, INTENT(OUT) :: epsnegative + CCTK_REAL, INTENT(OUT) :: retval + + CCTK_REAL :: sdetg + + sdetg = sqrt(det) + + call GRHydro_Con2PrimM_Polytype_pt(handle, gamma_eos, dens, sx, sy, sz, sc, & + Bconsx, Bconsy, Bconsz, rho, velx, vely, velz, & + epsilon, pressure, Bvecx, Bvecy, Bvecz, Bvecsq, & + w_lorentz, gxx, gxy, gxz, gyy, gyz, gzz, uxx, & + uxy, uxz, uyy, uyz, uzz, sdetg, epsnegative, & + retval) +end subroutine + +subroutine Prim2ConGenWrapper(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, dens,& + sx, sy, sz, tau, rho, velx, vely, velz, epsilon,& + press, w_lorentz) + + use Con2Prim_fortran_interfaces + + implicit none + + DECLARE_CCTK_FUNCTIONS + + CCTK_INT, INTENT(IN) :: handle + CCTK_REAL, INTENT(IN) :: gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, INTENT(IN) :: det + CCTK_REAL, INTENT(OUT) :: dens, sx, sy, sz, tau + CCTK_REAL, INTENT(IN) :: rho, velx, vely, velz, epsilon + CCTK_REAL, INTENT(OUT) :: press, w_lorentz + + CCTK_REAL :: sdetg + + sdetg = sqrt(det) + + call prim2con(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, dens,& + sx, sy, sz, tau, rho, velx, vely, velz, epsilon,& + press, w_lorentz) +end subroutine + +subroutine Prim2ConPolyWrapper(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, dens,& + sx, sy, sz, tau, rho, velx, vely, velz, epsilon,& + press, w_lorentz) + + use Con2Prim_fortran_interfaces + + implicit none + + DECLARE_CCTK_FUNCTIONS + + CCTK_INT, INTENT(IN) :: handle + CCTK_REAL, INTENT(IN) :: gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, INTENT(IN) :: det + CCTK_REAL, INTENT(OUT) :: dens, sx, sy, sz, tau + CCTK_REAL, INTENT(IN) :: rho, velx, vely, velz + CCTK_REAL, INTENT(OUT) :: epsilon, press, w_lorentz + + CCTK_REAL :: sdetg + + sdetg = sqrt(det) + + call prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, dens,& + sx, sy, sz, tau, rho, velx, vely, velz, epsilon,& + press, w_lorentz) +end subroutine + +subroutine Prim2ConGenMWrapper(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, dens,& + sx, sy, sz, tau, Bconsx, Bconsy, Bconsz, rho,& + velx, vely, velz, epsilon, press, Bvecx, Bvecy,& + Bvecz, w_lorentz) + + use Con2PrimM_fortran_interfaces + + implicit none + + DECLARE_CCTK_FUNCTIONS + + CCTK_INT, INTENT(IN) :: handle + CCTK_REAL, INTENT(IN) :: gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, INTENT(IN) :: det + CCTK_REAL, INTENT(OUT) :: dens, sx, sy, sz, tau + CCTK_REAL, INTENT(OUT) :: Bconsx, Bconsy, Bconsz + CCTK_REAL, INTENT(IN) :: rho, velx, vely, velz, epsilon + CCTK_REAL, INTENT(OUT) :: press + CCTK_REAL, INTENT(IN) :: Bvecx, Bvecy, Bvecz + CCTK_REAL, INTENT(OUT) :: w_lorentz + + CCTK_REAL :: sdetg + + sdetg = sqrt(det) + + call prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, sdetg, dens,& + sx, sy, sz, tau, Bconsx, Bconsy, Bconsz, rho,& + velx, vely, velz, epsilon, press, Bvecx, Bvecy,& + Bvecz, w_lorentz) +end subroutine + +subroutine Prim2ConGenM_hotWrapper(handle, GRHydro_reflevel, i, j, k, x, y, z,& + gxx, gxy, gxz, gyy, gyz, gzz, det, dens, sx,& + sy, sz, tau, Bconsx, Bconsy, Bconsz, rho,& + velx, vely, velz, epsilon, press, Bvecx,& + Bvecy, Bvecz, w_lorentz, temperature, Y_e) + + use Con2PrimM_fortran_interfaces + + implicit none + + DECLARE_CCTK_FUNCTIONS + + CCTK_INT, INTENT(IN) :: handle, GRHydro_reflevel + CCTK_INT, INTENT(IN) :: i, j, k + CCTK_REAL, INTENT(IN) :: x, y, z + CCTK_REAL, INTENT(IN) :: gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, INTENT(IN) :: det + CCTK_REAL, INTENT(OUT) :: dens, sx, sy, sz, tau + CCTK_REAL, INTENT(OUT) :: Bconsx, Bconsy, Bconsz + CCTK_REAL, INTENT(IN) :: rho, velx, vely, velz, epsilon + CCTK_REAL, INTENT(OUT) :: press + CCTK_REAL, INTENT(IN) :: Bvecx, Bvecy, Bvecz + CCTK_REAL, INTENT(OUT) :: w_lorentz + CCTK_REAL, INTENT(INOUT) :: temperature + CCTK_REAL, INTENT(IN) :: Y_e + + CCTK_REAL :: sdetg + + sdetg = sqrt(det) + + call prim2conM_hot(handle, GRHydro_reflevel, i, j, k, x, y, z, gxx,& + gxy, gxz, gyy, gyz, gzz, sdetg, dens, sx, sy, sz,& + tau, Bconsx, Bconsy, Bconsz, rho, velx, vely, velz,& + epsilon, press, Bvecx, Bvecy, Bvecz, w_lorentz,& + temperature, Y_e) +end subroutine + +subroutine Prim2ConPolyMWrapper(handle, gxx, gxy, gxz, gyy, gyz, gzz, det,& + dens, sx, sy, sz, tau, Bconsx, Bconsy, Bconsz,& + rho, velx, vely, velz, epsilon, press, Bvecx,& + Bvecy, Bvecz, w_lorentz) + + use Con2PrimM_fortran_interfaces + + implicit none + + DECLARE_CCTK_FUNCTIONS + + CCTK_INT, INTENT(IN) :: handle + CCTK_REAL, INTENT(IN) :: gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, INTENT(IN) :: det + CCTK_REAL, INTENT(OUT) :: dens, sx, sy, sz, tau + CCTK_REAL, INTENT(OUT) :: Bconsx, Bconsy, Bconsz + CCTK_REAL, INTENT(IN) :: rho, velx, vely, velz + CCTK_REAL, INTENT(OUT) :: epsilon + CCTK_REAL, INTENT(OUT) :: press + CCTK_REAL, INTENT(IN) :: Bvecx, Bvecy, Bvecz + CCTK_REAL, INTENT(OUT) :: w_lorentz + + CCTK_REAL :: sdetg + + sdetg = sqrt(det) + + call prim2conpolytypeM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det,& + dens, sx, sy, sz, tau, Bconsx, Bconsy, Bconsz,& + rho, velx, vely, velz, epsilon, press, Bvecx,& + Bvecy, Bvecz, w_lorentz) +end subroutine + +subroutine Con2PrimGenMeeWrapper(handle, keytemp, prec, gamma_eos, dens, sx, sy,& + sz, tau, Bconsx, Bconsy, Bconsz, entropycons,& + y_e, temp, rho, velx, vely, velz, epsilon,& + pressure, Bvecx, Bvecy, Bvecz, Bvecsq,& + w_lorentz, gxx, gxy, gxz, gyy, gyz, gzz, uxx,& + uxy, uxz, uyy, uyz, uzz, det, epsnegative,& + retval) + + use Con2Prim_fortran_interfaces + + implicit none + + DECLARE_CCTK_FUNCTIONS + + CCTK_INT, INTENT(IN) :: handle, keytemp + CCTK_REAL, INTENT(IN) :: prec, gamma_eos + CCTK_REAL, INTENT(INOUT) :: dens, sx, sy, sz, tau + CCTK_REAL, INTENT(IN) :: Bconsx, Bconsy, Bconsz + CCTK_REAL, INTENT(INOUT) :: entropycons, y_e, temp + CCTK_REAL, INTENT(INOUT) :: rho, velx, vely, velz, epsilon, pressure + CCTK_REAL, INTENT(OUT) :: Bvecx, Bvecy, Bvecz + CCTK_REAL, INTENT(OUT) :: Bvecsq + CCTK_REAL, INTENT(INOUT) :: w_lorentz + CCTK_REAL, INTENT(IN) :: gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, INTENT(IN) :: uxx, uxy, uxz, uyy, uyz, uzz + CCTK_REAL, INTENT(IN) :: det + CCTK_INT, INTENT(OUT) :: epsnegative + CCTK_REAL, INTENT(OUT) :: retval + + CCTK_REAL :: sdetg + + sdetg = sqrt(det) + + call GRHydro_Con2PrimM_ptee(handle, keytemp, prec, gamma_eos, dens, sx, sy, sz,& + tau, Bconsx, Bconsy, Bconsz, entropycons, y_e, temp,& + rho, velx, vely, velz, epsilon, pressure, Bvecx,& + Bvecy, Bvecz, Bvecsq, w_lorentz, gxx, gxy, gxz, gyy,& + gyz, gzz, uxx, uxy, uxz, uyy, uyz, uzz, sdetg,& + epsnegative, retval) +end subroutine File [modified]: make.code.defn Delta lines: +3 -4 =================================================================== --- trunk/src/make.code.defn 2013-07-06 18:11:18 UTC (rev 556) +++ trunk/src/make.code.defn 2013-07-06 18:11:22 UTC (rev 557) @@ -1,5 +1,4 @@ -# Main make.code.defn file for thorn Hydra -# $Header$ +# Main make.code.defn file for thorn GRHydro # Source files in this directory @@ -82,9 +81,9 @@ GRHydro_MP5Reconstruct.F90 \ Con2Prim_fortran_interfaces.F90 \ Con2PrimM_fortran_interfaces.F90 \ - GRHydro_Prim2Con.c + GRHydro_Prim2Con.c \ + GRHydro_Wrappers.F90 - # Subdirectories containing source files SUBDIRS = Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +57 -97 =================================================================== --- trunk/interface.ccl 2013-07-06 18:11:18 UTC (rev 556) +++ trunk/interface.ccl 2013-07-06 18:11:22 UTC (rev 557) @@ -34,130 +34,92 @@ CCTK_REAL IN gxz, CCTK_REAL IN gyy, CCTK_REAL IN gyz, \ CCTK_REAL IN gzz) -#void FUNCTION Con2Prim(CCTK_INT IN handle, \ -# CCTK_REAL OUT dens, \ -# CCTK_REAL OUT sx, CCTK_REAL OUT sy, CCTK_REAL OUT sz, \ -# CCTK_REAL OUT tau, \ -# CCTK_REAL OUT rho, \ -# CCTK_REAL OUT velx, CCTK_REAL OUT vely, CCTK_REAL OUT velz, \ -# CCTK_REAL OUT epsilon, CCTK_REAL OUT press, \ -# CCTK_REAL OUT w_lorentz, \ -# CCTK_REAL IN uxx, CCTK_REAL IN uxy, CCTK_REAL IN uxz, \ -# CCTK_REAL IN uyy, CCTK_REAL IN uyz, CCTK_REAL IN uzz, \ -# CCTK_REAL IN det, \ -# CCTK_REAL IN x, CCTK_REAL IN y, CCTK_REAL IN z, CCTK_REAL IN r, \ -# CCTK_INT OUT epsnegative, \ -# CCTK_REAL IN GRHydro_rho_min, CCTK_REAL IN pmin, \ -# CCTK_INT IN GRHydro_reflevel, CCTK_REAL IN GRHydro_C2P_failed) - -void FUNCTION Con2PrimGen(CCTK_INT INOUT handle, CCTK_REAL INOUT dens, \ +void FUNCTION Con2PrimGen(CCTK_INT IN handle, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ CCTK_REAL INOUT tau, CCTK_REAL INOUT rho, \ CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ CCTK_REAL INOUT w_lorentz, \ - CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ - CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ - CCTK_REAL INOUT det, \ - CCTK_REAL INOUT x, CCTK_REAL INOUT y, CCTK_REAL INOUT z, \ - CCTK_REAL INOUT r, \ + CCTK_REAL IN uxx, CCTK_REAL IN uxy, CCTK_REAL IN uxz, \ + CCTK_REAL IN uyy, CCTK_REAL IN uyz, CCTK_REAL IN uzz, \ + CCTK_REAL IN det, \ + CCTK_REAL IN x, CCTK_REAL IN y, CCTK_REAL IN z, \ + CCTK_REAL IN r, \ CCTK_INT OUT epsnegative, \ - CCTK_REAL INOUT GRHydro_rho_min, CCTK_REAL INOUT pmin, \ - CCTK_REAL INOUT epsmin, \ - CCTK_INT INOUT GRHydro_reflevel, \ + CCTK_REAL IN GRHydro_rho_min, CCTK_REAL IN pmin, \ + CCTK_REAL IN epsmin, \ + CCTK_INT IN GRHydro_reflevel, \ CCTK_REAL OUT retval) -void FUNCTION Con2PrimGenM2(CCTK_INT IN handle, CCTK_INT IN rl, CCTK_INT IN ii, CCTK_INT IN jj, \ - CCTK_INT IN kk, CCTK_REAL IN x, CCTK_REAL IN y, CCTK_REAL IN z, \ - CCTK_INT IN keytemp, CCTK_REAL IN eos_prec, CCTK_REAL IN prec, \ - CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ - CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ - CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ - CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ - CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ - CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ - CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ - CCTK_REAL INOUT Bvecsq, \ - CCTK_REAL INOUT w_lorentz, \ - CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ - CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ - CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ - CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ - CCTK_REAL INOUT det, \ - CCTK_INT OUT epsnegative, \ - CCTK_REAL OUT retval) - -void FUNCTION Con2PrimGenMee(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, \ - CCTK_REAL INOUT prec, CCTK_REAL INOUT gamma_eos, \ +void FUNCTION Con2PrimGenMee(CCTK_INT IN handle, CCTK_INT IN keytemp, \ + CCTK_REAL IN prec, CCTK_REAL IN gamma_eos, \ CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ CCTK_REAL INOUT tau, \ - CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ + CCTK_REAL IN Bconsx, CCTK_REAL IN Bconsy, CCTK_REAL IN Bconsz, \ CCTK_REAL INOUT entropycons, \ CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ - CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ - CCTK_REAL INOUT Bvecsq, \ + CCTK_REAL OUT Bvecx, CCTK_REAL OUT Bvecy, CCTK_REAL OUT Bvecz, \ + CCTK_REAL OUT Bvecsq, \ CCTK_REAL INOUT w_lorentz, \ - CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ - CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ - CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ - CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ - CCTK_REAL INOUT det, \ + CCTK_REAL IN gxx, CCTK_REAL IN gxy, CCTK_REAL IN gxz, \ + CCTK_REAL IN gyy, CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ + CCTK_REAL IN uxx, CCTK_REAL IN uxy, CCTK_REAL IN uxz, \ + CCTK_REAL IN uyy, CCTK_REAL IN uyz, CCTK_REAL IN uzz, \ + CCTK_REAL IN det, \ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) -void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec,CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ +void FUNCTION Con2PrimGenM(CCTK_INT IN handle, CCTK_INT IN keytemp, CCTK_REAL IN prec,CCTK_REAL IN gamma_eos, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ - CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ + CCTK_REAL INOUT tau, CCTK_REAL IN Bconsx, CCTK_REAL IN Bconsy, CCTK_REAL IN Bconsz, \ CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ - CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ - CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ - CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ - CCTK_REAL INOUT Bvecsq, \ + CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ + CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ + CCTK_REAL OUT Bvecx, CCTK_REAL OUT Bvecy, CCTK_REAL OUT Bvecz, \ + CCTK_REAL OUT bvecsq, \ CCTK_REAL INOUT w_lorentz, \ - CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ - CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ - CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ - CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ - CCTK_REAL INOUT det, \ + CCTK_REAL IN gxx, CCTK_REAL IN gxy, CCTK_REAL IN gxz, \ + CCTK_REAL IN gyy, CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ + CCTK_REAL IN uxx, CCTK_REAL IN uxy, CCTK_REAL IN uxz, \ + CCTK_REAL IN uyy, CCTK_REAL IN uyz, CCTK_REAL IN uzz, \ + CCTK_REAL IN det, \ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) -void FUNCTION Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL OUT dens, \ - CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ - CCTK_REAL OUT sz, CCTK_REAL OUT tau, \ - CCTK_REAL OUT rho, CCTK_REAL OUT velx, \ - CCTK_REAL OUT vely, CCTK_REAL OUT velz, \ - CCTK_REAL OUT epsilon, CCTK_REAL OUT press, \ - CCTK_REAL OUT w_lorentz, CCTK_REAL IN uxx, \ +void FUNCTION Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL INOUT dens, \ + CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, \ + CCTK_REAL INOUT sz, CCTK_REAL OUT tau, \ + CCTK_REAL INOUT rho, CCTK_REAL INOUT velx, \ + CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ + CCTK_REAL INOUT epsilon, CCTK_REAL INOUT press, \ + CCTK_REAL INOUT w_lorentz, CCTK_REAL IN uxx, \ CCTK_REAL IN uxy, CCTK_REAL IN uxz, \ CCTK_REAL IN uyy, CCTK_REAL IN uyz, \ CCTK_REAL IN uzz, CCTK_REAL IN det, \ CCTK_REAL IN x, CCTK_REAL IN y, CCTK_REAL IN z, \ CCTK_REAL IN r, CCTK_REAL IN rho_min, \ - CCTK_INT IN GRHydro_reflevel, CCTK_REAL IN GRHydro_C2P_failed) + CCTK_INT IN GRHydro_reflevel, CCTK_REAL OUT GRHydro_C2P_failed) -void FUNCTION Con2PrimPolyM(CCTK_INT INOUT handle, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ +void FUNCTION Con2PrimPolyM(CCTK_INT IN handle, CCTK_REAL IN gamma_eos, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ - CCTK_REAL INOUT sc, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ + CCTK_REAL INOUT sc, CCTK_REAL IN Bconsx, CCTK_REAL IN Bconsy, CCTK_REAL IN Bconsz, \ CCTK_REAL INOUT rho, \ CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ - CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ - CCTK_REAL INOUT Bvecsq, \ + CCTK_REAL OUT Bvecx, CCTK_REAL OUT Bvecy, CCTK_REAL OUT Bvecz, \ + CCTK_REAL OUT Bvecsq, \ CCTK_REAL INOUT w_lorentz, \ - CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ - CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ - CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ - CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ - CCTK_REAL INOUT det, \ + CCTK_REAL IN gxx, CCTK_REAL IN gxy, CCTK_REAL IN gxz, \ + CCTK_REAL IN gyy, CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ + CCTK_REAL IN uxx, CCTK_REAL IN uxy, CCTK_REAL IN uxz, \ + CCTK_REAL IN uyy, CCTK_REAL IN uyz, CCTK_REAL IN uzz, \ + CCTK_REAL IN det, \ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) -#void FUNCTION Prim2Con CCTK_INT handle, CCTK_REAL gxx, CCTK_REAL gxy, CCTK_REAL gxz, CCTK_REAL gyy, CCTK_REAL gyz, CCTK_REAL gzz, CCTK_REAL det, CCTK_REAL dens, CCTK_REAL sx, CCTK_REAL sy, CCTK_REAL sz, CCTK_REAL tau, CCTK_REAL rho, CCTK_REAL velx, CCTK_REAL vely, CCTK_REAL velz, CCTK_REAL epsilon, CCTK_REAL press, CCTK_REAL w_lorentz - void FUNCTION Prim2ConGen(CCTK_INT IN handle, \ CCTK_REAL IN gxx, CCTK_REAL IN gxy, \ CCTK_REAL IN gxz, CCTK_REAL IN gyy, \ @@ -228,18 +190,16 @@ PROVIDES FUNCTION SpatialDet WITH SpatialDeterminant LANGUAGE Fortran PROVIDES FUNCTION UpperMet WITH UpperMetric LANGUAGE Fortran -#PROVIDES FUNCTION Con2Prim WITH Con2Prim_pt LANGUAGE Fortran -PROVIDES FUNCTION Con2PrimPoly WITH Con2Prim_ptPolytype LANGUAGE Fortran -PROVIDES FUNCTION Con2PrimGenM WITH GRHydro_Con2PrimM_ptold LANGUAGE Fortran -PROVIDES FUNCTION Con2PrimGenM2 WITH GRHydro_Con2PrimM_pt LANGUAGE Fortran -PROVIDES FUNCTION Con2PrimGenMee WITH GRHydro_Con2PrimM_ptee LANGUAGE Fortran -PROVIDES FUNCTION Con2PrimGen WITH Con2Prim_pt LANGUAGE Fortran -PROVIDES FUNCTION Con2PrimPolyM WITH GRHydro_Con2PrimM_Polytype_pt LANGUAGE Fortran -PROVIDES FUNCTION Prim2ConGen WITH prim2con LANGUAGE Fortran -PROVIDES FUNCTION Prim2ConPoly WITH prim2conpolytype LANGUAGE Fortran -PROVIDES FUNCTION Prim2ConGenM WITH prim2conM LANGUAGE Fortran -PROVIDES FUNCTION Prim2ConGenM_hot WITH prim2conM_hot LANGUAGE Fortran -PROVIDES FUNCTION Prim2ConPolyM WITH prim2conpolytypeM LANGUAGE Fortran +PROVIDES FUNCTION Con2PrimPoly WITH Con2PrimPolyWrapper LANGUAGE Fortran +PROVIDES FUNCTION Con2PrimGenM WITH Con2PrimGenMWrapper LANGUAGE Fortran +PROVIDES FUNCTION Con2PrimGenMee WITH Con2PrimGenMeeWrapper LANGUAGE Fortran +PROVIDES FUNCTION Con2PrimGen WITH Con2PrimGenWrapper LANGUAGE Fortran +PROVIDES FUNCTION Con2PrimPolyM WITH Con2PrimPolyMWrapper LANGUAGE Fortran +PROVIDES FUNCTION Prim2ConGen WITH Prim2ConGenWrapper LANGUAGE Fortran +PROVIDES FUNCTION Prim2ConPoly WITH Prim2ConPolyWrapper LANGUAGE Fortran +PROVIDES FUNCTION Prim2ConGenM WITH Prim2ConGenMWrapper LANGUAGE Fortran +PROVIDES FUNCTION Prim2ConGenM_hot WITH Prim2ConGenM_hotWrapper LANGUAGE Fortran +PROVIDES FUNCTION Prim2ConPolyM WITH Prim2ConPolyMWrapper LANGUAGE Fortran #################################################### ### Functions provided by MoL for registration ###

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 556)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:11 PM Added: /trunk/src/ GRHydro_Prim2Con.c Modified: /trunk/src/ GRHydro_Con2Prim.F90, GRHydro_ENOReconstruct_drv.F90, GRHydro_MP5Reconstruct_drv.F90, GRHydro_PPMMReconstruct_drv.F90, GRHydro_PPMReconstruct_drv.F90, GRHydro_Reconstruct.F90, GRHydro_ReconstructPoly.F90, GRHydro_TVDReconstruct_drv.F90, GRHydro_WENOReconstruct_drv.F90, make.code.defn Log: GRHydro: re-implement prim2con routine in C * new experimental prim2con routine that is about twice as fast as the old one. * Right now tested only for simple EOS. * Right now handles only prim2con call after reconstruction and must be enabled by changing comments in GRHydro_Reconstruct.F90 * remove duplicate code in all reconstruction *_drv.F90 routines and just have the initialization of plus and minus variables in the main GRHydro_Reconstruct.F90 routine From: Christian Ott <cott(a)bethe.tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: <Unavailable> File [modified]: GRHydro_ENOReconstruct_drv.F90 Delta lines: <Unavailable> File [modified]: GRHydro_MP5Reconstruct_drv.F90 Delta lines: <Unavailable> File [modified]: GRHydro_PPMMReconstruct_drv.F90 Delta lines: <Unavailable> File [modified]: GRHydro_PPMReconstruct_drv.F90 Delta lines: <Unavailable> File [added]: GRHydro_Prim2Con.c Delta lines: <Unavailable> File [modified]: GRHydro_Reconstruct.F90 Delta lines: <Unavailable> File [modified]: GRHydro_ReconstructPoly.F90 Delta lines: <Unavailable> File [modified]: GRHydro_TVDReconstruct_drv.F90 Delta lines: <Unavailable> File [modified]: GRHydro_WENOReconstruct_drv.F90 Delta lines: <Unavailable> File [modified]: make.code.defn Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 555)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:11 PM Modified: /trunk/ interface.ccl, schedule.ccl /trunk/src/ GRHydro_BvecfromAvec.F90, GRHydro_Con2Prim.F90, GRHydro_Con2PrimAM.F90, GRHydro_Con2PrimHot.F90, GRHydro_Con2PrimM.F90, GRHydro_Con2PrimM_polytype_pt.c, GRHydro_Con2PrimM_pt.c, GRHydro_Con2PrimM_pt_EOSOmni.c, GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c, GRHydro_Con2PrimM_pt_EOSOmniold.c, GRHydro_Con2PrimM_ptee.c, GRHydro_EoSChangeGamma.F90, GRHydro_FluxSplit.F90, GRHydro_Macros.h, GRHydro_Prim2Con.F90, GRHydro_Prim2ConAM.F90, GRHydro_Prim2ConM.F90, GRHydro_Source.F90, GRHydro_SourceAM.F90, GRHydro_SourceM.F90, GRHydro_UpdateMask.F90, GRHydro_UpdateMaskM.F90, Utils.F90 Log: GRHydro: add grid function for sqrt(detg) * add new 1 tl grid function sdetg that stores the sqrt of the determinent of the 3-metric. * replace lots of re-computation of det by use of this grid function From: Christian Ott <cott(a)bethe.tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BvecfromAvec.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2Prim.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimAM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimHot.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM_polytype_pt.c Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM_pt.c Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM_pt_EOSOmni.c Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM_pt_EOSOmniold.c Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM_ptee.c Delta lines: <Unavailable> File [modified]: GRHydro_EoSChangeGamma.F90 Delta lines: <Unavailable> File [modified]: GRHydro_FluxSplit.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Macros.h Delta lines: <Unavailable> File [modified]: GRHydro_Prim2Con.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Prim2ConAM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Prim2ConM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Source.F90 Delta lines: <Unavailable> File [modified]: GRHydro_SourceAM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_SourceM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_UpdateMask.F90 Delta lines: <Unavailable> File [modified]: GRHydro_UpdateMaskM.F90 Delta lines: <Unavailable> File [modified]: Utils.F90 Delta lines: <Unavailable> Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: <Unavailable> File [modified]: schedule.ccl Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 554)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:11 PM Modified: /trunk/src/ GRHydro_Prim2Con.F90, GRHydro_Prim2ConM.F90 Log: GRHydro: pass correct temperature o prim2con call on coarse levels Fix a bug in prim2con_hot that was causing temp to go negative near the shock in 3D CCSN simulations with WENO reconstruction. From: abdik <abdik(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.F90 Delta lines: +1 -0 =================================================================== --- trunk/src/GRHydro_Prim2Con.F90 2013-07-06 18:11:00 UTC (rev 553) +++ trunk/src/GRHydro_Prim2Con.F90 2013-07-06 18:11:04 UTC (rev 554) @@ -378,6 +378,7 @@ ! away), we use the average temperature between cells and call ! the EOS with keytemp=1 lkeytemp=1 + temp = temp0 call EOS_Omni_press(handle,lkeytemp,GRHydro_eos_rf_prec,n,& drho,deps,temp,ye,dpress,keyerr,anyerr) lkeytemp=0 File [modified]: GRHydro_Prim2ConM.F90 Delta lines: +1 -0 =================================================================== --- trunk/src/GRHydro_Prim2ConM.F90 2013-07-06 18:11:00 UTC (rev 553) +++ trunk/src/GRHydro_Prim2ConM.F90 2013-07-06 18:11:04 UTC (rev 554) @@ -294,6 +294,7 @@ ! away), we use the average temperature between cells and call ! the EOS with keytemp=1 keytemp=1 + temp = temp0 call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& drho,deps,temp,ye,dpress,keyerr,anyerr) keytemp=0

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 553)
by rhaas＠tapir.caltech.edu 06 Jul '13

06 Jul '13

User: rhaas Date: 2013/07/06 01:11 PM Modified: /trunk/src/ GRHydro_PPM.F90, GRHydro_Source.F90 Log: GRHydro: write floating point operations in left/right symmetric manner * write DIFF_X_4 in an explicitly symmetric manner this ensures that the order of evaluation is such that symmetry in i+/-1 is preserved by the numerical derivative * add code to make ePPM and Source more symmetric this is an attempt to better preserve the CoM in long runs by eliminating asymmetries between left and right. This is fragile in that it does not make the code more robust against existing roundoff-level asymmetries, it only prevents them from developping. From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.F90 Delta lines: +87 -89 =================================================================== --- trunk/src/GRHydro_PPM.F90 2013-07-06 18:10:57 UTC (rev 552) +++ trunk/src/GRHydro_PPM.F90 2013-07-06 18:11:00 UTC (rev 553) @@ -145,7 +145,7 @@ dvely(i) = 0.5d0 * (vy(i+1) - vy(i-1)) dvelz(i) = 0.5d0 * (vz(i+1) - vz(i-1)) dpress(i) = press(i+1) - press(i-1) - d2rho(i) = (rho(i+1) - 2.d0 * rho(i) + rho(i-1))! / 6.d0 / dx / dx + d2rho(i) = ((rho(i+1) + rho(i-1)) - 2.d0 * rho(i))! / 6.d0 / dx / dx ! since we use d2rho only for the condition d2rho(i+1)*d2rhoi(i-1)<0 ! the denominator is not necessary end do @@ -209,9 +209,9 @@ if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then &&\ alim = ah &&\ else &&\ - D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ - D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ - D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ + D2a = 3.0d0 * ((a(i) + a(i+1)) - 2.0d0*ah ) &&\ + D2aL = ((a(i-1) + a(i+1)) - 2.0d0*a(i) ) &&\ + D2aR = ((a(i) + a(i+2)) - 2.0d0*a(i+1)) &&\ D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then &&\ alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ @@ -224,9 +224,9 @@ if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then &&\ alim = ah &&\ else &&\ - D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ - D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ - D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ + D2a = 3.0d0 * ((a(i) + a(i+1)) - 2.0d0*ah ) &&\ + D2aL = ((a(i-1) + a(i+1)) - 2.0d0*a(i) ) &&\ + D2aR = ((a(i) + a(i+2)) - 2.0d0*a(i+1)) &&\ D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then &&\ alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ @@ -319,7 +319,7 @@ do i = 3, nx - 2 if ( (d2rho(i+1)*d2rho(i-1) < 0.d0).and.(abs(rho(i+1)-rho(i-1)) - & ppm_epsilon_shock * min(abs(rho(i+1)), abs(rho(i-1))) > 0.d0) ) then - etatilde = (rho(i-2) - rho(i+2) + 4.d0 * drho(i)) / (drho(i) * 12.d0) + etatilde = ((rho(i-2) - rho(i+2)) + 4.d0 * drho(i)) / (drho(i) * 12.d0) else etatilde = 0.d0 end if @@ -540,12 +540,12 @@ daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) &&\ + D2aC = (a(i-1) + a(i+1)) - 2.0d0*a(i) &&\ + D2aL = (a(i-2) + a(i) ) - 2.0d0*a(i-1) &&\ + D2aLL = (a(i-3) + a(i-1)) - 2.0d0*a(i-2) &&\ + D2aR = (a(i) + a(i+2)) - 2.0d0*a(i+1) &&\ + D2aRR = (a(i+1) + a(i+3)) - 2.0d0*a(i+2) &&\ D3a = D2aR - D2aC &&\ D3aL = D2aC - D2aL &&\ D3aR = D2aRR - D2aR &&\ @@ -568,9 +568,9 @@ aplus(i) = a(i) + daplus * rhi &&\ aminus(i) = a(i) - daminus * rhi &&\ else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + aminus(i) = a(i) - (2.0d0*(1.0d0-rhi)*daplus + rhi*daminus) &&\ else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + aplus(i) = a(i) + (2.0d0*(1.0d0-rhi)*daminus + rhi*daplus) &&\ endif &&\ endif &&\ endif &&\ @@ -594,12 +594,12 @@ daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) &&\ + D2aC = (a(i-1) + a(i+1)) - 2.0d0*a(i) &&\ + D2aL = (a(i-2) + a(i) ) - 2.0d0*a(i-1) &&\ + D2aLL = (a(i-3) + a(i-1)) - 2.0d0*a(i-2) &&\ + D2aR = (a(i) + a(i+2)) - 2.0d0*a(i+1) &&\ + D2aRR = (a(i+1) + a(i+3)) - 2.0d0*a(i+2) &&\ D3a = D2aR - D2aC &&\ D3aL = D2aC - D2aL &&\ D3aR = D2aRR - D2aR &&\ @@ -623,9 +623,9 @@ aplus(i) = a(i) + daplus * rhi &&\ aminus(i) = a(i) - daminus * rhi &&\ else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + aminus(i) = a(i) - (2.0d0*(1.0d0-rhi)*daplus + rhi*daminus) &&\ else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + aplus(i) = a(i) + (2.0d0*(1.0d0-rhi)*daminus + rhi*daplus) &&\ endif &&\ else &&\ if (daplus*daminus .lt. 0) then &&\ @@ -667,10 +667,10 @@ daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) &&\ + D2aC = (a(i-1) + a(i+1)) - 2.0d0*a(i) &&\ + D2aL = (a(i-2) + a(i) ) - 2.0d0*a(i-1) &&\ + D2aR = (a(i) + a(i+2)) - 2.0d0*a(i+1) &&\ if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ else &&\ @@ -686,9 +686,9 @@ aplus(i) = a(i) + daplus * rhi &&\ aminus(i) = a(i) - daminus * rhi &&\ else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + aminus(i) = a(i) - (2.0d0*(1.0d0-rhi)*daplus + rhi*daminus) &&\ else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + aplus(i) = a(i) + (2.0d0*(1.0d0-rhi)*daminus + rhi*daplus) &&\ endif &&\ endif &&\ trivial_rp(i-1) = .false. &&\ @@ -712,10 +712,10 @@ daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = (a(i-1) + a(i+1)) - 2.0d0*a(i) &&\ + D2aL = (a(i-2) + a(i) ) - 2.0d0*a(i-1) &&\ + D2aR = (a(i) + a(i+2)) - 2.0d0*a(i+1) &&\ if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ else &&\ @@ -732,9 +732,9 @@ aplus(i) = a(i) + daplus * rhi &&\ aminus(i) = a(i) - daminus * rhi &&\ else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + aminus(i) = a(i) - (2.0d0*(1.0d0-rhi)*daplus + rhi*daminus) &&\ else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + aplus(i) = a(i) + (2.0d0*(1.0d0-rhi)*daminus + rhi*daplus) &&\ endif &&\ else &&\ if (daplus*daminus .lt. 0) then &&\ @@ -1045,7 +1045,6 @@ logical :: cond - #define STEEP(x,dx,dmx) \ if ( (x(i+1) - x(i)) * (x(i) - x(i-1)) > 0.d0 ) then &&\ dmx(i) = sign(one, dx(i)) * \ @@ -1099,9 +1098,9 @@ if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then &&\ alim = ah &&\ else &&\ - D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ - D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ - D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ + D2a = 3.0d0 * ((a(i) + a(i+1)) - 2.0d0*ah ) &&\ + D2aL = ((a(i-1) + a(i+1)) - 2.0d0*a(i) ) &&\ + D2aR = ((a(i) + a(i+2)) - 2.0d0*a(i+1)) &&\ D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then &&\ alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ @@ -1114,9 +1113,9 @@ if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then &&\ alim = ah &&\ else &&\ - D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ - D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ - D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ + D2a = 3.0d0 * ((a(i) + a(i+1)) - 2.0d0*ah ) &&\ + D2aL = ((a(i-1) + a(i+1)) - 2.0d0*a(i) ) &&\ + D2aR = ((a(i) + a(i+2)) - 2.0d0*a(i+1)) &&\ D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then &&\ alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ @@ -1221,12 +1220,12 @@ daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) &&\ + D2aC = (a(i-1) + a(i+1)) - 2.0d0*a(i) &&\ + D2aL = (a(i-2) + a(i) ) - 2.0d0*a(i-1) &&\ + D2aLL = (a(i-3) + a(i-1)) - 2.0d0*a(i-2) &&\ + D2aR = (a(i) + a(i+2)) - 2.0d0*a(i+1) &&\ + D2aRR = (a(i+1) + a(i+3)) - 2.0d0*a(i+2) &&\ D3a = D2aR - D2aC &&\ D3aL = D2aC - D2aL &&\ D3aR = D2aRR - D2aR &&\ @@ -1249,9 +1248,9 @@ aplus(i) = a(i) + daplus * rhi &&\ aminus(i) = a(i) - daminus * rhi &&\ else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + aminus(i) = a(i) - (2.0d0*(1.0d0-rhi)*daplus + rhi*daminus) &&\ else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + aplus(i) = a(i) + (2.0d0*(1.0d0-rhi)*daminus + rhi*daplus) &&\ endif &&\ endif &&\ endif &&\ @@ -1271,12 +1270,12 @@ daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) &&\ + D2aC = (a(i-1) + a(i+1)) - 2.0d0*a(i) &&\ + D2aL = (a(i-2) + a(i) ) - 2.0d0*a(i-1) &&\ + D2aLL = (a(i-3) + a(i-1)) - 2.0d0*a(i-2) &&\ + D2aR = (a(i) + a(i+2)) - 2.0d0*a(i+1) &&\ + D2aRR = (a(i+1) + a(i+3)) - 2.0d0*a(i+2) &&\ D3a = D2aR - D2aC &&\ D3aL = D2aC - D2aL &&\ D3aR = D2aRR - D2aR &&\ @@ -1300,9 +1299,9 @@ aplus(i) = a(i) + daplus * rhi &&\ aminus(i) = a(i) - daminus * rhi &&\ else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + aminus(i) = a(i) - (2.0d0*(1.0d0-rhi)*daplus + rhi*daminus) &&\ else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + aplus(i) = a(i) + (2.0d0*(1.0d0-rhi)*daminus + rhi*daplus) &&\ endif &&\ else &&\ if (daplus*daminus .lt. 0) then &&\ @@ -1339,10 +1338,10 @@ daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) &&\ + D2aC = (a(i-1) + a(i+1)) - 2.0d0*a(i) &&\ + D2aL = (a(i-2) + a(i) ) - 2.0d0*a(i-1) &&\ + D2aR = (a(i) + a(i+2)) - 2.0d0*a(i+1) &&\ if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ else &&\ @@ -1358,9 +1357,9 @@ aplus(i) = a(i) + daplus * rhi &&\ aminus(i) = a(i) - daminus * rhi &&\ else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + aminus(i) = a(i) - (2.0d0*(1.0d0-rhi)*daplus + rhi*daminus) &&\ else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + aplus(i) = a(i) + (2.0d0*(1.0d0-rhi)*daminus + rhi*daplus) &&\ endif &&\ endif &&\ else &&\ @@ -1380,10 +1379,10 @@ daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) &&\ + D2aC = (a(i-1) + a(i+1)) - 2.0d0*a(i) &&\ + D2aL = (a(i-2) + a(i) ) - 2.0d0*a(i-1) &&\ + D2aR = (a(i) + a(i+2)) - 2.0d0*a(i+1) &&\ if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ else &&\ @@ -1400,9 +1399,9 @@ aplus(i) = a(i) + daplus * rhi &&\ aminus(i) = a(i) - daminus * rhi &&\ else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + aminus(i) = a(i) - (2.0d0*(1.0d0-rhi)*daplus + rhi*daminus) &&\ else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + aplus(i) = a(i) + (2.0d0*(1.0d0-rhi)*daminus + rhi*daplus) &&\ endif &&\ else &&\ if (daplus*daminus .lt. 0) then &&\ @@ -1772,7 +1771,6 @@ CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax - do i = 2, nx - 1 dpress(i) = press(i+1) - press(i-1) end do @@ -1836,9 +1834,9 @@ if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then &&\ alim = ah &&\ else &&\ - D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ - D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ - D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ + D2a = 3.0d0 * ((a(i) + a(i+1)) - 2.0d0*ah ) &&\ + D2aL = ((a(i-1) + a(i+1)) - 2.0d0*a(i) ) &&\ + D2aR = ((a(i) + a(i+2)) - 2.0d0*a(i+1)) &&\ if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then &&\ alim = 0.5d0*(a(i)+a(i+1)) - sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ else &&\ @@ -1928,12 +1926,12 @@ daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) &&\ + D2aC = (a(i-1) + a(i+1)) - 2.0d0*a(i) &&\ + D2aL = (a(i-2) + a(i) ) - 2.0d0*a(i-1) &&\ + D2aLL = (a(i-3) + a(i-1)) - 2.0d0*a(i-2) &&\ + D2aR = (a(i) + a(i+2)) - 2.0d0*a(i+1) &&\ + D2aRR = (a(i+1) + a(i+3)) - 2.0d0*a(i+2) &&\ D3a = D2aR - D2aC &&\ D3aL = D2aC - D2aL &&\ D3aR = D2aRR - D2aR &&\ @@ -1956,9 +1954,9 @@ aplus(i) = a(i) + daplus * rhi &&\ aminus(i) = a(i) - daminus * rhi &&\ else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + aminus(i) = a(i) - (2.0d0*(1.0d0-rhi)*daplus + rhi*daminus) &&\ else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + aplus(i) = a(i) + (2.0d0*(1.0d0-rhi)*daminus + rhi*daplus) &&\ endif &&\ endif &&\ endif &&\ @@ -1981,10 +1979,10 @@ daplus = aplus(i)-a(i) &&\ daminus = a(i)-aminus(i) &&\ if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ - D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ - D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ - D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ - D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) &&\ + D2aC = (a(i-1) + a(i+1)) - 2.0d0*a(i) &&\ + D2aL = (a(i-2) + a(i) ) - 2.0d0*a(i-1) &&\ + D2aR = (a(i) + a(i+2)) - 2.0d0*a(i+1) &&\ if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ else &&\ @@ -2000,9 +1998,9 @@ aplus(i) = a(i) + daplus * rhi &&\ aminus(i) = a(i) - daminus * rhi &&\ else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + aminus(i) = a(i) - (2.0d0*(1.0d0-rhi)*daplus + rhi*daminus) &&\ else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + aplus(i) = a(i) + (2.0d0*(1.0d0-rhi)*daminus + rhi*daplus) &&\ endif &&\ endif &&\ else &&\ File [modified]: GRHydro_Source.F90 Delta lines: +6 -0 =================================================================== --- trunk/src/GRHydro_Source.F90 2013-07-06 18:10:57 UTC (rev 552) +++ trunk/src/GRHydro_Source.F90 2013-07-06 18:11:00 UTC (rev 553) @@ -15,12 +15,18 @@ ! Fourth order f.d. +#if(1) +#define DIFF_X_4(q) ((q(i-2,j,k)-q(i+2,j,k)) + 8.d0 * (q(i+1,j,k) - q(i-1,j,k))) / 12.d0 * ida +#define DIFF_Y_4(q) ((q(i,j-2,k)-q(i,j+2,k)) + 8.d0 * (q(i,j+1,k) - q(i,j-1,k))) / 12.d0 * idb +#define DIFF_Z_4(q) ((q(i,j,k-2)-q(i,j,k+2)) + 8.d0 * (q(i,j,k+1) - q(i,j,k-1))) / 12.d0 * idc +#else #define DIFF_X_4(q) (-q(i+2,j,k) + 8.d0 * q(i+1,j,k) - 8.d0 * q(i-1,j,k) + \ q(i-2,j,k)) / 12.d0 * ida #define DIFF_Y_4(q) (-q(i,j+2,k) + 8.d0 * q(i,j+1,k) - 8.d0 * q(i,j-1,k) + \ q(i,j-2,k)) / 12.d0 * idb #define DIFF_Z_4(q) (-q(i,j,k+2) + 8.d0 * q(i,j,k+1) - 8.d0 * q(i,j,k-1) + \ q(i,j,k-2)) / 12.d0 * idc +#endif #include "cctk.h" #include "cctk_Parameters.h"

1 0

← Newer
1
2
3
4
5
6
7
8
9
Older →

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits July 2013