Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 197)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:48 PM Modified: /trunk/ param.ccl, schedule.ccl /trunk/src/ GRHydro_BondiM.c Log: GRHydro: add parameters to set Bondi solution only in certain range. From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +17 -0 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-01-11 15:04:15 UTC (rev 196) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:53 UTC (rev 197) @@ -84,6 +84,7 @@ //some math helpers #define SQR(x) ((x)*(x)) +#define CUBE(x) ((x)*(x)*(x)) static CCTK_REAL Mdot, rs, vs_sq, vs, cs_sq, cs, rhos, hs, K, Qdot, gamma_eos, r_sol; @@ -902,6 +903,19 @@ ***********************************************************************************/ void GRHydro_BondiM(CCTK_ARGUMENTS) { + GRHydro_BondiM_Internal(CCTK_PASS_CTOC, 1e100, -1e100); +} + +void GRHydro_BondiM_Range(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_PARAMETERS; + GRHydro_BondiM_Internal(CCTK_PASS_CTOC, bondi_freeze_inner_radius, + bondi_freeze_outer_radius); +} + + +void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max) +{ DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; CCTK_REAL det; @@ -1134,6 +1148,9 @@ rspher = x_spher[RR]; + if(rspher > range_min && rspher < range_max) + continue; + /* Find nearest point in the Bondi solution : */ j = (int) ( 0.5 + (log10(rspher) - logrmin) / dlogr ) ; Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +14 -0 =================================================================== --- trunk/param.ccl 2013-01-11 15:04:15 UTC (rev 196) +++ trunk/param.ccl 2013-03-28 01:48:53 UTC (rev 197) @@ -395,6 +395,20 @@ (0:* :: "positive" } 1.0 +CCTK_BOOLEAN bondi_evolve_only_annulus "reset to initial data outside of bondi_freeze_inner_radius and bondi_freeze_outer_radius" +{ +} "no" + +CCTK_REAL bondi_freeze_inner_radius "reset to initial at radii below this" +{ + *:* :: "any value" +} -1. + +CCTK_REAL bondi_freeze_outer_radius "reset to initial at radii above this" +{ + *:* :: "any value" +} 1e300 + # For Poloidal Magnetic field test: CCTK_REAL poloidal_A_b "Vector potential strength" File [modified]: schedule.ccl Delta lines: +7 -0 =================================================================== --- trunk/schedule.ccl 2013-01-11 15:04:15 UTC (rev 196) +++ trunk/schedule.ccl 2013-03-28 01:48:53 UTC (rev 197) @@ -250,6 +250,13 @@ } "setup GRHydro vars for the magnetized Bondi solution" } +if(bondi_evolve_only_annulus) { +SCHEDULE GRHydro_BondiM_Range IN HydroBase_Con2Prim BEFORE Con2Prim +{ + LANG: C +} "force analytic solution outside anulus" +} + 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

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 501)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:47 PM Modified: /trunk/src/ GRHydro_TVDReconstruct_drv.F90 Log: GRHydro: support reconstruct_Wv in TVD reconstructor From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_TVDReconstruct_drv.F90 Delta lines: +107 -8 =================================================================== --- trunk/src/GRHydro_TVDReconstruct_drv.F90 2013-03-28 01:47:14 UTC (rev 500) +++ trunk/src/GRHydro_TVDReconstruct_drv.F90 2013-03-28 01:47:16 UTC (rev 501) @@ -50,6 +50,8 @@ ! 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 vel, lvel TARGET Bvec, lBvec @@ -67,14 +69,32 @@ CCTK_INT :: ierr + ! variables used when reconstruct_Wv == true + CCTK_REAL :: inv_w ! 1/Lorentz factor + CCTK_REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: wvel ! vel*w_lorentz + CCTK_REAL :: agxx,agxy,agxz,agyy,agyz,agzz ! metric components + ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc vup => lvel Bprim => lBvec else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz vup => vel Bprim => Bvec end if @@ -83,6 +103,29 @@ if (ierr .ne. 0) then call CCTK_WARN(0, "Allocation problems with trivial_rp") end if + + !!$ reconstruct w^i = w_lorentz*v^i to ensure slower than light + !!$ speeds? + if (reconstruct_Wv .ne. 0) then + ! all velocity like quantities are now w_lorentz*velocity. We will convert + ! back to ordinary velocity at the end, using w_lorentz = sqrt(1 + g_{ij} + ! w^i w^j). + allocate(wvel(cctk_ash(1),cctk_ash(2),cctk_ash(3),3),STAT=ierr) + if (ierr .ne. 0) then + call CCTK_WARN(0, "Allocation problems with wvel") + end if + !$OMP PARALLEL DO PRIVATE(i,j,k) + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + wvel(i,j,k,1) = vup(i,j,k,1) * w_lorentz(i,j,k) + wvel(i,j,k,2) = vup(i,j,k,2) * w_lorentz(i,j,k) + wvel(i,j,k,3) = vup(i,j,k,3) * w_lorentz(i,j,k) + end do + end do + end do + !$OMP END PARALLEL DO + end if call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX") call SpaceMask_GetStateBits(trivialx, "Hydro_RiemannProblemX", & @@ -184,14 +227,69 @@ call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vup(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vup(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vup(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) - call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - eps, epsplus, epsminus, trivial_rp, hydro_excision_mask) + if (reconstruct_Wv.ne.0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + wvel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + wvel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + wvel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) + ! divide out the Loretnz factor obtained from w_lorentz = + ! sqrt(1+g_{ij} w^i w^j) for both the + ! plus and minus quantities this should by construction ensure + ! that any Lorentz factor calculated from them later on is + ! physical (ie. > 1.d0) + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces + !$OMP PARALLEL DO PRIVATE(i,j,k,inv_w,agxx,agxy,agxz,agyy,agyz,agzz) + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + transport_constraints*(1-xoffset) + agxx = 0.5d0*( g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset) ) + agxy = 0.5d0*( g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset) ) + agxz = 0.5d0*( g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset) ) + agyy = 0.5d0*( g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset) ) + agyz = 0.5d0*( g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset) ) + agzz = 0.5d0*( g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset) ) + inv_w = 1d0/sqrt( 1.d0 + agxx*velxminus(i,j,k)*velxminus(i,j,k) & + + agyy*velyminus(i,j,k)*velyminus(i,j,k) & + + agzz*velzminus(i,j,k)*velzminus(i,j,k) & + + 2.d0*agxy*velxminus(i,j,k)*velyminus(i,j,k) & + + 2.d0*agxz*velxminus(i,j,k)*velzminus(i,j,k) & + + 2.d0*agyz*velyminus(i,j,k)*velzminus(i,j,k) ) + velxminus(i,j,k) = velxminus(i,j,k)*inv_w + velyminus(i,j,k) = velyminus(i,j,k)*inv_w + velzminus(i,j,k) = velzminus(i,j,k)*inv_w + + agxx = 0.5d0*( g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset) ) + agxy = 0.5d0*( g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset) ) + agxz = 0.5d0*( g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset) ) + agyy = 0.5d0*( g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset) ) + agyz = 0.5d0*( g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset) ) + agzz = 0.5d0*( g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset) ) + inv_w = 1d0/sqrt( 1.d0 + agxx*velxplus(i,j,k)*velxplus(i,j,k) & + + agyy*velyplus(i,j,k)*velyplus(i,j,k) & + + agzz*velzplus(i,j,k)*velzplus(i,j,k) & + + 2.d0*agxy*velxplus(i,j,k)*velyplus(i,j,k) & + + 2.d0*agxz*velxplus(i,j,k)*velzplus(i,j,k) & + + 2.d0*agyz*velyplus(i,j,k)*velzplus(i,j,k) ) + velxplus(i,j,k) = velxplus(i,j,k)*inv_w + velyplus(i,j,k) = velyplus(i,j,k)*inv_w + velzplus(i,j,k) = velzplus(i,j,k)*inv_w + end do + end do + end do + !$OMP END PARALLEL DO + else + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + vup(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + vup(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + vup(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) + end if + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + eps, epsplus, epsminus, trivial_rp, hydro_excision_mask) if(evolve_mhd.ne.0) then call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & Bprim(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask) @@ -270,6 +368,7 @@ !!$ TVD ends. deallocate(trivial_rp) + ! Fortran 90 will automatically deallocate wvel for us end subroutine GRHydro_TVDReconstruct_drv

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 500)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:47 PM Added: /trunk/src/ GRHydro_Con2PrimM_ptee.c Modified: /trunk/ interface.ccl /trunk/src/ GRHydro_Con2PrimM.F90, make.code.defn Log: Finish the implementation of con2prim using the entropy equation. * Tested successfully pointwisely. It still needs to be tested on evolution. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +34 -10 =================================================================== --- trunk/src/GRHydro_Con2PrimM.F90 2013-03-28 01:47:09 UTC (rev 499) +++ trunk/src/GRHydro_Con2PrimM.F90 2013-03-28 01:47:14 UTC (rev 500) @@ -334,7 +334,11 @@ Bvecz_tmp = Bprim(i,j,k,3) keytemp = 0 - + !Watch out for the values returned to b2. Here b2 is the Bprim^2 + !while inside the point-wise con2prim routines it is the square + !of the comoving B-field, b^{\mu} b_{\mu}. It is overwritten + !in this routine, but we may need to find a better notation + !avoid future confusions. call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & @@ -347,7 +351,13 @@ epsnegative,GRHydro_C2P_failed(i,j,k)) if(evolve_entropy.ne.0) then - entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho(i,j,k) + if(GRHydro_C2P_failed(i,j,k).ne.0) then + !Use previous time step for rho: + entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho(i,j,k) + else + !Use the current correct value of rho returned by con2prim: + entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho_tmp + endif endif if(GRHydro_C2P_failed(i,j,k).ne.0) then @@ -430,14 +440,28 @@ Bvecy_tmp = Bprim(i,j,k,2) Bvecz_tmp = Bprim(i,j,k,3) - call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, & - dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & - Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3),rho_tmp,& - velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& - Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& - g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & - epsnegative,GRHydro_C2P_failed(i,j,k)) + if(evolve_entropy.ne.0) then + call GRHydro_Con2PrimM_ptee(GRHydro_eos_handle, keytemp, & + GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & + scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & + Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), & + entropycons(i,j,k), xye(1), & + xtemp(1),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& + eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& + w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + else + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, & + dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & + Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3),rho_tmp,& + velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& + Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + end if rho(i,j,k) = rho_tmp press(i,j,k) = press_tmp File [added]: GRHydro_Con2PrimM_ptee.c Delta lines: +592 -0 =================================================================== --- trunk/src/GRHydro_Con2PrimM_ptee.c (rev 0) +++ trunk/src/GRHydro_Con2PrimM_ptee.c 2013-03-28 01:47:14 UTC (rev 500) @@ -0,0 +1,592 @@ +/*********************************************************************************** + Copyright 2006 Scott C. Noble, Charles F. Gammie, Jonathan C. McKinney, + and Luca Del Zanna. + + PVS_GRMHD + + This file was derived from PVS_GRMHD. The authors of PVS_GRMHD include + Scott C. Noble, Charles F. Gammie, Jonathan C. McKinney, and Luca Del Zanna. + PVS_GRMHD is available under the GPL from: + http://rainman.astro.uiuc.edu/codelib/ + + You are morally obligated to cite the following paper in his/her + scientific literature that results from use of this file: + + [1] Noble, S. C., Gammie, C. F., McKinney, J. C., \& Del Zanna, L. \ 2006, + Astrophysical Journal, 641, 626. + + PVS_GRMHD is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + PVS_GRMHD is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with PVS_GRMHD; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + + If the user has any questions, please direct them to Scott C. Noble at + scn(a)astro.rit.edu . + +***********************************************************************************/ + + + +#include <stdio.h> +#include <stdlib.h> +#include <stdarg.h> +#include <string.h> +#include <math.h> +#include <float.h> +#include <complex.h> + +#include "cctk.h" +#include "cctk_Parameters.h" + +/* Set this to be 1 if you want debug output */ +#define DEBUG_CON2PRIMM (0) + + +/* Adiabatic index used for the state equation */ + +#define MAX_NEWT_ITER (30) /* Max. # of Newton-Raphson iterations for find_root_2D(); */ +#define NEWT_TOL (1.0e-10) /* Min. of tolerance allowed for Newton-Raphson iterations */ +#define MIN_NEWT_TOL (1.0e-10) /* Max. of tolerance allowed for Newton-Raphson iterations */ +#define EXTRA_NEWT_ITER (2) + +#define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ +#define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ + +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed + to always be smaller than this. This + is used to detect solver failures */ + +#define FAIL_VAL (1.e30) /* Generic value to which we set variables when a problem arises */ + +/************************************************** + The following functions assume a Gamma-law EOS: +***************************************************/ + +/* Local Globals */ +struct LocGlob { + CCTK_REAL Bsq, QdotBsq, Qtsq, Qdotn, D, half_Bsq, Sc, g_o_gm1, + W_for_gnr2, rho_for_gnr2, W_for_gnr2_old, rho_for_gnr2_old, drho_dW ; +} ; + +// Declarations: + + +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( + CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, + CCTK_REAL *gamma_eos, + CCTK_REAL *dens_in, + CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, + CCTK_REAL *tau_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, + CCTK_REAL *entropycons_in, + CCTK_REAL *y_e_in, CCTK_REAL* temp_in, + CCTK_REAL *rho, + CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, + CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, + CCTK_REAL *w_lorentz, + CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, + CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, + CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, + CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, + CCTK_REAL *det, + CCTK_INT *epsnegative, + CCTK_REAL *retval) ; + +static CCTK_INT general_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, + struct LocGlob *lgp, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); + +static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, + CCTK_REAL gammaeos, struct LocGlob *lgp); + +/**********************************************************************/ +/********************************************************************************** + + Con2PrimM_pt(): + ----------------------------- + + -- Attempts an inversion from GRMHD conserved variables to primitive variables assuming a guess. + + -- Uses the 2D method of Noble et al. (2006): + -- Solves for two independent variables (W,v^2) via a 2D + Newton-Raphson method + -- Can be used (in principle) with a general equation of state. + + -- Minimizes two residual functions using a homemade Newton-Raphson routine. + -- It is homemade so that it can catch exceptions and handle them correctly, plus it is + optimized for this particular problem. + + -- Note that the notation used herein is that of Noble et al. (2006) except for the argument + list. + + +INPUT: (using GRHydro variable defintions) + + s[x,y,z] = scons[0,1,2] = \alpha \sqrt(\gamma) T^0_i + dens, tau = as defined in GRHydro and are assumed to be densitized (i.e. with sqrt(\gamma)) + dens = D = \sqrt(\gamma) W \rho + tau = \alpha^2 \sqrt(\gamma) T^{00} - D + g[x,y,z][x,y,x] = spatial metric corresponding to \gamma + u[x,y,z][x,y,z] = inverse of the spatial metric, g[x,y,z][x,y,x] + det = sqrt(\gamma) + B[x,y,z] = Bvec[0,1,2] + bsq = b^\mu b_\mu + + epsnegative = (integer) + = 0 if rho and epsilon are positive + != 0 otherwise + + + -- (users should set B[x,y,z] = 0 for hydrodynamic runs) + + +OUTPUT: (using GRHydro variable defintions) + rho, eps = as defined in GRHydro, primitive variables + vel[x,y,z] = as defined in GRHydro, primitive variables + + +RETURN VALUE: of retval = (i*100 + j) where + i = 0 -> Newton-Raphson solver either was not called (yet or not used) + or returned successfully; + 1 -> Newton-Raphson solver did not converge to a solution with the + given tolerances; + 2 -> Newton-Raphson procedure encountered a numerical divergence + (occurrence of "nan" or "+/-inf" ; + + j = 0 -> success + 1 -> failure: some sort of failure in Newton-Raphson; + 2 -> failure: unphysical vsq = v^2 value at initial guess; + 3 -> failure: W<0 or W>W_TOO_BIG + 4 -> failure: v^2 > 1 + ( used to be 5 -> failure: rho,uu <= 0 but now sets epsnegative to non-zero ) + +**********************************************************************************/ + +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( + CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, + CCTK_REAL *gamma_eos, + CCTK_REAL *dens_in, + CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, + CCTK_REAL *tau_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, + CCTK_REAL *entropycons_in, + CCTK_REAL *y_e_in, CCTK_REAL* temp_in, + CCTK_REAL *rho, + CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, + CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, + CCTK_REAL *w_lorentz, + CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, + CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, + CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, + CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, + CCTK_REAL *det, + CCTK_INT *epsnegative, + CCTK_REAL *retval) + +{ + CCTK_REAL x_1d[1]; + CCTK_REAL sx, sy, sz; + CCTK_REAL usx, usy, usz; + CCTK_REAL tau, dens, gammaeos; + CCTK_REAL QdotB; + CCTK_REAL rho0,u,p,w,gammasq,gamma,W,vsq; + CCTK_REAL g_o_WBsq, QdB_o_W; + CCTK_REAL detg = (*det); + CCTK_REAL sqrt_detg = sqrt(detg); + CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; + CCTK_REAL rho_gm1; + CCTK_REAL utsq; + + DECLARE_CCTK_PARAMETERS; + + struct LocGlob lg; + + gammaeos = *gamma_eos; + + /* Assume ok initially: */ + // BCM: But let the driver function take care of its initialization + //*retval = 0.; + *epsnegative = 0; + +#if(DEBUG_CON2PRIMM) + fprintf(stdout," *dens = %26.16e \n", *dens_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); + fprintf(stdout," *entropycons = %26.16e \n", *entropycons_in ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + fprintf(stdout," *vely = %26.16e \n", *vely ); + fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); + fprintf(stdout," *pressure = %26.16e \n", *pressure ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); + fprintf(stdout," *retval = %26.16e \n", *retval ); + fflush(stdout); +#endif + + /* First undensitize all conserved variables : */ + sx = ( *sx_in) * inv_sqrt_detg; + sy = ( *sy_in) * inv_sqrt_detg; + sz = ( *sz_in) * inv_sqrt_detg; + tau = ( *tau_in) * inv_sqrt_detg; + dens = (*dens_in) * inv_sqrt_detg; + + usx = (*uxx)*sx + (*uxy)*sy + (*uxz)*sz; + usy = (*uxy)*sx + (*uyy)*sy + (*uyz)*sz; + usz = (*uxz)*sx + (*uyz)*sy + (*uzz)*sz; + + *Bx = (*Bconsx_in) * inv_sqrt_detg; + *By = (*Bconsy_in) * inv_sqrt_detg; + *Bz = (*Bconsz_in) * inv_sqrt_detg; + + lg.Sc = (*entropycons_in) * inv_sqrt_detg; + // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: + + lg.g_o_gm1 = gammaeos/(gammaeos-1.0); + + lg.Bsq = + (*gxx) * (*Bx) * (*Bx) + + (*gyy) * (*By) * (*By) + + (*gzz) * (*Bz) * (*Bz) + + 2*( + (*gxy) * (*Bx) * (*By) + + (*gxz) * (*Bx) * (*Bz) + + (*gyz) * (*By) * (*Bz) ); + + QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ; + lg.QdotBsq = QdotB*QdotB ; + + lg.Qdotn = -(tau + dens) ; + + lg.Qtsq = (usx * sx + usy * sy + usz * sz) ; + + lg.D = dens; + + lg.half_Bsq = 0.5*lg.Bsq; + + /* calculate W from last timestep and use for guess */ + vsq = + (*gxx) * (*velx) * (*velx) + + (*gyy) * (*vely) * (*vely) + + (*gzz) * (*velz) * (*velz) + + 2*( + (*gxy) * (*velx) * (*vely) + + (*gxz) * (*velx) * (*velz) + + (*gyz) * (*vely) * (*velz) ); + + if( (vsq < 0.) && (fabs(vsq) < 1.0e-13) ) { + vsq = fabs(vsq); + } + if(vsq < 0. || vsq > 1. ) { + *retval = 2.; + return; + } + + gammasq = 1. / (1. - vsq); + gamma = sqrt(gammasq); + + // Always calculate rho from D and gamma so that using D in EOS remains consistent + // i.e. you don't get positive values for dP/d(vsq) . + rho0 = lg.D / gamma ; + rho_gm1 = pow(rho0,(gammaeos-1.)); + p = lg.Sc * rho_gm1 / gamma; + u = p / (gammaeos-1.); + w = rho0 + u + p ; + +// W_last = w*gammasq ; + + // Calculate W and vsq: + x_1d[0] = rho0; +// *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, &lg, func_vsq ) ; + *retval = general_newton_raphson( x_1d, gammaeos, &lg, func_rho ) ; + rho0 = x_1d[0]; + + /* Problem with solver, so return denoting error before doing anything further */ + if( ((*retval) != 0.) || (rho0 == FAIL_VAL) ) { + *retval = *retval*100.+1.; + return; + } + else{ + if( rho0 > W_TOO_BIG) { + *retval = 3.; + return; + } + } + + // Calculate v^2: + utsq = (lg.D-rho0)*(lg.D+rho0)/(rho0*rho0); + gammasq = 1.+utsq; + gamma = sqrt(gammasq); + + if( utsq < 0. ) { + *retval = 4.; + return; + } + + + // Recover the primitive variables from the scalars and conserved variables: + rho0 = lg.D / gamma; + rho_gm1 = pow(rho0,(gammaeos-1.)); + p = lg.Sc * rho_gm1 / gamma; + u = p / (gammaeos-1.); + w = rho0 + u + p ; + W = w * gammasq; + + // User may want to handle this case differently, e.g. do NOT return upon + // a negative rho/u, calculate v^i so that rho/u can be floored by other routine: + if( (rho0 <= 0.) || (u <= 0.) ) { + *epsnegative = 1; + return; + } + + *rho = rho0; + *epsilon = u / rho0; + *w_lorentz = gamma; + *pressure = p ; + + g_o_WBsq = 1./(W+lg.Bsq); + QdB_o_W = QdotB / W; + *bsq = lg.Bsq / gammasq + QdB_o_W*QdB_o_W; + + *velx = g_o_WBsq * ( usx + QdB_o_W*(*Bx) ) ; + *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ; + *velz = g_o_WBsq * ( usz + QdB_o_W*(*Bz) ) ; + + if (*rho <= rho_abs_min*(1.0+GRHydro_atmo_tolerance) ) { + *rho = rho_abs_min; + *velx = 0.0; + *vely = 0.0; + *velz = 0.0; + *w_lorentz = 1.0; + } + + +#if(DEBUG_CON2PRIMM) + fprintf(stdout,"rho = %26.16e \n",*rho ); + fprintf(stdout,"epsilon = %26.16e \n",*epsilon ); + fprintf(stdout,"pressure = %26.16e \n",*pressure ); + fprintf(stdout,"w_lorentz = %26.16e \n",*w_lorentz); + fprintf(stdout,"bsq = %26.16e \n",*bsq ); + fprintf(stdout,"velx = %26.16e \n",*velx ); + fprintf(stdout,"vely = %26.16e \n",*vely ); + fprintf(stdout,"velz = %26.16e \n",*velz ); + fprintf(stdout,"gammaeos = %26.16e \n",gammaeos ); + fflush(stdout); +#endif + + /* done! */ + return; + +} + +/**********************************************************************/ +/************************************************************ + + general_newton_raphson(): + + -- performs Newton-Rapshon method on an arbitrary system. + + -- inspired in part by Num. Rec.'s routine newt(); + + Arguements: + + -- x[] = set of independent variables to solve for; + -- n = number of independent variables and residuals; + -- funcd = name of function that calculates residuals, etc.; + +*****************************************************************/ +static CCTK_INT general_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, + struct LocGlob *lgp, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) +{ + CCTK_REAL f, df, dx[1], x_old[1], resid[1], + jac[1][1]; + CCTK_REAL errx, x_orig[1]; + CCTK_INT n_iter, i_extra, doing_extra; + CCTK_REAL W,W_old; + + int keep_iterating; + + + // Initialize various parameters and variables: + errx = 1. ; + df = f = 1.; + i_extra = doing_extra = 0; + x_old[0] = x_orig[0] = x[0] ; + + W = W_old = 0.; + + n_iter = 0; + + + /* Start the Newton-Raphson iterations : */ + keep_iterating = 1; + while( keep_iterating ) { + + (*funcd) (x, dx, resid, jac, &f, &df, gammaeos, lgp); /* returns with new dx, f, df */ + + /* Save old values before calculating the new: */ + errx = 0.; + x_old[0] = x[0] ; + /* don't use line search : */ + x[0] += dx[0] ; + + /****************************************/ + /* Calculate the convergence criterion */ + /****************************************/ + + /* For the new criterion, always look at error in "W" : */ + // METHOD specific: + errx = (x[0]==0.) ? fabs(dx[0]) : fabs(dx[0]/x[0]); + + + /****************************************/ + /* Make sure that the new x[] is physical : */ + /****************************************/ + x[0] = fabs(x[0]); + + + /*****************************************************************************/ + /* If we've reached the tolerance level, then just do a few extra iterations */ + /* before stopping */ + /*****************************************************************************/ + + if( (fabs(errx) <= NEWT_TOL) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) { + doing_extra = 1; + } + + if( doing_extra == 1 ) i_extra++ ; + + if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) || + (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + keep_iterating = 0; + } + + n_iter++; + + } // END of while(keep_iterating) + + + /* Check for bad untrapped divergences : */ + if( (!finite(f)) || (!finite(df)) || (!finite(x[0])) ) { + return(2); + } + + + if( fabs(errx) > MIN_NEWT_TOL){ + return(1); + } + if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){ + return(0); + } + if( fabs(errx) <= NEWT_TOL ){ + return(0); + } + + return(0); + +} + +/**********************************************************************/ +/********************************************************************************* + func_rho(): + + -- residual/jacobian routine to calculate rho Qtsq equation with + the definition of W + W = ( 1 + GAMMA * K_atm * rho^(GAMMA-1)/(GAMMA-1) ) D^2 / rho + substituted in. + + Arguments: + x = current value of independent var's (on input & output); + dx = Newton-Raphson step (on output); + resid = residuals based on x (on output); + jac = Jacobian matrix based on x (on output); + f = resid.resid/2 (on output) + df = -2*f; (on output) + n = dimension of x[]; + *********************************************************************************/ +static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, + CCTK_REAL gammaeos, struct LocGlob *lgp) +{ + + CCTK_REAL t1,t2; + CCTK_REAL t12; + CCTK_REAL t17; + CCTK_REAL t3, t100,rhosq,t200,rho,W,dWdrho,dvsqdrho,vsq; + + rho = x[0]; + rhosq = rho*rho; + t200 = 1./(lgp->D*lgp->D); + t100 = lgp->g_o_gm1*lgp->Sc*pow(rho,(gammaeos-1.)); + W = lgp->D*( lgp->D + t100 ) / rho ; + dWdrho = lgp->D * ( -lgp->D + t100*(gammaeos-2.) ) / rhosq; + t1 = W*W; + t2 = lgp->Bsq+W; + // t3 = pow(Bsq+W,2.0); + t3 = t2*t2; + vsq = (lgp->D-rho)*(lgp->D+rho)*t200; + dvsqdrho = -2*rho*t200; + resid[0] = t1*(lgp->Qtsq-vsq*t3)+lgp->QdotBsq*(t2+W); + t12 = lgp->Bsq*lgp->Bsq; + t17 = dWdrho*vsq; + jac[0][0] = 2*lgp->QdotBsq*dWdrho + +((lgp->Qtsq-vsq*t12)*2*dWdrho+(-6*t17*lgp->Bsq-dvsqdrho*t12 + +(-2*dvsqdrho*lgp->Bsq-4*t17-dvsqdrho*W)*W)*W)*W; + + dx[0] = -resid[0]/jac[0][0]; + *f = 0.5*resid[0]*resid[0]; + *df = -2. * (*f); + + return; + +} + + +/****************************************************************************** + END + ******************************************************************************/ + + +#undef DEBUG_CON2PRIMM File [modified]: make.code.defn Delta lines: +1 -0 =================================================================== --- trunk/src/make.code.defn 2013-03-28 01:47:09 UTC (rev 499) +++ trunk/src/make.code.defn 2013-03-28 01:47:14 UTC (rev 500) @@ -46,6 +46,7 @@ GRHydro_Tmunu.F90 \ GRHydro_Con2PrimM.F90 \ GRHydro_Con2PrimM_pt.c \ + GRHydro_Con2PrimM_ptee.c \ GRHydro_Con2PrimM_polytype_pt.c \ GRHydro_EigenproblemM.F90 \ GRHydro_FluxM.F90 \ Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +22 -0 =================================================================== --- trunk/interface.ccl 2013-03-28 01:47:09 UTC (rev 499) +++ trunk/interface.ccl 2013-03-28 01:47:14 UTC (rev 500) @@ -84,6 +84,27 @@ 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, \ + 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 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 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 Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT 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, \ @@ -207,6 +228,7 @@ #PROVIDES FUNCTION Con2Prim WITH Con2Prim_pt LANGUAGE Fortran PROVIDES FUNCTION Con2PrimPoly WITH Con2Prim_ptPolytype LANGUAGE Fortran PROVIDES FUNCTION Con2PrimGenM 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

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 499)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:47 PM Modified: /trunk/src/ GRHydro_Prim2ConM.F90 Log: GRHydro: move critical section in prim2conM into if statement otherwise the code will enter and leave the critical section during each loop iteration From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2ConM.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/GRHydro_Prim2ConM.F90 2013-03-28 01:47:00 UTC (rev 498) +++ trunk/src/GRHydro_Prim2ConM.F90 2013-03-28 01:47:09 UTC (rev 499) @@ -159,12 +159,12 @@ xtemp(1) = 0.5d0*(temperature(i,j,k) + & temperature(i-xoffset,j-yoffset,k-zoffset)) - !$OMP CRITICAL if (y_e_minus(i,j,k) .le. 0.0d0 .or. y_e_plus(i,j,k) .le. 0.0d0) then + !$OMP CRITICAL write(NaN_WarnLine,'(a100,7g15.6)') '(y_e_minus,y_e_plus,x,y,z,rho)', y_e(i,j,k), y_e_minus(i,j,k), y_e_plus(i,j,k), x(i,j,k),y(i,j,k),z(i,j,k),rho(i,j,k) call CCTK_WARN(1, NaN_WarnLine) + !$OMP END CRITICAL endif - !$OMP END CRITICAL call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & avg_detl,densminus(i,j,k),sxminus(i,j,k),&

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 498)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:47 PM Modified: /trunk/src/ GRHydro_PPM.F90, GRHydro_PPMM.F90 Log: GRHydro: explicitly cast from CCTK_REAL to int to avoid warning From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.F90 Delta lines: +7 -7 =================================================================== --- trunk/src/GRHydro_PPM.F90 2013-03-28 01:46:46 UTC (rev 497) +++ trunk/src/GRHydro_PPM.F90 2013-03-28 01:47:00 UTC (rev 498) @@ -494,7 +494,7 @@ end do else !!$ Really implement C&W, page 197; which requires stencil 4. do i = 4, nx - 3 - s=sign(one, -dpress(i)) + s=int(sign(one, -dpress(i))) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) if (abs(1.d0 - flatten) > 0.d0) then trivial_rp(i-1) = .false. @@ -829,7 +829,7 @@ end do else !!$ Really implement C&W, page 197; which requires stencil 4. do i = 4, nx - 3 - s=sign(one, -dpress(i)) + s=int(sign(one, -dpress(i))) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) if (abs(1.d0 - flatten) > 0.d0) then trivial_rp(i-1) = .false. @@ -1196,7 +1196,7 @@ end do else !!$ Really implement C&W, page 197; which requires stencil 4. do i = 4, nx - 3 - s=sign(one, -dpress(i)) + s=int(sign(one, -dpress(i))) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) tempplus(i) = flatten * tempplus(i) + (1.d0 - flatten) * temperature(i) tempminus(i) = flatten * tempminus(i) + (1.d0 - flatten) * temperature(i) @@ -1459,7 +1459,7 @@ end do else !!$ Really implement C&W, page 197; which requires stencil 4. do i = 4, nx - 3 - s=sign(one, -dpress(i)) + s=int(sign(one, -dpress(i))) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) tempplus(i) = flatten * tempplus(i) + (1.d0 - flatten) * temperature(i) tempminus(i) = flatten * tempminus(i) + (1.d0 - flatten) * temperature(i) @@ -1679,7 +1679,7 @@ else !!$ Really implement C&W, page 197; which requires stencil 4. do itracer=1,number_of_tracers do i = 4, nx - 3 - s=sign(one, -dpress(i)) + s=int(sign(one, -dpress(i))) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) tracerplus(i,itracer) = flatten * tracerplus(i,itracer) + & (1.d0 - flatten) * tracer(i,itracer) @@ -1910,7 +1910,7 @@ end do else !!$ Really implement C&W, page 197; which requires stencil 4. do i = 4, nx - 3 - s=sign(one, -dpress(i)) + s=int(sign(one, -dpress(i))) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) Y_e_plus(i) = flatten * Y_e_plus(i) + & (1.d0 - flatten) * Y_e(i) @@ -2040,7 +2040,7 @@ end do else !!$ Really implement C&W, page 197; which requires stencil 4. do i = 4, nx - 3 - s=sign(one, -dpress(i)) + s=int(sign(one, -dpress(i))) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) Y_e_plus(i) = flatten * Y_e_plus(i) + & (1.d0 - flatten) * Y_e(i) File [modified]: GRHydro_PPMM.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_PPMM.F90 2013-03-28 01:46:46 UTC (rev 497) +++ trunk/src/GRHydro_PPMM.F90 2013-03-28 01:47:00 UTC (rev 498) @@ -411,7 +411,7 @@ end do else !!$ Really implement C&W, page 197; which requires stencil 4. do i = 4, nx - 3 - s=sign(1.d0, -dpress(i)) + s=int(sign(1.d0, -dpress(i))) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) if (abs(1.d0 - flatten) > 0.d0) then trivial_rp(i-1) = .false.

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 497)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:46 PM Modified: /trunk/src/ GRHydro_HLLE.F90, GRHydro_HLLEM.F90, GRHydro_HLLE_AM.F90 Log: GRHydro: change CCTK_WARN(0,...) to CCTK_ERROR this should tell the compiler that the function never returns and help optimizing the code by indicating eg slow branches of if statements From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_HLLE.F90 Delta lines: +10 -10 =================================================================== --- trunk/src/GRHydro_HLLE.F90 2013-03-28 01:46:42 UTC (rev 496) +++ trunk/src/GRHydro_HLLE.F90 2013-03-28 01:46:46 UTC (rev 497) @@ -216,7 +216,7 @@ call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", & &"trivial") else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if if(evolve_temper.eq.1.and.reconstruct_temper.eq.1) then @@ -272,7 +272,7 @@ avg_beta = 0.5d0 * (beta3(i+xoffset,j+yoffset,k+zoffset) + & beta3(i,j,k)) else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) @@ -311,7 +311,7 @@ velzplus(i,j,k),pressplus(i,j,k),& avg_det,avg_alp,avg_beta) else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if else !!! The end of this branch is right at the bottom of the routine @@ -327,7 +327,7 @@ else if (flux_direction == 3) then usendh = uzzh else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if !!$ Calculate the jumps in the conserved variables @@ -498,7 +498,7 @@ pressminus(i+xoffset,j+yoffset,k+zoffset),& avg_det,avg_alp,avg_beta) else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if !!$ Compute H viscosity if requested @@ -536,7 +536,7 @@ etaY(i,j-1,k, vup, eos_c),& etaY(i,j-1,k+1, vup, eos_c)) else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if ! modify eigenvalues of Roe's matrix by computed H viscosity @@ -678,7 +678,7 @@ call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", & &"trivial") else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil @@ -714,7 +714,7 @@ avg_beta = 0.5d0 * (beta3(i+xoffset,j+yoffset,k+zoffset) + & beta3(i,j,k)) else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) @@ -740,7 +740,7 @@ else if (flux_direction == 3) then usendh = uzzh else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if !!$ Calculate the jumps in the conserved variables @@ -813,7 +813,7 @@ fminus(:) = (velzminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * & cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if !!$ Find minimum and maximum wavespeeds File [modified]: GRHydro_HLLEM.F90 Delta lines: +9 -9 =================================================================== --- trunk/src/GRHydro_HLLEM.F90 2013-03-28 01:46:42 UTC (rev 496) +++ trunk/src/GRHydro_HLLEM.F90 2013-03-28 01:46:46 UTC (rev 497) @@ -119,7 +119,7 @@ call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", & &"trivial") else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if ! constraint transport needs to be able to average fluxes in the directions @@ -229,7 +229,7 @@ else if (flux_direction == 3) then avg_beta=avg_betaz else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) @@ -337,7 +337,7 @@ endif else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if else !!! The end of this branch is right at the bottom of the routine @@ -349,7 +349,7 @@ else if (flux_direction == 3) then usendh = uzzh else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if !!$ Calculate the jumps in the conserved variables @@ -541,7 +541,7 @@ endif else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if @@ -779,7 +779,7 @@ call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", & &"trivial") else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil @@ -840,7 +840,7 @@ avg_beta = 0.5d0 * (beta3(i+xoffset,j+yoffset,k+zoffset) + & beta3(i,j,k)) else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) @@ -865,7 +865,7 @@ else if (flux_direction == 3) then usendh = uzzh else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if !!$ b^2 = (1-v^2)B^2+(B dot v)^2 @@ -929,7 +929,7 @@ fminus(:) = (velzminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * & cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if !!$ Find minimum and maximum wavespeeds File [modified]: GRHydro_HLLE_AM.F90 Delta lines: +9 -9 =================================================================== --- trunk/src/GRHydro_HLLE_AM.F90 2013-03-28 01:46:42 UTC (rev 496) +++ trunk/src/GRHydro_HLLE_AM.F90 2013-03-28 01:46:46 UTC (rev 497) @@ -117,7 +117,7 @@ call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", & &"trivial") else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if ! constraint transport needs to be able to average fluxes in the directions @@ -217,7 +217,7 @@ else if (flux_direction == 3) then avg_beta=avg_betaz else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) @@ -316,7 +316,7 @@ endif else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if else !!! The end of this branch is right at the bottom of the routine @@ -328,7 +328,7 @@ else if (flux_direction == 3) then usendh = uzzh else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if !!$ Calculate the jumps in the conserved variables @@ -505,7 +505,7 @@ endif else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if @@ -736,7 +736,7 @@ call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", & &"trivial") else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil @@ -797,7 +797,7 @@ avg_beta = 0.5d0 * (beta3(i+xoffset,j+yoffset,k+zoffset) + & beta3(i,j,k)) else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) @@ -822,7 +822,7 @@ else if (flux_direction == 3) then usendh = uzzh else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if !!$ b^2 = (1-v^2)B^2+(B dot v)^2 @@ -886,7 +886,7 @@ fminus(:) = (velzminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * & cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) else - call CCTK_WARN(0, "Flux direction not x,y,z") + call CCTK_ERROR("Flux direction not x,y,z") end if !!$ Find minimum and maximum wavespeeds

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 496)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:46 PM Modified: /trunk/src/ GRHydro_EigenproblemM.F90, GRHydro_HLLEM.F90, GRHydro_HLLE_AM.F90, GRHydro_PPMM.F90 Log: GRHydro: pass pressure into eigenvalue routine for MHD bugfix for hot eos case, just more efficient for cold eos From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_EigenproblemM.F90 Delta lines: +16 -21 =================================================================== --- trunk/src/GRHydro_EigenproblemM.F90 2013-03-28 01:46:38 UTC (rev 495) +++ trunk/src/GRHydro_EigenproblemM.F90 2013-03-28 01:46:42 UTC (rev 496) @@ -41,22 +41,22 @@ CONTAINS -subroutine eigenvaluesM(handle,rho,velx,vely,velz,eps,w_lorentz,& +subroutine eigenvaluesM(handle,rho,velx,vely,velz,eps,press,w_lorentz,& Bvcx,Bvcy,Bvcz,lam,gxx,gxy,gxz,gyy,gyz,gzz,u,alp,beta) implicit none DECLARE_CCTK_PARAMETERS - CCTK_REAL rho,velx,vely,velz,eps,w_lorentz - CCTK_REAL Bvcx,Bvcy,Bvcz - CCTK_REAL lam(5) - CCTK_REAL gxx,gxy,gxz,gyy,gyz,gzz - CCTK_REAL alp,beta,u + CCTK_REAL, intent(in) :: rho,velx,vely,velz,eps,w_lorentz + CCTK_REAL, intent(in) :: Bvcx,Bvcy,Bvcz + CCTK_REAL, intent(out) :: lam(5) + CCTK_REAL, intent(in) :: gxx,gxy,gxz,gyy,gyz,gzz + CCTK_REAL, intent(in) :: alp,beta,u CCTK_REAL cs2,one,two,U2 CCTK_REAL vlowx,vlowy,vlowz,v2,w CCTK_REAL lam1,lam2,lam3,lamm,lamp,lamm_nobeta,lamp_nobeta - CCTK_INT handle + CCTK_INT, intent(in) :: handle CCTK_REAL dpdrho,dpdeps,press CCTK_REAL Bvcxlow,Bvcylow,Bvczlow,Bvc2,rhohstar,va2 @@ -75,11 +75,6 @@ !!$ Set required fluid quantities -! call EOS_Omni_cs2(handle,keytemp,GRHydro_eos_rf_prec,n,& -! rho,eps,xtemp,xye,cs2,keyerr,anyerr) - call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& - rho,eps,xtemp,xye,press,keyerr,anyerr) - call EOS_Omni_DPressByDEps(handle,keytemp,GRHydro_eos_rf_prec,n,& rho,eps,xtemp,xye,dpdeps,keyerr,anyerr) @@ -138,24 +133,24 @@ end subroutine eigenvaluesM -subroutine eigenvaluesM_hot(handle,ii,jj,kk,rho,velx,vely,velz,eps,w_lorentz,& +subroutine eigenvaluesM_hot(handle,ii,jj,kk,rho,velx,vely,velz,eps,press,w_lorentz,& Bvcx,Bvcy,Bvcz,temp,ye,lam,gxx,gxy,gxz,gyy,gyz,gzz,u,alp,beta) implicit none DECLARE_CCTK_PARAMETERS - CCTK_REAL rho,velx,vely,velz,eps,w_lorentz - CCTK_REAL Bvcx,Bvcy,Bvcz - CCTK_REAL lam(5) - CCTK_REAL gxx,gxy,gxz,gyy,gyz,gzz - CCTK_REAL alp,beta,u - CCTK_REAL temp,ye + CCTK_REAL, intent(in) :: rho,velx,vely,velz,eps,w_lorentz + CCTK_REAL, intent(in) :: Bvcx,Bvcy,Bvcz + CCTK_REAL, intent(out) :: lam(5) + CCTK_REAL, intent(in) :: gxx,gxy,gxz,gyy,gyz,gzz + CCTK_REAL, intent(in) :: alp,beta,u + CCTK_REAL, intent(in) :: temp,ye + CCTK_REAL, intent(in) :: press CCTK_REAL cs2,one,two,U2 CCTK_REAL vlowx,vlowy,vlowz,v2,w CCTK_REAL lam1,lam2,lam3,lamm,lamp,lamm_nobeta,lamp_nobeta - CCTK_INT handle, ii,jj,kk - CCTK_REAL press + CCTK_INT, intent(in) :: handle, ii,jj,kk CCTK_REAL Bvcxlow,Bvcylow,Bvczlow,Bvc2,rhohstar,va2 CCTK_REAL Bdotv,b2 File [modified]: GRHydro_HLLEM.F90 Delta lines: +28 -28 =================================================================== --- trunk/src/GRHydro_HLLEM.F90 2013-03-28 01:46:38 UTC (rev 495) +++ trunk/src/GRHydro_HLLEM.F90 2013-03-28 01:46:42 UTC (rev 496) @@ -375,26 +375,26 @@ if (flux_direction == 1) then if(evolve_temper.ne.1) then call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & + prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(6),prim_m(7), & prim_m(8),prim_m(9),prim_m(10),& lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& - usendh,avg_alp,avg_beta) + usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & - prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & + prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(6),prim_p(7), & prim_p(8),prim_p(9),prim_p(10),& lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) else xtemp = temperature(i,j,k) call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& - prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & + prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(6),prim_m(7), & prim_m(8),prim_m(9),prim_m(10),& xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) xtemp = temperature(i,j,k) call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& - prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & + prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(6),prim_p(7), & prim_p(8),prim_p(9),prim_p(10),& xtemp,y_e_plus(i,j,k),& lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& @@ -431,26 +431,26 @@ else if (flux_direction == 2) then if(evolve_temper.ne.1) then call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & + prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(6),prim_m(7), & prim_m(9),prim_m(10),prim_m(8),& lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& - usendh,avg_alp,avg_beta) + usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & - prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & + prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(6),prim_p(7), & prim_p(9),prim_p(10),prim_p(8),& lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) else xtemp = temperature(i,j,k) call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& - prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & + prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(6),prim_m(7), & prim_m(9),prim_m(10),prim_m(8),& xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& - usendh,avg_alp,avg_beta) - xtemp = temperature(i,j,k) + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& - prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & + prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(6),prim_p(7), & prim_p(9),prim_p(10),prim_p(8),& xtemp,y_e_plus(i,j,k),& lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& @@ -487,26 +487,26 @@ else if (flux_direction == 3) then if(evolve_temper.ne.1) then call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & + prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(6),prim_m(7), & prim_m(10),prim_m(8),prim_m(9),& lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& - usendh,avg_alp,avg_beta) + usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & - prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & + prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(6),prim_p(7), & prim_p(10),prim_p(8),prim_p(9),& lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) else xtemp = temperature(i,j,k) call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& - prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & + prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(6),prim_m(7), & prim_m(10),prim_m(8),prim_m(9),& xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& - usendh,avg_alp,avg_beta) - xtemp = temperature(i,j,k) + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& - prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & + prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(6),prim_p(7), & prim_p(10),prim_p(8),prim_p(9),& xtemp,y_e_plus(i,j,k),& lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& @@ -885,12 +885,12 @@ if (flux_direction == 1) then call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & + prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(6),prim_m(7), & mag_m(1),mag_m(2),mag_m(3),& lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& - usendh,avg_alp,avg_beta) + usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & - prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & + prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(6),prim_p(7), & mag_p(1),mag_p(2),mag_p(3),& lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) @@ -900,12 +900,12 @@ cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) else if (flux_direction == 2) then call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & + prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(6),prim_m(7), & mag_m(2),mag_m(3),mag_m(1),& lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& - usendh,avg_alp,avg_beta) + usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & - prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & + prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(6),prim_p(7), & mag_p(2),mag_p(3),mag_p(1),& lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) @@ -915,12 +915,12 @@ cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) else if (flux_direction == 3) then call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & + prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(6),prim_m(7), & mag_m(3),mag_m(1),mag_m(2),& lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& - usendh,avg_alp,avg_beta) + usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle,& - prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & + prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(6),prim_p(7), & mag_p(3),mag_p(1),mag_p(2),& lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) File [modified]: GRHydro_HLLE_AM.F90 Delta lines: +28 -28 =================================================================== --- trunk/src/GRHydro_HLLE_AM.F90 2013-03-28 01:46:38 UTC (rev 495) +++ trunk/src/GRHydro_HLLE_AM.F90 2013-03-28 01:46:42 UTC (rev 496) @@ -351,26 +351,26 @@ if (flux_direction == 1) then if(evolve_temper.ne.1) then call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & + prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(6),prim_m(7), & prim_m(8),prim_m(9),prim_m(10),& lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& - usendh,avg_alp,avg_beta) + usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & - prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & + prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(6),prim_p(7), & prim_p(8),prim_p(9),prim_p(10),& lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) else xtemp = temperature(i,j,k) call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& - prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & + prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(6),prim_m(7), & prim_m(8),prim_m(9),prim_m(10),& xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) xtemp = temperature(i,j,k) call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& - prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & + prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(6),prim_p(7), & prim_p(8),prim_p(9),prim_p(10),& xtemp,y_e_plus(i,j,k),& lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& @@ -403,26 +403,26 @@ else if (flux_direction == 2) then if(evolve_temper.ne.1) then call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & + prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(6),prim_m(7), & prim_m(9),prim_m(10),prim_m(8),& lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& - usendh,avg_alp,avg_beta) + usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & - prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & + prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(6),prim_p(7), & prim_p(9),prim_p(10),prim_p(8),& lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) else xtemp = temperature(i,j,k) call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& - prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & + prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(6),prim_m(7), & prim_m(9),prim_m(10),prim_m(8),& xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& - usendh,avg_alp,avg_beta) - xtemp = temperature(i,j,k) + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& - prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & + prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(6),prim_p(7), & prim_p(9),prim_p(10),prim_p(8),& xtemp,y_e_plus(i,j,k),& lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& @@ -455,26 +455,26 @@ else if (flux_direction == 3) then if(evolve_temper.ne.1) then call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & + prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(6),prim_m(7), & prim_m(10),prim_m(8),prim_m(9),& lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& - usendh,avg_alp,avg_beta) + usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & - prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & + prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(6),prim_p(7), & prim_p(10),prim_p(8),prim_p(9),& lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) else xtemp = temperature(i,j,k) call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& - prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & + prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(6),prim_m(7), & prim_m(10),prim_m(8),prim_m(9),& xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& - usendh,avg_alp,avg_beta) - xtemp = temperature(i,j,k) + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& - prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & + prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(6),prim_p(7), & prim_p(10),prim_p(8),prim_p(9),& xtemp,y_e_plus(i,j,k),& lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& @@ -842,12 +842,12 @@ if (flux_direction == 1) then call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & + prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(6),prim_m(7), & mag_m(1),mag_m(2),mag_m(3),& lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& - usendh,avg_alp,avg_beta) + usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & - prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & + prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(6),prim_p(7), & mag_p(1),mag_p(2),mag_p(3),& lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) @@ -857,12 +857,12 @@ cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) else if (flux_direction == 2) then call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & + prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(6),prim_m(7), & mag_m(2),mag_m(3),mag_m(1),& lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& - usendh,avg_alp,avg_beta) + usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & - prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & + prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(6),prim_p(7), & mag_p(2),mag_p(3),mag_p(1),& lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) @@ -872,12 +872,12 @@ cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) else if (flux_direction == 3) then call eigenvaluesM(GRHydro_eos_handle,& - prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & + prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(6),prim_m(7), & mag_m(3),mag_m(1),mag_m(2),& lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& - usendh,avg_alp,avg_beta) + usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle,& - prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & + prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(6),prim_p(7), & mag_p(3),mag_p(1),mag_p(2),& lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) File [modified]: GRHydro_PPMM.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_PPMM.F90 2013-03-28 01:46:38 UTC (rev 495) +++ trunk/src/GRHydro_PPMM.F90 2013-03-28 01:46:42 UTC (rev 496) @@ -247,7 +247,7 @@ det, agxx, agxy, agxz, agyy, agyz, agzz) call eigenvaluesM(handle,& D_UPW(rho), D_UPW(velx), D_UPW(vely), D_UPW(velz), & - D_UPW(eps), D_UPW(w_lorentz),& + D_UPW(eps), D_UPW(press), D_UPW(w_lorentz),& D_UPW(Bvcx), D_UPW(Bvcy), D_UPW(Bvcz), lam, & agxx, agxy, agxz, agyy, agyz, agzz, & uxx, D_UPW(alp), D_UPW(beta))

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 495)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:46 PM Modified: /trunk/src/ GRHydro_Boundaries.F90, GRHydro_Con2Prim.F90, GRHydro_Differences.F90, GRHydro_ENOReconstruct.F90, GRHydro_Eigenproblem.F90, GRHydro_EigenproblemM.F90, GRHydro_Eigenproblem_Marquina.F90, GRHydro_EoSChangeGamma.F90, GRHydro_Flux.F90, GRHydro_FluxSplit.F90, GRHydro_HLLE.F90, GRHydro_HLLEM.F90, GRHydro_HLLE_AM.F90, GRHydro_MP5Reconstruct.F90, GRHydro_Marquina.F90, GRHydro_PPM.F90, GRHydro_PPMM.F90, GRHydro_PPMMReconstruct_drv.F90, GRHydro_PPMReconstruct_drv.F90, GRHydro_Particle.F90, GRHydro_Prim2Con.F90, GRHydro_Prim2ConAM.F90, GRHydro_Prim2ConM.F90, GRHydro_Reconstruct.F90, GRHydro_RiemannSolve.F90, GRHydro_RiemannSolveAM.F90, GRHydro_RiemannSolveM.F90, GRHydro_RoeSolver.F90, GRHydro_WENOReconstruct.F90 Log: GRHydro: remove unused variables From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Boundaries.F90 Delta lines: +0 -3 =================================================================== --- trunk/src/GRHydro_Boundaries.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_Boundaries.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -376,7 +376,6 @@ integer, dimension(3) :: sw integer :: ierr = 0 - integer :: i,j,k CCTK_INT :: GRHydro_UseGeneralCoordinates, general_coordinates CCTK_INT, parameter :: faces=CCTK_ALL_FACES @@ -515,7 +514,6 @@ integer, dimension(3) :: sw integer :: ierr = 0 - integer :: i,j,k CCTK_INT :: GRHydro_UseGeneralCoordinates, general_coordinates CCTK_INT, parameter :: faces=CCTK_ALL_FACES @@ -651,7 +649,6 @@ integer, dimension(3) :: sw integer :: ierr = 0 - integer :: i,j,k CCTK_INT, parameter :: faces=CCTK_ALL_FACES CCTK_INT, parameter :: ione=1 File [modified]: GRHydro_Con2Prim.F90 Delta lines: +8 -9 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_Con2Prim.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -54,7 +54,6 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - integer :: CCTK_MyProc integer :: i, j, k, itracer, nx, ny, nz CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, pmin, epsmin, dummy1, dummy2 logical :: epsnegative @@ -416,11 +415,11 @@ press, uxx, uxy, uxz, uyy, uyz, uzz, det, w_lorentz, x, & y, z, r, GRHydro_rho_min - CCTK_REAL s2, c0, c1, c2, c3, c4, f, df, ftol, v2, w, vlowx, vlowy, vlowz + CCTK_REAL s2, f, df, vlowx, vlowy, vlowz CCTK_INT count, i, handle, GRHydro_reflevel CCTK_REAL GRHydro_C2P_failed, dummy1, dummy2 - CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, epsold, epsnew, w2, & - w2mhalf, temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin + CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, & + temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin character(len=200) warnline logical epsnegative, mustbisect @@ -778,11 +777,11 @@ press, uxx, uxy, uxz, uyy, uyz, uzz, det, w_lorentz, x, & y, z, r, GRHydro_rho_min CCTK_REAL temp, ye - CCTK_REAL s2, c0, c1, c2, c3, c4, f, df, ftol, v2, w, vlowx, vlowy, vlowz + CCTK_REAL s2, f, df, vlowx, vlowy, vlowz CCTK_INT cctk_iteration, ii,jj,kk,count, i, handle, GRHydro_reflevel CCTK_REAL GRHydro_C2P_failed - CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, epsold, epsnew, w2, & - w2mhalf, temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin + CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, & + temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin CCTK_REAL epsminl,pminl,plow,tmp, dummy1, dummy2 CCTK_REAL local_perc_ptol @@ -1248,7 +1247,7 @@ w_lorentz = 1.d0 end if -51 press = pnew + press = pnew vlowx = usx / ( (rho + rho*epsilon + press) * w_lorentz**2) vlowy = usy / ( (rho + rho*epsilon + press) * w_lorentz**2) vlowz = usz / ( (rho + rho*epsilon + press) * w_lorentz**2) @@ -1666,7 +1665,7 @@ CCTK_REAL s2, f, df, vlowx, vlowy, vlowz CCTK_INT count, handle, GRHydro_reflevel - CCTK_REAL udens, usx, usy, usz, rhoold, rhonew, epsold, epsnew, & + CCTK_REAL udens, usx, usy, usz, rhoold, rhonew, & enthalpy, denthalpy, sqrtdet, invsqrtdet, invfac, GRHydro_C2P_failed, dummy1, dummy2 character(len=200) warnline File [modified]: GRHydro_Differences.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_Differences.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_Differences.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -27,7 +27,7 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k, nx, ny, nz, kk, ll + CCTK_INT :: i, j, k, nx, ny, nz CCTK_REAL :: idx, idy, idz CCTK_REAL :: dr_eps = 1.d-2 File [modified]: GRHydro_ENOReconstruct.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_ENOReconstruct.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_ENOReconstruct.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -40,7 +40,7 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - CCTK_INT :: i, j, k, l, m, q, allocstat + CCTK_INT :: i, j, l, m, q, allocstat CCTK_REAL :: denominator, numerator, numerator_product if(.not.coeffs_allocated) then File [modified]: GRHydro_Eigenproblem.F90 Delta lines: +1 -2 =================================================================== --- trunk/src/GRHydro_Eigenproblem.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_Eigenproblem.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -145,7 +145,6 @@ CCTK_REAL vlowx,vlowy,vlowz,v2,w CCTK_REAL lam1,lam2,lam3,lamm,lamp,lamm_nobeta,lamp_nobeta CCTK_INT handle,ii,jj,kk - CCTK_REAL dpdrho,dpdeps,press ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) @@ -598,7 +597,7 @@ CCTK_REAL lam1,lam2,lam3,lamm,lamp,lamm_nobeta,lamp_nobeta CCTK_REAL cs2,one,two CCTK_REAL vlowx,vlowy,vlowz,v2,w - CCTK_REAL press,dpdrho,dpdeps,enthalpy,kappa + CCTK_REAL press,dpdeps,enthalpy,kappa CCTK_REAL axp,axm,vxp,vxm,cxx,cxy,cxz,gam,xsi,dlt CCTK_INT handle File [modified]: GRHydro_EigenproblemM.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_EigenproblemM.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_EigenproblemM.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -155,7 +155,7 @@ CCTK_REAL vlowx,vlowy,vlowz,v2,w CCTK_REAL lam1,lam2,lam3,lamm,lamp,lamm_nobeta,lamp_nobeta CCTK_INT handle, ii,jj,kk - CCTK_REAL dpdrho,dpdeps,press + CCTK_REAL press CCTK_REAL Bvcxlow,Bvcylow,Bvczlow,Bvc2,rhohstar,va2 CCTK_REAL Bdotv,b2 File [modified]: GRHydro_Eigenproblem_Marquina.F90 Delta lines: +10 -14 =================================================================== --- trunk/src/GRHydro_Eigenproblem_Marquina.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_Eigenproblem_Marquina.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -49,8 +49,8 @@ DECLARE_CCTK_PARAMETERS - CCTK_REAL rhor,velxr,velyr,velzr,epsr,w_lorentzr - CCTK_REAL rhol,velxl,velyl,velzl,epsl,w_lorentzl + CCTK_REAL rhor,velxr,velyr,velzr,epsr + CCTK_REAL rhol,velxl,velyl,velzl,epsl CCTK_REAL densl,sxl,syl,szl,taul CCTK_REAL densr,sxr,syr,szr,taur @@ -63,8 +63,8 @@ !!$ LOCAL VARS - CCTK_REAL du(5),aa(5,5),qdiff1,qdiff2,qdiff3,qdiff4,qdiff5 - integer i,j,k,l + CCTK_REAL du(5),aa(5,5) + integer i,j,k CCTK_REAL paug(5,6),tmp1,tmp2,sump,summ,f_du(5) CCTK_REAL leivec1l(5),leivec2l(5),leivec3l(5),leivecpl(5),leivecml(5) CCTK_REAL leivec1r(5),leivec2r(5),leivec3r(5),leivecpr(5),leivecmr(5) @@ -78,8 +78,6 @@ CCTK_REAL vlowxr,vlowyr,vlowzr,v2r,wr CCTK_REAL vlowxl,vlowyl,vlowzl,v2l,wl - CCTK_REAL lamp_nobeta,lamm_nobeta - CCTK_REAL pressl,dpdrhol,dpdepsl,enthalpyl,kappal CCTK_REAL pressr,dpdrhor,dpdepsr,enthalpyr,kappar CCTK_REAL axpl,axml,vxpl,vxml,xsil,dltl @@ -774,8 +772,8 @@ DECLARE_CCTK_PARAMETERS - CCTK_REAL rhor,velxr,velyr,velzr,epsr,w_lorentzr - CCTK_REAL rhol,velxl,velyl,velzl,epsl,w_lorentzl + CCTK_REAL rhor,velxr,velyr,velzr,epsr + CCTK_REAL rhol,velxl,velyl,velzl,epsl CCTK_REAL densl,sxl,syl,szl,taul CCTK_REAL densr,sxr,syr,szr,taur CCTK_REAL templ,tempr,ye_l,ye_r @@ -789,8 +787,8 @@ !!$ LOCAL VARS - CCTK_REAL du(5),aa(5,5),qdiff1,qdiff2,qdiff3,qdiff4,qdiff5 - integer i,j,k,l + CCTK_REAL du(5),aa(5,5) + integer i,j,k CCTK_REAL paug(5,6),tmp1,tmp2,sump,summ,f_du(5) CCTK_REAL leivec1l(5),leivec2l(5),leivec3l(5),leivecpl(5),leivecml(5) CCTK_REAL leivec1r(5),leivec2r(5),leivec3r(5),leivecpr(5),leivecmr(5) @@ -804,10 +802,8 @@ CCTK_REAL vlowxr,vlowyr,vlowzr,v2r,wr CCTK_REAL vlowxl,vlowyl,vlowzl,v2l,wl - CCTK_REAL lamp_nobeta,lamm_nobeta - - CCTK_REAL pressl,dpdrhol,dpdepsl,enthalpyl,kappal - CCTK_REAL pressr,dpdrhor,dpdepsr,enthalpyr,kappar + CCTK_REAL pressl,dpdepsl,enthalpyl,kappal + CCTK_REAL pressr,dpdepsr,enthalpyr,kappar CCTK_REAL axpl,axml,vxpl,vxml,xsil,dltl CCTK_REAL axpr,axmr,vxpr,vxmr,xsir,dltr CCTK_REAL cxx,cxy,cxz,gam,vxa,vxb File [modified]: GRHydro_EoSChangeGamma.F90 Delta lines: +0 -4 =================================================================== --- trunk/src/GRHydro_EoSChangeGamma.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_EoSChangeGamma.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -291,10 +291,6 @@ CCTK_REAL :: local_Gamma, local_k, eos_k_initial - CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: Q - - character(len=100) infoline - ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 File [modified]: GRHydro_Flux.F90 Delta lines: +0 -1 =================================================================== --- trunk/src/GRHydro_Flux.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_Flux.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -36,7 +36,6 @@ CCTK_REAL :: dens,sx,sy,sz,tau,velx,vely,velz,det CCTK_REAL :: densf,sxf,syf,szf,tauf CCTK_REAL :: alp,beta,press - integer :: i,j,k integer :: xoffset,yoffset,zoffset densf = dens * ( (velx*xoffset + vely*yoffset + & File [modified]: GRHydro_FluxSplit.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/GRHydro_FluxSplit.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_FluxSplit.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -40,7 +40,7 @@ DECLARE_CCTK_FUNCTIONS integer :: nx, ny, nz, i, j, k, ierr, max_handle - CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, psi4pt, beta + CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, beta CCTK_REAL, dimension(5) :: lambda CCTK_REAL :: alpha1_local, alpha2_local, alpha3_local, alpha4_local, & alpha5_local @@ -287,7 +287,7 @@ integer :: nx, ny, nz, i, j, k CCTK_REAL, dimension(:), allocatable :: upper, det, dummy - CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, psi4pt + CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz nx = cctk_lsh(1) ny = cctk_lsh(2) File [modified]: GRHydro_HLLE.F90 Delta lines: +4 -4 =================================================================== --- trunk/src/GRHydro_HLLE.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_HLLE.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -166,9 +166,9 @@ CCTK_REAL, dimension(5) :: qdiff CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det, etabar CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & - uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh - CCTK_INT :: type_bits, trivial, not_trivial + CCTK_INT :: type_bits, trivial ! sign requires its arguments to be of identical KIND CCTK_REAL, parameter :: one = 1d0 @@ -634,9 +634,9 @@ CCTK_REAL, dimension(number_of_tracers) :: qdiff CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & - uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh - CCTK_INT :: type_bits, trivial, not_trivial + CCTK_INT :: type_bits, trivial ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates File [modified]: GRHydro_HLLEM.F90 Delta lines: +4 -6 =================================================================== --- trunk/src/GRHydro_HLLEM.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_HLLEM.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -53,8 +53,7 @@ CCTK_REAL, dimension(5) :: lamminus,lamplus CCTK_REAL :: charmin, charmax, charpm,chartop,avg_alp,avg_det, sdet CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & - uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & - cs2_p, cs2_m, dpdeps_p, dpdeps_m + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh CCTK_REAL :: rhoenth_p, rhoenth_m, avg_betax, avg_betay, avg_betaz CCTK_REAL :: vxtp,vytp,vztp,vxtm,vytm,vztm,ab0p,ab0m,b2p,b2m,bdotvp,bdotvm CCTK_REAL :: wp,wm,v2p,v2m,bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,vA2m,vA2p @@ -66,7 +65,7 @@ CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm CCTK_REAL :: charmax_dc, charmin_dc, charpm_dc - CCTK_INT :: type_bits, trivial, not_trivial + CCTK_INT :: type_bits, trivial CCTK_REAL :: xtemp ! save memory when MP is not used @@ -724,11 +723,10 @@ CCTK_REAL, dimension(3) :: mag_p, mag_m CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & - uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & - cs2_p, cs2_m, dpdeps_p, dpdeps_m + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh CCTK_REAL :: b2p,b2m,vA2m,vA2p - CCTK_INT :: type_bits, trivial, not_trivial + CCTK_INT :: type_bits, trivial ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates File [modified]: GRHydro_HLLE_AM.F90 Delta lines: +4 -6 =================================================================== --- trunk/src/GRHydro_HLLE_AM.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_HLLE_AM.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -53,8 +53,7 @@ CCTK_REAL, dimension(5) :: lamminus,lamplus CCTK_REAL :: charmin, charmax, charpm,chartop,avg_alp,avg_det, sdet CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & - uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & - cs2_p, cs2_m, dpdeps_p, dpdeps_m + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh CCTK_REAL :: rhoenth_p, rhoenth_m, avg_betax, avg_betay, avg_betaz CCTK_REAL :: vxtp,vytp,vztp,vxtm,vytm,vztm,ab0p,ab0m,b2p,b2m,bdotvp,bdotvm CCTK_REAL :: wp,wm,v2p,v2m,bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,vA2m,vA2p @@ -64,7 +63,7 @@ CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm CCTK_REAL :: charmax_dc, charmin_dc, charpm_dc - CCTK_INT :: type_bits, trivial, not_trivial + CCTK_INT :: type_bits, trivial CCTK_REAL :: xtemp ! save memory when MP is not used @@ -681,11 +680,10 @@ CCTK_REAL, dimension(3) :: mag_p, mag_m CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & - uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & - cs2_p, cs2_m, dpdeps_p, dpdeps_m + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh CCTK_REAL :: b2p,b2m,vA2m,vA2p - CCTK_INT :: type_bits, trivial, not_trivial + CCTK_INT :: type_bits, trivial ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates File [modified]: GRHydro_MP5Reconstruct.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/GRHydro_MP5Reconstruct.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_MP5Reconstruct.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -40,7 +40,7 @@ DECLARE_CCTK_PARAMETERS - CCTK_INT :: nx, i, j, k, r + CCTK_INT :: nx, i CCTK_REAL, dimension(nx) :: v, vplus, vminus CCTK_INT, dimension(nx) :: hydro_excision_mask @@ -95,7 +95,7 @@ #define MINMOD4(w,x,y,z) \ 0.125d0*( sign(one,w)+sign(one,x) )*abs( (sign(one,w)+sign(one,y)) * (sign(one,w)+sign(one,z)) )*min(abs(w), abs(x), abs(y), abs(z)) -#define MP5(am2, am1, a, ap1, ap2, arecon) \ +#define MP5(am2, am1, a, ap1, ap2, arecon) \ vl = (2.0d0*am2 - 13.0d0*am1 + 47.0d0*a + 27.0d0*ap1 - 3.0d0*ap2)/60.0d0 &&\ vmp = a + MINMOD( ap1-a, mp5_alpha*(a-am1) ) &&\ if ((vl-a)*(vl-vmp) .le. mp5_eps) then &&\ File [modified]: GRHydro_Marquina.F90 Delta lines: +1 -2 =================================================================== --- trunk/src/GRHydro_Marquina.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_Marquina.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -45,13 +45,12 @@ consp,consm_i,fplus,fminus,f_marquina,primp,primm_i CCTK_REAL :: avg_alp,avg_beta,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & avg_det,uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,& - pressp,pressm_i, & tmp_w_lorentzp, tmp_w_lorentzm_i, w_lorentzp,w_lorentzm_i, usendh, psi4h integer :: m integer :: i,j,k integer :: keytemp - CCTK_INT :: type_bits, trivial, not_trivial + CCTK_INT :: type_bits, trivial if(evolve_temper.eq.1.and.reconstruct_temper.eq.1) then keytemp = 1 File [modified]: GRHydro_PPM.F90 Delta lines: +2 -18 =================================================================== --- trunk/src/GRHydro_PPM.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_PPM.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -73,13 +73,6 @@ CCTK_REAL, dimension(nx) :: rho,velx,vely,velz,eps CCTK_REAL, dimension(nx) :: rhominus,velxminus,velyminus,velzminus,epsminus CCTK_REAL, dimension(nx) :: rhoplus,velxplus,velyplus,velzplus,epsplus - CCTK_REAL, dimension(nx) :: rhominusl,velxminusl,velyminusl,velzminusl - CCTK_REAL, dimension(nx) :: epsminusl - CCTK_REAL, dimension(nx) :: rhoplusl,velxplusl,velyplusl,velzplusl,epsplusl - CCTK_REAL, dimension(nx) :: rhominusr,velxminusr,velyminusr,velzminusr - CCTK_REAL, dimension(nx) :: epsminusr - CCTK_REAL, dimension(nx) :: rhoplusr,velxplusr,velyplusr,velzplusr,epsplusr - CCTK_REAL, dimension(nx) :: atmosphere_mask CCTK_INT :: i,s CCTK_REAL, dimension(nx) :: drho,dvelx,dvely,dvelz,deps @@ -98,7 +91,6 @@ CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det CCTK_REAL, dimension(5) :: lam - CCTK_REAL :: dupw, dloc, delta CCTK_REAL :: agxx, agxy, agxz, agyy, agyz, agzz CCTK_REAL, dimension(nx) :: xwind, l_ev_l, l_ev_r @@ -1046,17 +1038,12 @@ CCTK_REAL, dimension(nx) :: temperature, velx CCTK_REAL, dimension(nx) :: tempminus CCTK_REAL, dimension(nx) :: tempplus - CCTK_REAL, dimension(nx) :: tempminusl - CCTK_REAL, dimension(nx) :: tempplusl - CCTK_REAL, dimension(nx) :: tempminusr - CCTK_REAL, dimension(nx) :: tempplusr - CCTK_REAL, dimension(nx) :: atmosphere_mask CCTK_INT :: i,s CCTK_REAL, dimension(nx) :: dtemp CCTK_REAL, dimension(nx) :: dmtemp CCTK_REAL, dimension(nx) :: press,dpress,tilde_flatten - CCTK_REAL :: dpress2,dvel,w,flatten,eta,etatilde + CCTK_REAL :: dpress2,dvel,w,flatten CCTK_INT, dimension(nx) :: hydro_excision_mask @@ -1569,7 +1556,6 @@ CCTK_REAL, dimension(nx) :: press,dpress,tilde_flatten CCTK_REAL, dimension(nx,number_of_tracers) :: dmtracer, dtracer, tracerflat!, d2tracer CCTK_REAL :: dpress2,w,flatten,dvel - CCTK_REAL :: eta, etatilde !!$ Average slopes delta_m(a). See (1.7) of Colella and Woodward, p.178 !!$ This is the expression for an even grid. @@ -1777,14 +1763,12 @@ CCTK_REAL :: dx CCTK_REAL, dimension(nx) :: rho,velx,vely,velz CCTK_REAL, dimension(nx) :: Y_e,Y_e_minus,Y_e_plus - CCTK_REAL :: Y_eflatomega CCTK_INT :: i,s CCTK_REAL, dimension(nx) :: press,dpress,tilde_flatten - CCTK_REAL, dimension(nx) :: dmY_e, dY_e, Y_eflat!, d2tracer + CCTK_REAL, dimension(nx) :: dmY_e, dY_e CCTK_REAL :: dpress2,w,flatten,dvel - CCTK_REAL :: eta, etatilde CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax File [modified]: GRHydro_PPMM.F90 Delta lines: +0 -1 =================================================================== --- trunk/src/GRHydro_PPMM.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_PPMM.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -73,7 +73,6 @@ CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det CCTK_REAL, dimension(5) :: lam - CCTK_REAL :: dupw, dloc, delta CCTK_REAL :: agxx, agxy, agxz, agyy, agyz, agzz CCTK_REAL, dimension(nx) :: xwind, l_ev_l, l_ev_r File [modified]: GRHydro_PPMMReconstruct_drv.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_PPMMReconstruct_drv.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_PPMMReconstruct_drv.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -60,7 +60,7 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - integer :: nx, ny, nz, i, j, k, itracer + integer :: nx, ny, nz, i, j, k logical, dimension(:,:,:), allocatable :: trivial_rp File [modified]: GRHydro_PPMReconstruct_drv.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_PPMReconstruct_drv.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -53,7 +53,7 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - integer :: nx, ny, nz, i, j, k, itracer + integer :: nx, ny, nz, i, j, k logical, dimension(:,:,:), allocatable :: trivial_rp File [modified]: GRHydro_Particle.F90 Delta lines: +1 -2 =================================================================== --- trunk/src/GRHydro_Particle.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_Particle.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -36,14 +36,13 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, l + CCTK_INT :: i CCTK_INT :: interp_handle, table_handle, status, coord_system_handle character(len=200) :: particle_interp character(len=128) :: warn_message CCTK_INT :: particle_interp_len character(len=7) :: particle_order - character(len=15) :: vname CCTK_INT, dimension(1) :: lsh CCTK_POINTER, dimension(3) :: interp_coords File [modified]: GRHydro_Prim2Con.F90 Delta lines: +0 -1 =================================================================== --- trunk/src/GRHydro_Prim2Con.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_Prim2Con.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -441,7 +441,6 @@ CCTK_INT :: i, j, k CCTK_INT :: keytemp CCTK_REAL :: det - character(len=512) :: warnline ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates File [modified]: GRHydro_Prim2ConAM.F90 Delta lines: +1 -2 =================================================================== --- trunk/src/GRHydro_Prim2ConAM.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_Prim2ConAM.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -438,7 +438,6 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - CCTK_REAL :: xtemp(1) CCTK_INT :: i, j, k CCTK_REAL :: det CCTK_REAL :: maxtau0 @@ -631,7 +630,7 @@ CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, & drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, & - w_tmp, w, vlowx, vlowy, vlowz, sqrtdet + w_tmp, w, vlowx, vlowy, vlowz CCTK_INT :: handle CCTK_REAL :: Bdotv,ab0,b2,blowx,blowy,blowz character(len=256) NaN_WarnLine File [modified]: GRHydro_Prim2ConM.F90 Delta lines: +1 -2 =================================================================== --- trunk/src/GRHydro_Prim2ConM.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_Prim2ConM.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -455,7 +455,6 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - CCTK_REAL :: xtemp(1) CCTK_INT :: i, j, k CCTK_REAL :: det CCTK_REAL :: maxtau0 @@ -654,7 +653,7 @@ CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, dBconsx, dBconsy, dBconsz, & drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, & - w_tmp, w, vlowx, vlowy, vlowz, sqrtdet + w_tmp, w, vlowx, vlowy, vlowz CCTK_INT :: handle CCTK_REAL :: Bdotv,ab0,b2,blowx,blowy,blowz character(len=256) NaN_WarnLine File [modified]: GRHydro_Reconstruct.F90 Delta lines: +0 -1 =================================================================== --- trunk/src/GRHydro_Reconstruct.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_Reconstruct.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -49,7 +49,6 @@ CCTK_REAL :: local_min_tracer, dummy1, dummy2 CCTK_INT :: type_bits, not_trivial CCTK_REAL :: agxx, agxy, agxz, agyy, agyz, agzz, w - character*256 :: warnline ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates File [modified]: GRHydro_RiemannSolve.F90 Delta lines: +0 -4 =================================================================== --- trunk/src/GRHydro_RiemannSolve.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_RiemannSolve.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -36,8 +36,6 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i,j,k - if (CCTK_EQUALS(riemann_solver,"HLLE")) then call GRHydro_HLLE(CCTK_PASS_FTOF) @@ -93,8 +91,6 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i,j,k - if (CCTK_EQUALS(riemann_solver,"HLLE")) then call GRHydro_HLLE(CCTK_PASS_FTOF) File [modified]: GRHydro_RiemannSolveAM.F90 Delta lines: +0 -4 =================================================================== --- trunk/src/GRHydro_RiemannSolveAM.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_RiemannSolveAM.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -36,8 +36,6 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i,j,k - if (CCTK_EQUALS(riemann_solver,"HLLE").or.CCTK_EQUALS(riemann_solver,"LLF")) then call GRHydro_HLLE_AM(CCTK_PASS_FTOF) @@ -100,8 +98,6 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i,j,k - if (CCTK_EQUALS(riemann_solver,"HLLE")) then call GRHydro_HLLE_AM(CCTK_PASS_FTOF) File [modified]: GRHydro_RiemannSolveM.F90 Delta lines: +0 -4 =================================================================== --- trunk/src/GRHydro_RiemannSolveM.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_RiemannSolveM.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -36,8 +36,6 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i,j,k - if (CCTK_EQUALS(riemann_solver,"HLLE").or.CCTK_EQUALS(riemann_solver,"LLF")) then call GRHydro_HLLEM(CCTK_PASS_FTOF) @@ -100,8 +98,6 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i,j,k - if (CCTK_EQUALS(riemann_solver,"HLLE")) then call GRHydro_HLLEM(CCTK_PASS_FTOF) File [modified]: GRHydro_RoeSolver.F90 Delta lines: +4 -4 =================================================================== --- trunk/src/GRHydro_RoeSolver.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_RoeSolver.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -38,15 +38,15 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS CCTK_REAL, dimension(5) :: roeflux,roeave,qdiff,consp,consm_i,& - fplus,fminus,f_roe,primp,primm_i,consh + fplus,fminus,f_roe,primp,primm_i CCTK_REAL :: avg_alp,avg_beta,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & avg_det,uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,& - rhoave, velxave, velyave, velzave, pressave, epsave, & - w_lorentzave, usendh, alp_l, alp_r, psi4h + rhoave, velxave, velyave, velzave, epsave, & + w_lorentzave, usendh integer :: m integer :: i,j,k - CCTK_INT :: type_bits, trivial, not_trivial + CCTK_INT :: type_bits, trivial if (flux_direction == 1) then call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX") File [modified]: GRHydro_WENOReconstruct.F90 Delta lines: +3 -4 =================================================================== --- trunk/src/GRHydro_WENOReconstruct.F90 2013-03-28 01:46:34 UTC (rev 494) +++ trunk/src/GRHydro_WENOReconstruct.F90 2013-03-28 01:46:38 UTC (rev 495) @@ -40,8 +40,7 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - CCTK_INT :: i, j, k, l, m, q, allocstat - CCTK_REAL :: denominator, numerator, numerator_product + CCTK_INT :: allocstat if(.not.coeffs_allocated) then ! Right now we hardcode to 5th order @@ -156,7 +155,7 @@ DECLARE_CCTK_PARAMETERS - CCTK_INT :: order, nx, i, j, k, r + CCTK_INT :: order, nx, i, j CCTK_REAL, dimension(nx) :: v, vplus, vminus CCTK_INT, dimension(nx) :: hydro_excision_mask @@ -164,7 +163,7 @@ logical, dimension(nx) :: excise logical :: normal_weno - CCTK_REAL :: large = 1.d10, gamma1, gamma2, gamma3, beta1, beta2, beta3, vnorm, betanorm + CCTK_REAL :: beta1, beta2, beta3, vnorm, betanorm CCTK_REAL :: wplus1, wplus2, wplus3, wbarplus1, wbarplus2, wbarplus3 CCTK_REAL :: wminus1, wminus2, wminus3, wbarminus1, wbarminus2, wbarminus3

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 494)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:46 PM Modified: /trunk/src/ GRHydro_Eigenproblem.F90 Log: GRHydro: reduce size of critical section in eigenbalues routine From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Eigenproblem.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/GRHydro_Eigenproblem.F90 2013-03-28 01:46:32 UTC (rev 493) +++ trunk/src/GRHydro_Eigenproblem.F90 2013-03-28 01:46:34 UTC (rev 494) @@ -86,16 +86,16 @@ (1.0d0 + eps + press/rho) if(cs2.lt.0.0d0) then - !$OMP CRITICAL if (abs(cs2) .gt. 1.0d-4) then + !$OMP CRITICAL write(warnline,'(a60,6g16.7)') 'abs(cs2), rho, dpdrho, press*dpdeps/rho**2, eps, press/rho: ', abs(cs2), rho, dpdrho, press * dpdeps / (rho**2), eps, press/rho call CCTK_WARN(1,warnline) call CCTK_WARN(1,"cs2 < 0! Check speed of sound calculation!") + !$OMP END CRITICAL cs2 = 0.0d0 else cs2 = 0.0d0 endif - !$OMP END CRITICAL endif

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 493)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:46 PM Modified: /trunk/src/ GRHydro_HLLE.F90 Log: GRHydro: convert more tabs to (8) spaces. No other changes. From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_HLLE.F90 Delta lines: +4 -4 =================================================================== --- trunk/src/GRHydro_HLLE.F90 2013-03-28 01:46:29 UTC (rev 492) +++ trunk/src/GRHydro_HLLE.F90 2013-03-28 01:46:32 UTC (rev 493) @@ -114,13 +114,13 @@ if (evolve_temper.ne.1) then call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& - rho(i,j,k),eps(i,j,k),xtemp,xye,xpress,keyerr,anyerr) + rho(i,j,k),eps(i,j,k),xtemp,xye,xpress,keyerr,anyerr) - call EOS_Omni_DPressByDEps(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& - rho(i,j,k),eps(i,j,k),xtemp,xye,dpdeps,keyerr,anyerr) + call EOS_Omni_DPressByDEps(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,dpdeps,keyerr,anyerr) call EOS_Omni_DPressByDRho(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& - rho(i,j,k),eps(i,j,k),xtemp,xye,dpdrho,keyerr,anyerr) + rho(i,j,k),eps(i,j,k),xtemp,xye,dpdrho,keyerr,anyerr) cs2 = (dpdrho + xpress(1) * dpdeps / (rho(i,j,k)**2))/ & (1.0d0 + eps(i,j,k) + xpress(1)/rho(i,j,k))

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits