Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[McLachlan]McLachlan branch, dgfe, updated. 1adad32dceef2f27eda583cfd56f0034aab5ec34
by Carpet git repositories 28 Aug '12

28 Aug '12

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "McLachlan". The branch, dgfe has been updated via 1adad32dceef2f27eda583cfd56f0034aab5ec34 (commit) from f005f555f8b1108310d5a49e1636edb6b6299355 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 1adad32dceef2f27eda583cfd56f0034aab5ec34 Author: Erik Schnetter <schnetter(a)gmail.com> Date: Tue Aug 28 12:08:12 2012 -0400 Remove temporary files from repository ----------------------------------------------------------------------- Summary of changes: m/CL_BSSN_DG_P5/configuration.ccl | 7 - m/CL_BSSN_DG_P5/interface.ccl | 580 ---- m/CL_BSSN_DG_P5/param.ccl | 2320 ------------- m/CL_BSSN_DG_P5/schedule.ccl | 1286 ------- m/CL_BSSN_DG_P5/src/Boundaries.cc | 3629 -------------------- m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_Constraints.cc | 964 ------ .../src/CL_BSSN_DG_P5_ConstraintsBoundary.cc | 395 --- .../src/CL_BSSN_DG_P5_ConvertFromADMBase.cc | 378 -- .../src/CL_BSSN_DG_P5_ConvertFromADMBaseDerivs.cc | 568 --- .../src/CL_BSSN_DG_P5_ConvertToADMBase.cc | 350 -- m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_Enforce.cc | 399 --- m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_EvolFlux.cc | 854 ----- m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_EvolRHS.cc | 2403 ------------- m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_RHSBoundary.cc | 1070 ------ m/CL_BSSN_DG_P5/src/Differencing.h | 144 - m/CL_BSSN_DG_P5/src/RegisterMoL.cc | 75 - m/CL_BSSN_DG_P5/src/RegisterSymmetries.cc | 1369 -------- m/CL_BSSN_DG_P5/src/Startup.cc | 10 - m/CL_BSSN_DG_P5/src/make.code.defn | 3 - m/CL_BSSN_DG_P5_Helper/configuration.ccl | 1 - m/CL_BSSN_DG_P5_Helper/interface.ccl | 38 - m/CL_BSSN_DG_P5_Helper/param.ccl | 28 - m/CL_BSSN_DG_P5_Helper/schedule.ccl | 86 - m/CL_BSSN_DG_P5_Helper/src/RegisterConstrained.c | 32 - m/CL_BSSN_DG_P5_Helper/src/RegisterSlicing.c | 10 - m/CL_BSSN_DG_P5_Helper/src/SelectBCsADMBase.c | 31 - m/CL_BSSN_DG_P5_Helper/src/SetGroupTags.c | 102 - m/CL_BSSN_DG_P5_Helper/src/make.code.defn | 2 - 28 files changed, 0 insertions(+), 17134 deletions(-) delete mode 100644 m/CL_BSSN_DG_P5/configuration.ccl delete mode 100644 m/CL_BSSN_DG_P5/interface.ccl delete mode 100644 m/CL_BSSN_DG_P5/param.ccl delete mode 100644 m/CL_BSSN_DG_P5/schedule.ccl delete mode 100644 m/CL_BSSN_DG_P5/src/Boundaries.cc delete mode 100644 m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_Constraints.cc delete mode 100644 m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_ConstraintsBoundary.cc delete mode 100644 m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_ConvertFromADMBase.cc delete mode 100644 m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_ConvertFromADMBaseDerivs.cc delete mode 100644 m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_ConvertToADMBase.cc delete mode 100644 m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_Enforce.cc delete mode 100644 m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_EvolFlux.cc delete mode 100644 m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_EvolRHS.cc delete mode 100644 m/CL_BSSN_DG_P5/src/CL_BSSN_DG_P5_RHSBoundary.cc delete mode 100644 m/CL_BSSN_DG_P5/src/Differencing.h delete mode 100644 m/CL_BSSN_DG_P5/src/RegisterMoL.cc delete mode 100644 m/CL_BSSN_DG_P5/src/RegisterSymmetries.cc delete mode 100644 m/CL_BSSN_DG_P5/src/Startup.cc delete mode 100644 m/CL_BSSN_DG_P5/src/make.code.defn delete mode 100644 m/CL_BSSN_DG_P5_Helper/configuration.ccl delete mode 100644 m/CL_BSSN_DG_P5_Helper/interface.ccl delete mode 100644 m/CL_BSSN_DG_P5_Helper/param.ccl delete mode 100644 m/CL_BSSN_DG_P5_Helper/schedule.ccl delete mode 100644 m/CL_BSSN_DG_P5_Helper/src/RegisterConstrained.c delete mode 100644 m/CL_BSSN_DG_P5_Helper/src/RegisterSlicing.c delete mode 100644 m/CL_BSSN_DG_P5_Helper/src/SelectBCsADMBase.c delete mode 100644 m/CL_BSSN_DG_P5_Helper/src/SetGroupTags.c delete mode 100644 m/CL_BSSN_DG_P5_Helper/src/make.code.defn hooks/post-receive -- McLachlan

1 0

[McLachlan]McLachlan branch, dgfe, updated. f005f555f8b1108310d5a49e1636edb6b6299355
by Carpet git repositories 28 Aug '12

28 Aug '12

1 0

[McLachlan]McLachlan branch, dgfe, updated. 20ab5c82397b96b5c3b2a439e1622c14a2ff2f6b
by Carpet git repositories 28 Aug '12

28 Aug '12

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 152)
by rhaas＠tapir.caltech.edu 27 Aug '12

27 Aug '12

User: rhaas Date: 2012/08/27 04:19 PM Modified: /trunk/src/ GRHydro_Bondi.c, GRHydro_BondiM.c Log: remove unused GRHydro_reflevel declaration that causes compilation to fail File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Bondi.c Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_Bondi.c 2012-08-09 20:02:21 UTC (rev 151) +++ trunk/src/GRHydro_Bondi.c 2012-08-27 21:19:38 UTC (rev 152) @@ -808,7 +808,7 @@ CCTK_INT *bad_bondi; CCTK_REAL dlogr,logrmin; CCTK_REAL rmin_bondi,rmax_bondi; - CCTK_INT GRHydro_reflevel=0; +// CCTK_INT GRHydro_reflevel=0; CCTK_INT nbondi; CCTK_INT ileft, iright, nans_exist; File [modified]: GRHydro_BondiM.c Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_BondiM.c 2012-08-09 20:02:21 UTC (rev 151) +++ trunk/src/GRHydro_BondiM.c 2012-08-27 21:19:38 UTC (rev 152) @@ -807,7 +807,7 @@ CCTK_REAL *r_bondi, *logr_bondi, *rho_bondi, *u_bondi, *v_bondi; CCTK_REAL dlogr,logrmin; CCTK_REAL rmin_bondi,rmax_bondi; - CCTK_INT GRHydro_reflevel=0; +// CCTK_INT GRHydro_reflevel=0; int nbondi;

1 0

New changesets in Carpet Mercurial
by Carpet mercurial repositories 27 Aug '12

27 Aug '12

changeset: 3633:63f147d03fdd user: Erik Schnetter <schnetter(a)gmail.com> date: Mon Aug 27 16:05:14 2012 -0400 summary: Carpet: Improve layout of grid structure files changeset: 3634:23af8058871c user: Erik Schnetter <schnetter(a)gmail.com> date: Mon Aug 27 16:05:42 2012 -0400 summary: Carpet: Initialise outer_boundaries during recursive load balancing changeset: 3635:63e9e192816a tag: tip user: Erik Schnetter <schnetter(a)gmail.com> date: Mon Aug 27 16:06:00 2012 -0400 summary: CarpetInterp2: White space change

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 419)
by rhaas＠tapir.caltech.edu 27 Aug '12

27 Aug '12

User: rhaas Date: 2012/08/27 02:19 PM Modified: /trunk/src/ GRHydro_UpdateMaskM.F90 Log: GRHydro: clear atmosphere_mask_real in MHD AtmosphereReset From: Roland Haas <roland.haas(a)physics.gatech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_UpdateMaskM.F90 Delta lines: +1 -0 =================================================================== --- trunk/src/GRHydro_UpdateMaskM.F90 2012-08-27 19:19:44 UTC (rev 418) +++ trunk/src/GRHydro_UpdateMaskM.F90 2012-08-27 19:19:46 UTC (rev 419) @@ -153,6 +153,7 @@ Bprim(i,j,k,1),Bprim(i,j,k,2),Bprim(i,j,k,3),w_lorentz(i,j,k)) if (wk_atmosphere .eq. 0) then atmosphere_mask(i, j, k) = 0 + atmosphere_mask_real(i, j, k) = 0 end if end if

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 418)
by rhaas＠tapir.caltech.edu 27 Aug '12

27 Aug '12

User: rhaas Date: 2012/08/27 02:19 PM Modified: /trunk/src/ GRHydro_UpdateMask.F90, GRHydro_UpdateMaskM.F90 Log: GRHydro: add OpenMP parallelization to IntialAtmosphereReset From: Roland Haas <roland.haas(a)physics.gatech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_UpdateMask.F90 Delta lines: +7 -0 =================================================================== --- trunk/src/GRHydro_UpdateMask.F90 2012-08-27 19:19:40 UTC (rev 417) +++ trunk/src/GRHydro_UpdateMask.F90 2012-08-27 19:19:44 UTC (rev 418) @@ -201,6 +201,7 @@ CCTK_INT :: i,j,k + !$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) @@ -211,6 +212,7 @@ end do end do end do + !$OMP END PARALLEL DO call CCTK_INFO("Setting up the atmosphere mask: points with rho<rho_min are set to be atmosphere") @@ -445,6 +447,7 @@ rho_min = rho_min * initial_atmosphere_factor endif +!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) @@ -477,6 +480,7 @@ temperature(i,j,k),y_e(i,j,k)) y_e_con(i,j,k) = dens(i,j,k) * y_e(i,j,k) else + keytemp = 0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& @@ -518,6 +522,7 @@ temperature_p(i,j,k),y_e_p(i,j,k)) y_e_con_p(i,j,k) = dens_p(i,j,k) * y_e_p(i,j,k) else + keytemp = 0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho_p(i,j,k),eps_p(i,j,k),xtemp,xye,press_p(i,j,k),keyerr,anyerr) call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& @@ -559,6 +564,7 @@ temperature_p_p(i,j,k),y_e_p_p(i,j,k)) y_e_con_p_p(i,j,k) = dens_p_p(i,j,k) * y_e_p_p(i,j,k) else + keytemp = 0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho_p_p(i,j,k),eps_p_p(i,j,k),xtemp,xye,press_p_p(i,j,k),keyerr,anyerr) call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& @@ -576,6 +582,7 @@ end do end do end do +!$OMP END PARALLEL DO end subroutine GRHydro_InitialAtmosphereReset File [modified]: GRHydro_UpdateMaskM.F90 Delta lines: +6 -0 =================================================================== --- trunk/src/GRHydro_UpdateMaskM.F90 2012-08-27 19:19:40 UTC (rev 417) +++ trunk/src/GRHydro_UpdateMaskM.F90 2012-08-27 19:19:44 UTC (rev 418) @@ -161,6 +161,7 @@ end do end do end do +!$OMP END PARALLEL DO !!$ call GRHydro_BoundariesM(CCTK_PASS_FTOF) #undef faulty_gxx @@ -316,6 +317,7 @@ rho_min = rho_min * initial_atmosphere_factor endif +!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) @@ -350,6 +352,7 @@ else + keytemp = 0 call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& @@ -396,6 +399,7 @@ else + keytemp = 0 call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& rho_p(i,j,k),eps_p(i,j,k),xtemp,xye,press_p(i,j,k),keyerr,anyerr) call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& @@ -444,6 +448,7 @@ else + keytemp = 0 call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& rho_p_p(i,j,k),eps_p_p(i,j,k),xtemp,xye,press_p_p(i,j,k),keyerr,anyerr) call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& @@ -464,6 +469,7 @@ end do end do end do +!$OMP END PARALLEL DO write(*,*) " GRHydro_InitialAtmosphereReset" !!$ call GRHydro_BoundariesM(CCTK_PASS_FTOF)

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 417)
by rhaas＠tapir.caltech.edu 27 Aug '12

27 Aug '12

User: rhaas Date: 2012/08/27 02:19 PM Removed: /trunk/src/ GRHydro_RegisterVarsM.cc Log: GRHydro: remove unused file From: Roland Haas <roland.haas(a)physics.gatech.edu> File Changes: Directory: /trunk/src/ ====================== File [removed]: GRHydro_RegisterVarsM.cc Delta lines: +0 -129 =================================================================== --- trunk/src/GRHydro_RegisterVarsM.cc 2012-08-27 19:19:38 UTC (rev 416) +++ trunk/src/GRHydro_RegisterVarsM.cc 2012-08-27 19:19:40 UTC (rev 417) @@ -1,129 +0,0 @@ -// GRHydro_RegisterVars.cc -// -// converted from F90 to improve readability and maintainability -// -// Frank Loeffler - -#include <cstdio> -#include <string> - -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" - -using namespace std; - -// Utility functions to register variables with MoL -// Note: We could check for the return value here, but MoL issues a -// level 0 warning in that case anyway. If that changes in the -// future, a check can simply be inserted here. -static void register_evolved(string v1, string v2) -{ - MoLRegisterEvolvedGroup(CCTK_GroupIndex(v1.c_str()), CCTK_GroupIndex(v2.c_str())); -} -static void register_constrained(string v1) -{ - MoLRegisterConstrainedGroup(CCTK_GroupIndex(v1.c_str())); -} -static void register_saveandrestore(string v1) -{ - MoLRegisterSaveAndRestoreGroup(CCTK_GroupIndex(v1.c_str())); -} - -// Main function called by Cactus to register variables with MoL - -extern "C"void GRHydro_RegisterM(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - // We need some aliased functions, so we first check if they are available - string needed_funs[5] = {"MoLRegisterEvolvedGroup", - "MoLRegisterConstrainedGroup", - "MoLRegisterSaveAndRestoreGroup", - "MoLRegisterEvolved", - "MoLRegisterConstrained"}; - for (int i = 0; i < 5; i++) - if (!CCTK_IsFunctionAliased(needed_funs[i].c_str())) - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "The function \"%s\" has not been aliased!", - needed_funs[i].c_str()); - - // Now we can set which variables have to be registered as which type with MoL - register_constrained("HydroBase::rho"); - register_constrained("HydroBase::press"); - register_constrained("HydroBase::eps"); - register_constrained("HydroBase::vel"); - register_constrained("HydroBase::w_lorentz"); - - if (CCTK_EQUALS(evolution_method, "GRHydro")) - { - // dens and scon - register_evolved("GRHydro::dens", "GRHydro::densrhs"); - register_evolved("GRHydro::scon", "GRHydro::srhs"); - register_evolved("HydroBase::Bvec", "GRHydro::Bvecrhs"); - - // tau - if (CCTK_EQUALS(GRHydro_eos_type, "General")) - register_evolved("GRHydro::tau" , "GRHydro::taurhs"); - else if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) - register_constrained("GRHydro::tau"); - else - CCTK_WARN(0, "Don't recognize the type of EOS!"); - - if(clean_divergence) { - register_evolved("GRHydro::psidc" , "GRHydro::psidcrhs"); - } - - // lapse, metric, curv - register_saveandrestore("admbase::lapse"); - register_saveandrestore("admbase::metric"); - register_saveandrestore("admbase::curv"); - - // shift - if (!CCTK_EQUALS(initial_shift, "none")) - { - if (CCTK_EQUALS(shift_evolution_method, "Comoving")) - { - register_constrained("admbase::shift"); - register_evolved("GRHydro::GRHydro_coords", "GRHydro::GRHydro_coords_rhs"); - } - else - register_saveandrestore("admbase::shift"); - } - - // tracer - if (evolve_tracer != 0) - register_evolved("GRHydro::GRHydro_cons_tracers", - "GRHydro::GRHydro_tracer_rhs"); - - // note that we set in pararamcheck - // evolve_y_e and evolve_temper, but MoLRegister - // happens before paramcheck... - if (CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) { - register_constrained("HydroBase::Y_e"); - register_evolved("GRHydro::Y_e_con", - "GRHydro::Y_e_con_rhs"); - } - - if (CCTK_EQUALS(temperature_evolution_method,"GRHydro")) { - register_constrained("HydroBase::temperature"); - register_constrained("HydroBase::entropy"); - } - // particles - if (number_of_particles > 0) - register_evolved("GRHydro::particles", "GRHydro::particle_rhs"); - } - else if (CCTK_EQUALS(evolution_method, "none")) - { - register_constrained("GRHydro::dens"); - register_constrained("GRHydro::scon"); - register_constrained("GRHydro::tau"); - register_constrained("HydroBase::Bvec"); - if(clean_divergence) { - register_constrained("GRHydro::psidc"); - } - } -} - -

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 416)
by rhaas＠tapir.caltech.edu 27 Aug '12

27 Aug '12

User: rhaas Date: 2012/08/27 02:19 PM Modified: /trunk/ param.ccl /trunk/src/ GRHydro_PPM.F90, GRHydro_Reconstruct.F90 Log: GRHydro: Option to reconstruct W*vel (W = Lorentz factor) instead of just the primitive velocity 'vel'. This ensures that the Lorentz factor can never become unphysical during reconstruction. This patch is necessary for NS collapse when ePPM+Refluxing is used. The idea is due to Roland. From: Christian Reisswig <reisswig(a)scriwalker.(none)> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.F90 Delta lines: +147 -123 =================================================================== --- trunk/src/GRHydro_PPM.F90 2012-08-27 19:19:35 UTC (rev 415) +++ trunk/src/GRHydro_PPM.F90 2012-08-27 19:19:38 UTC (rev 416) @@ -106,7 +106,25 @@ logical :: cond + CCTK_REAL, dimension(nx) :: vx,vy,vz + + + !!$ 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). + vx = w_lorentz*velx + vy = w_lorentz*vely + vz = w_lorentz*velz + else + vx = velx + vy = vely + vz = velz + end if + !!$ Initially, all the Riemann problems will be trivial trivial_rp = .true. @@ -129,9 +147,9 @@ do i = 2, nx - 1 drho(i) = 0.5d0 * (rho(i+1) - rho(i-1)) - dvelx(i) = 0.5d0 * (velx(i+1) - velx(i-1)) - dvely(i) = 0.5d0 * (vely(i+1) - vely(i-1)) - dvelz(i) = 0.5d0 * (velz(i+1) - velz(i-1)) + dvelx(i) = 0.5d0 * (vx(i+1) - vx(i-1)) + dvely(i) = 0.5d0 * (vy(i+1) - vy(i-1)) + dvelz(i) = 0.5d0 * (vz(i+1) - vz(i-1)) dpress(i) = press(i+1) - press(i-1) d2rho(i) = (rho(i+1) - 2.d0 * rho(i) + rho(i-1))! / 6.d0 / dx / dx ! since we use d2rho only for the condition d2rho(i+1)*d2rhoi(i-1)<0 @@ -147,9 +165,9 @@ do i = 2, nx - 1 STEEP(rho, drho, dmrho) - STEEP(velx, dvelx, dmvelx) - STEEP(vely, dvely, dmvely) - STEEP(velz, dvelz, dmvelz) + STEEP(vx, dvelx, dmvelx) + STEEP(vy, dvely, dmvely) + STEEP(vz, dvelz, dmvelz) end do if (poly .eq. 0) then do i = 2, nx - 1 @@ -163,13 +181,13 @@ rhoplus(i) = 0.5d0 * (rho(i) + rho(i+1)) + & (dmrho(i) - dmrho(i+1)) / 6.d0 rhominus(i+1) = rhoplus(i) - velxplus(i) = 0.5d0 * (velx(i) + velx(i+1)) + & + velxplus(i) = 0.5d0 * (vx(i) + vx(i+1)) + & (dmvelx(i) - dmvelx(i+1)) / 6.d0 velxminus(i+1) = velxplus(i) - velyplus(i) = 0.5d0 * (vely(i) + vely(i+1)) + & + velyplus(i) = 0.5d0 * (vy(i) + vy(i+1)) + & (dmvely(i) - dmvely(i+1)) / 6.d0 velyminus(i+1) = velyplus(i) - velzplus(i) = 0.5d0 * (velz(i) + velz(i+1)) + & + velzplus(i) = 0.5d0 * (vz(i) + vz(i+1)) + & (dmvelz(i) - dmvelz(i+1)) / 6.d0 velzminus(i+1) = velzplus(i) end do @@ -235,16 +253,16 @@ LIMIT(rho, rhoplus(i), enhanced_ppm_C2, rhoplus(i)) rhominus(i+1) = rhoplus(i) - APPROX_AT_CELL_INTERFACE(velx, velxplus(i)) - LIMIT(velx, velxplus(i), enhanced_ppm_C2, velxplus(i)) + APPROX_AT_CELL_INTERFACE(vx, velxplus(i)) + LIMIT(vx, velxplus(i), enhanced_ppm_C2, velxplus(i)) velxminus(i+1) = velxplus(i) - APPROX_AT_CELL_INTERFACE(vely, velyplus(i)) - LIMIT(vely, velyplus(i), enhanced_ppm_C2, velyplus(i)) + APPROX_AT_CELL_INTERFACE(vy, velyplus(i)) + LIMIT(vy, velyplus(i), enhanced_ppm_C2, velyplus(i)) velyminus(i+1) = velyplus(i) - APPROX_AT_CELL_INTERFACE(velz, velzplus(i)) - LIMIT(velz, velzplus(i), enhanced_ppm_C2, velzplus(i)) + APPROX_AT_CELL_INTERFACE(vz, velzplus(i)) + LIMIT(vz, velzplus(i), enhanced_ppm_C2, velzplus(i)) velzminus(i+1) = velzplus(i) end do if (poly .eq. 0) then @@ -264,16 +282,16 @@ LIMIT(rho, rhoplus(i), enhanced_ppm_C2, rhoplus(i)) rhominus(i+1) = rhoplus(i) - APPROX_AT_CELL_INTERFACE_STENCIL4(velx, velxplus(i)) - LIMIT(velx, velxplus(i), enhanced_ppm_C2, velxplus(i)) + APPROX_AT_CELL_INTERFACE_STENCIL4(vx, velxplus(i)) + LIMIT(vx, velxplus(i), enhanced_ppm_C2, velxplus(i)) velxminus(i+1) = velxplus(i) - APPROX_AT_CELL_INTERFACE_STENCIL4(vely, velyplus(i)) - LIMIT(vely, velyplus(i), enhanced_ppm_C2, velyplus(i)) + APPROX_AT_CELL_INTERFACE_STENCIL4(vy, velyplus(i)) + LIMIT(vy, velyplus(i), enhanced_ppm_C2, velyplus(i)) velyminus(i+1) = velyplus(i) - APPROX_AT_CELL_INTERFACE_STENCIL4(velz, velzplus(i)) - LIMIT(velz, velzplus(i), enhanced_ppm_C2, velzplus(i)) + APPROX_AT_CELL_INTERFACE_STENCIL4(vz, velzplus(i)) + LIMIT(vz, velzplus(i), enhanced_ppm_C2, velzplus(i)) velzminus(i+1) = velzplus(i) end do if (poly .eq. 0) then @@ -444,7 +462,7 @@ !!$ Zone flattening. See appendix of C&W, p. 197-8. do i = 3, nx - 2 dpress2 = press(i+2) - press(i-2) - dvel = velx(i+1) - velx(i-1) + dvel = vx(i+1) - vx(i-1) if ( (abs(dpress(i)) > ppm_epsilon * min(press(i-1),press(i+1))) .and. & (dvel < 0.d0) ) then w = 1.d0 @@ -471,12 +489,12 @@ end if rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i) rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i) - velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * velx(i) - velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * velx(i) - velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vely(i) - velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vely(i) - velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * velz(i) - velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * velz(i) + velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i) + velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) + velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) + velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) + velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) + velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i) if (poly .eq. 0) then epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i) epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i) @@ -492,12 +510,12 @@ end if rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i) rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i) - velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * velx(i) - velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * velx(i) - velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vely(i) - velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vely(i) - velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * velz(i) - velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * velz(i) + velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i) + velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) + velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) + velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) + velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) + velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i) if (poly .eq. 0) then epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i) epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i) @@ -767,50 +785,14 @@ endif -! Reconstruct quantity on plus face -#define AVGP(a) \ - 0.5*(a(i)+a(i+1)) - -! Reconstruct quantity on minus face -#define AVGM(a) \ - 0.5*(a(i)+a(i-1)) - - -!! Compute scalar product v^i v_i on plus face -#define VEL2P(dvelx,dvely,dvelz) \ - (AVGP(gxx)*dvelx*dvelx + AVGP(gyy)*dvely*dvely + AVGP(gzz) \ - *dvelz*dvelz + 2*AVGP(gxy)*dvelx*dvely + 2*AVGP(gxz)*dvelx *dvelz + 2*AVGP(gyz) \ - *dvely*dvelz) - -!! Compute scalar product v^i v_i on minus face -#define VEL2M(dvelx,dvely,dvelz) \ - (AVGM(gxx)*dvelx*dvelx + AVGM(gyy)*dvely*dvely + AVGM(gzz) \ - *dvelz*dvelz + 2*AVGM(gxy)*dvelx*dvely + 2*AVGM(gxz)*dvelx *dvelz + 2*AVGM(gyz) \ - *dvely*dvelz) - - -!! Check if velocity is below speed of light. If not, reduce reconstruction to first order! -#define VELCHECK(vxminus,vx,vxplus,vyminus,vy,vyplus,vzminus,vz,vzplus) \ - if (VEL2M(vxminus(i),vyminus(i),vzminus(i)) .ge. (1.0d0-enhanced_ppm_C4) .or. VEL2P(vxplus(i),vyplus(i),vzplus(i)) .ge. (1.0d0-enhanced_ppm_C4)) then &&\ - vxminus(i) = vx(i) &&\ - vxplus(i) = vx(i) &&\ - vyminus(i) = vy(i) &&\ - vyplus(i) = vy(i) &&\ - vzminus(i) = vz(i) &&\ - vzplus(i) = vz(i) &&\ - endif - - if (use_enhanced_ppm .eq. 1) then !! Constrain parabolic profiles, PPM 2011/2008 if (PPM3) then do i = 3, nx - 2 MON_WITH_LOCAL_EXTREMUM(rhominus,rho,rhoplus, enhanced_ppm_C2) - MON_WITH_LOCAL_EXTREMUM(velxminus,velx,velxplus, enhanced_ppm_C2) - MON_WITH_LOCAL_EXTREMUM(velyminus,vely,velyplus, enhanced_ppm_C2) - MON_WITH_LOCAL_EXTREMUM(velzminus,velz,velzplus, enhanced_ppm_C2) - ! Ensure that velocity is below speed of light! - VELCHECK(velxminus,velx,velxplus,velyminus,vely,velyplus,velzminus,velz,velzplus) + MON_WITH_LOCAL_EXTREMUM(velxminus,vx,velxplus, enhanced_ppm_C2) + MON_WITH_LOCAL_EXTREMUM(velyminus,vy,velyplus, enhanced_ppm_C2) + MON_WITH_LOCAL_EXTREMUM(velzminus,vz,velzplus, enhanced_ppm_C2) end do if (poly .eq. 0) then do i = 3, nx - 2 @@ -821,11 +803,9 @@ else do i = 4, nx - 3 MON_WITH_LOCAL_EXTREMUM_STENCIL4(rhominus,rho,rhoplus, enhanced_ppm_C2, enhanced_ppm_C3) - MON_WITH_LOCAL_EXTREMUM_STENCIL4(velxminus,velx,velxplus, enhanced_ppm_C2, enhanced_ppm_C3) - MON_WITH_LOCAL_EXTREMUM_STENCIL4(velyminus,vely,velyplus, enhanced_ppm_C2, enhanced_ppm_C3) - MON_WITH_LOCAL_EXTREMUM_STENCIL4(velzminus,velz,velzplus, enhanced_ppm_C2, enhanced_ppm_C3) - ! Ensure that velocity is below speed of light! - VELCHECK(velxminus,velx,velxplus,velyminus,vely,velyplus,velzminus,velz,velzplus) + MON_WITH_LOCAL_EXTREMUM_STENCIL4(velxminus,vx,velxplus, enhanced_ppm_C2, enhanced_ppm_C3) + MON_WITH_LOCAL_EXTREMUM_STENCIL4(velyminus,vy,velyplus, enhanced_ppm_C2, enhanced_ppm_C3) + MON_WITH_LOCAL_EXTREMUM_STENCIL4(velzminus,vz,velzplus, enhanced_ppm_C2, enhanced_ppm_C3) end do if (poly .eq. 0) then do i = 4, nx - 3 @@ -844,12 +824,12 @@ end if rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i) rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i) - velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * velx(i) - velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * velx(i) - velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vely(i) - velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vely(i) - velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * velz(i) - velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * velz(i) + velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i) + velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) + velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) + velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) + velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) + velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i) if (poly .eq. 0) then epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i) epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i) @@ -865,12 +845,12 @@ end if rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i) rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i) - velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * velx(i) - velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * velx(i) - velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vely(i) - velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vely(i) - velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * velz(i) - velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * velz(i) + velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * vx(i) + velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) + velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) + velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) + velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) + velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(i) if (poly .eq. 0) then epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i) epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i) @@ -885,9 +865,9 @@ do i = GRHydro_stencil, nx - GRHydro_stencil + 1 ! original Colella&Woodward monotonicity preservation MON(rhominus,rho,rhoplus) - MON(velxminus,velx,velxplus) - MON(velyminus,vely,velyplus) - MON(velzminus,velz,velzplus) + MON(velxminus,vx,velxplus) + MON(velyminus,vy,velyplus) + MON(velzminus,vz,velzplus) end do if (poly .eq. 0) then do i = GRHydro_stencil, nx - GRHydro_stencil + 1 @@ -919,20 +899,20 @@ if (cond) then rhominus(i)=rho(i) rhoplus(i)=rho(i) - velxminus(i)=velx(i) - velxplus(i)=velx(i) - velyminus(i)=vely(i) - velyplus(i)=vely(i) - velzminus(i)=velz(i) - velzplus(i)=velz(i) + velxminus(i)=vx(i) + velxplus(i)=vx(i) + velyminus(i)=vy(i) + velyplus(i)=vy(i) + velzminus(i)=vz(i) + velzplus(i)=vz(i) rhominus(i-1)=rho(i) rhoplus(i-1)=rho(i) - velxminus(i-1)=velx(i) - velxplus(i-1)=velx(i) - velyminus(i-1)=vely(i) - velyplus(i-1)=vely(i) - velzminus(i-1)=velz(i) - velzplus(i-1)=velz(i) + velxminus(i-1)=vx(i) + velxplus(i-1)=vx(i) + velyminus(i-1)=vy(i) + velyplus(i-1)=vy(i) + velzminus(i-1)=vz(i) + velzplus(i-1)=vz(i) if (poly .eq. 0) then epsminus(i)=eps(i) epsplus(i)=eps(i) @@ -946,9 +926,9 @@ end if if (cond) then call PPM_TVD(rho(i-1), rho(i), rho(i+1), rhominus(i), rhoplus(i)) - call PPM_TVD(velx(i-1), velx(i), velx(i+1), velxminus(i), velxplus(i)) - call PPM_TVD(vely(i-1), vely(i), vely(i+1), velyminus(i), velyplus(i)) - call PPM_TVD(velz(i-1), velz(i), velz(i+1), velzminus(i), velzplus(i)) + call PPM_TVD(vx(i-1), vx(i), vx(i+1), velxminus(i), velxplus(i)) + call PPM_TVD(vy(i-1), vy(i), vy(i+1), velyminus(i), velyplus(i)) + call PPM_TVD(vz(i-1), vz(i), vz(i+1), velzminus(i), velzplus(i)) if (poly .eq. 0) then call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i)) end if @@ -961,20 +941,20 @@ if (cond) then rhominus(i)=rho(i) rhoplus(i)=rho(i) - velxminus(i)=velx(i) - velxplus(i)=velx(i) - velyminus(i)=vely(i) - velyplus(i)=vely(i) - velzminus(i)=velz(i) - velzplus(i)=velz(i) + velxminus(i)=vx(i) + velxplus(i)=vx(i) + velyminus(i)=vy(i) + velyplus(i)=vy(i) + velzminus(i)=vz(i) + velzplus(i)=vz(i) rhominus(i+1)=rho(i) rhoplus(i+1)=rho(i) - velxminus(i+1)=velx(i) - velxplus(i+1)=velx(i) - velyminus(i+1)=vely(i) - velyplus(i+1)=vely(i) - velzminus(i+1)=velz(i) - velzplus(i+1)=velz(i) + velxminus(i+1)=vx(i) + velxplus(i+1)=vx(i) + velyminus(i+1)=vy(i) + velyplus(i+1)=vy(i) + velzminus(i+1)=vz(i) + velzplus(i+1)=vz(i) if (poly .eq. 0) then epsminus(i)=eps(i) epsplus(i)=eps(i) @@ -988,9 +968,9 @@ end if if (cond) then call PPM_TVD(rho(i-1), rho(i), rho(i+1), rhominus(i), rhoplus(i)) - call PPM_TVD(velx(i-1), velx(i), velx(i+1), velxminus(i), velxplus(i)) - call PPM_TVD(vely(i-1), vely(i), vely(i+1), velyminus(i), velyplus(i)) - call PPM_TVD(velz(i-1), velz(i), velz(i+1), velzminus(i), velzplus(i)) + call PPM_TVD(vx(i-1), vx(i), vx(i+1), velxminus(i), velxplus(i)) + call PPM_TVD(vy(i-1), vy(i), vy(i+1), velyminus(i), velyplus(i)) + call PPM_TVD(vz(i-1), vz(i), vz(i+1), velzminus(i), velzplus(i)) if (poly .eq. 0) then call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i)) end if @@ -998,6 +978,50 @@ end if end if end do + + !!$ transform back to vel if Wv was used! + if (reconstruct_Wv.ne.0) then + do i = grhydro_stencil, nx - grhydro_stencil + 1 + ! 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) + agxx = 0.5d0*( gxx(i) + gxx(i-1) ) + agxy = 0.5d0*( gxy(i) + gxy(i-1) ) + agxz = 0.5d0*( gxz(i) + gxz(i-1) ) + agyy = 0.5d0*( gyy(i) + gyy(i-1) ) + agyz = 0.5d0*( gyz(i) + gyz(i-1) ) + agzz = 0.5d0*( gzz(i) + gzz(i-1) ) + w = sqrt( 1.d0 + agxx*velxminus(i)*velxminus(i) + & + agyy*velyminus(i)*velyminus(i) & + + agzz*velzminus(i)*velzminus(i) + & + 2.d0*agxy*velxminus(i)*velyminus(i) & + + 2.d0*agxz*velxminus(i)*velzminus(i) + & + 2.d0*agyz*velyminus(i)*velzminus(i) ) + velxminus(i) = velxminus(i)/w + velyminus(i) = velyminus(i)/w + velzminus(i) = velzminus(i)/w + + agxx = 0.5d0*( gxx(i) + gxx(i+1) ) + agxy = 0.5d0*( gxy(i) + gxy(i+1) ) + agxz = 0.5d0*( gxz(i) + gxz(i+1) ) + agyy = 0.5d0*( gyy(i) + gyy(i+1) ) + agyz = 0.5d0*( gyz(i) + gyz(i+1) ) + agzz = 0.5d0*( gzz(i) + gzz(i+1) ) + w = sqrt( 1.d0 + agxx*velxplus(i)*velxplus(i) + & + agyy*velyplus(i)*velyplus(i) & + + agzz*velzplus(i)*velzplus(i) + & + 2.d0*agxy*velxplus(i)*velyplus(i) & + + 2.d0*agxz*velxplus(i)*velzplus(i) + & + 2.d0*agyz*velyplus(i)*velzplus(i) ) + velxplus(i) = velxplus(i)/w + velyplus(i) = velyplus(i)/w + velzplus(i) = velzplus(i)/w + end do + end if + + return end subroutine SimplePPM_1d File [modified]: GRHydro_Reconstruct.F90 Delta lines: +68 -7 =================================================================== --- trunk/src/GRHydro_Reconstruct.F90 2012-08-27 19:19:35 UTC (rev 415) +++ trunk/src/GRHydro_Reconstruct.F90 2012-08-27 19:19:38 UTC (rev 416) @@ -36,6 +36,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 lvel, vel DECLARE_CCTK_ARGUMENTS @@ -45,14 +47,28 @@ CCTK_INT :: i,j,k CCTK_REAL :: local_min_tracer CCTK_INT :: type_bits, not_trivial + CCTK_REAL :: agxx, agxy, agxz, agyy, agyz, agzz, w ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc vup => lvel else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz vup => vel end if @@ -101,7 +117,7 @@ call CCTK_WARN(0, "Flux direction not x,y,z") end if - !$OMP PARALLEL DO PRIVATE(i,j,k) + !$OMP PARALLEL DO PRIVATE(i,j,k, agxx, agxy, agxz, agyy, agyz, agzz, w) do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 @@ -112,12 +128,57 @@ rhominus(i,j,k) = rho(i,j,k) epsplus(i,j,k) = eps(i,j,k) epsminus(i,j,k) = eps(i,j,k) - velxplus(i,j,k) = vup(i,j,k,1) - velyplus(i,j,k) = vup(i,j,k,2) - velzplus(i,j,k) = vup(i,j,k,3) - velxminus(i,j,k) = vup(i,j,k,1) - velyminus(i,j,k) = vup(i,j,k,2) - velzminus(i,j,k) = vup(i,j,k,3) + if (reconstruct_Wv.ne.0) then + ! divide out the Loretnz factor 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) + velxplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,1) + velyplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,2) + velzplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,3) + 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)) + w = 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)/w + velyplus(i,j,k) = velyplus(i,j,k)/w + velzplus(i,j,k) = velzplus(i,j,k)/w + + velxminus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,1) + velyminus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,2) + velzminus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,3) + 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)) + w = 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)/w + velyminus(i,j,k) = velyminus(i,j,k)/w + velzminus(i,j,k) = velzminus(i,j,k)/w + else + ! This is the standard way of doing it + velxplus(i,j,k) = vup(i,j,k,1) + velyplus(i,j,k) = vup(i,j,k,2) + velzplus(i,j,k) = vup(i,j,k,3) + velxminus(i,j,k) = vup(i,j,k,1) + velyminus(i,j,k) = vup(i,j,k,2) + velzminus(i,j,k) = vup(i,j,k,3) + endif if(evolve_y_e.ne.0) then y_e_plus(i,j,k) = y_e(i,j,k) y_e_minus(i,j,k) = y_e(i,j,k) Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +3 -0 =================================================================== --- trunk/param.ccl 2012-08-27 19:19:35 UTC (rev 415) +++ trunk/param.ccl 2012-08-27 19:19:38 UTC (rev 416) @@ -231,6 +231,9 @@ 0:1 :: "must be greater than or equal to 0, and not larger than 1" } 0.0 +boolean reconstruct_Wv "Reconstruct the primitive velocity W_Lorentz*vel, rather than just vel." STEERABLE=ALWAYS +{ +} no int eno_order "The order of accuracy of the ENO reconstruction" {

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 415)
by rhaas＠tapir.caltech.edu 27 Aug '12

27 Aug '12

User: rhaas Date: 2012/08/27 02:19 PM Modified: /trunk/src/ GRHydro_Con2Prim.F90, GRHydro_Con2PrimM.F90, GRHydro_HLLEM.F90, GRHydro_PPMMReconstruct_drv.F90, GRHydro_PPMReconstruct_drv.F90, GRHydro_Prim2ConM.F90, GRHydro_Reconstruct.F90, GRHydro_ReconstructPoly.F90, GRHydro_UpdateMaskM.F90 Log: GRHydro: Changes fixing issues with hot EOS treatment: 1) Set correct keytemp in C2P 2) Fix OpenMP private variable declaration in HLLEM 3) Fix velocity arguments for Y_e 1D PPM reconstruction 4) If in atmosphere and if evolve a Y_e also reset Y_e with cell averaged value in Reconstruction 5) Fix index permutation bug in ReconstructPoly routines for divergence cleaning field 6) Remove extra comment character in UpdateMaskM 7) Fix a missing velocity permutation when calling 1D reconstruction for Y_e. 8) Insert missing sqrt(det) factor for pressure term in resetting tau when temperature got too cold with a hot EOS. From: Philipp Moesta <pmoesta(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2012-08-27 19:19:32 UTC (rev 414) +++ trunk/src/GRHydro_Con2Prim.F90 2012-08-27 19:19:35 UTC (rev 415) @@ -327,7 +327,7 @@ call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),press(i,j,k),keyerr,anyerr) tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k) +press(i,j,k))& - * w_lorentz(i,j,k)**2 - dens(i,j,k) - press(i,j,k) + * w_lorentz(i,j,k)**2 - dens(i,j,k) - sqrt(det)*press(i,j,k) keytemp = 0 endif File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +3 -13 =================================================================== --- trunk/src/GRHydro_Con2PrimM.F90 2012-08-27 19:19:32 UTC (rev 414) +++ trunk/src/GRHydro_Con2PrimM.F90 2012-08-27 19:19:35 UTC (rev 415) @@ -238,17 +238,7 @@ uxx,uxy,uxz,uyy,uyz,uzz,det, & epsnegative,GRHydro_C2P_failed(i,j,k)) else - keytemp = 1 - !call Con2Prim_pt_hot(cctk_iteration,i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),& - ! scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),vup(i,j,k,1),& - ! vup(i,j,k,2),vup(i,j,k,3),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& - ! press(i,j,k),w_lorentz(i,j,k), & - ! uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & - ! z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & - ! GRHydro_reflevel, GRHydro_C2P_failed(i,j,k), GRHydro_perc_ptol) - ! Bprim(i,j,k,1) = 0.0d0 - ! Bprim(i,j,k,2) = 0.0d0 - ! Bprim(i,j,k,3) = 0.0d0 + keytemp = 0 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), & Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), Y_e(i,j,k), temperature(i,j,k), rho(i,j,k),& @@ -257,7 +247,7 @@ 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(GRHydro_C2P_failed(i,j,k).eq.2) then + if(GRHydro_C2P_failed(i,j,k).ne.0) then ! this means c2p did not converge. ! In this case, we attempt to call c2p with a reduced ! accuracy requirement; if it fails again, we abort @@ -271,7 +261,7 @@ 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(GRHydro_C2P_failed(i,j,k).eq.2) then + if(GRHydro_C2P_failed(i,j,k).ne.0) then !$OMP CRITICAL if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then call CCTK_WARN(1,"Convergence problem in c2p") File [modified]: GRHydro_HLLEM.F90 Delta lines: +4 -4 =================================================================== --- trunk/src/GRHydro_HLLEM.F90 2012-08-27 19:19:32 UTC (rev 414) +++ trunk/src/GRHydro_HLLEM.F90 2012-08-27 19:19:35 UTC (rev 415) @@ -133,7 +133,7 @@ !$OMP velzlowp,velzlowm,Bveczlowp,Bveczlowm,& !$OMP bdotvp,bdotvm,b2p,b2m,v2p,v2m,wp,wm,& !$OMP bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,& - !$OMP rhoenth_p,rhoenth_m,ab0p,ab0m,vA2p,vA2m,pressstarp,pressstarm,usendh,psidcdiff,psidcfp,psidcfm,charmin,charmax,chartop,charpm,m) + !$OMP rhoenth_p,rhoenth_m,ab0p,ab0m,vA2p,vA2m,pressstarp,pressstarm,usendh,psidcdiff,psidcfp,psidcfm,charmin,charmax,chartop,charpm,m,xtemp) do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + transport_constraints*(1-zoffset) do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + transport_constraints*(1-yoffset) do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + transport_constraints*(1-xoffset) @@ -362,7 +362,7 @@ lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) else - xtemp = temperature(i+xoffset,j+yoffset,k+zoffset) + 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(8),prim_m(9),prim_m(10),& @@ -414,7 +414,7 @@ lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) else - xtemp = temperature(i+xoffset,j+yoffset,k+zoffset) + 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(9),prim_m(10),prim_m(8),& @@ -466,7 +466,7 @@ lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) else - xtemp = temperature(i+xoffset,j+yoffset,k+zoffset) + 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(10),prim_m(8),prim_m(9),& File [modified]: GRHydro_PPMMReconstruct_drv.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/GRHydro_PPMMReconstruct_drv.F90 2012-08-27 19:19:32 UTC (rev 414) +++ trunk/src/GRHydro_PPMMReconstruct_drv.F90 2012-08-27 19:19:35 UTC (rev 415) @@ -318,7 +318,7 @@ do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + transport_constraints do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + transport_constraints call SimplePPM_ye_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & - velx(j,:,k),vely(j,:,k),velz(j,:,k), & + vely(j,:,k),velz(j,:,k),velx(j,:,k), & Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), & press(j,:,k)) end do @@ -391,7 +391,7 @@ do k = GRHydro_stencil, ny - GRHydro_stencil + 1 + transport_constraints do j = GRHydro_stencil, nx - GRHydro_stencil + 1 + transport_constraints call SimplePPM_ye_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & - velx(j,k,:),vely(j,k,:),velz(j,k,:), & + velz(j,k,:),velx(j,k,:),vely(j,k,:), & Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(j,k,:), & press(j,k,:)) end do File [modified]: GRHydro_PPMReconstruct_drv.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv.F90 2012-08-27 19:19:32 UTC (rev 414) +++ trunk/src/GRHydro_PPMReconstruct_drv.F90 2012-08-27 19:19:35 UTC (rev 415) @@ -246,7 +246,7 @@ do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 call SimplePPM_ye_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & - velx(j,:,k),vely(j,:,k),velz(j,:,k), & + vely(j,:,k),velz(j,:,k),velx(j,:,k), & Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), & press(j,:,k)) end do @@ -298,7 +298,7 @@ do k = GRHydro_stencil, ny - GRHydro_stencil + 1 do j = GRHydro_stencil, nx - GRHydro_stencil + 1 call SimplePPM_ye_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & - velx(j,k,:),vely(j,k,:),velz(j,k,:), & + velz(j,k,:),velx(j,k,:),vely(j,k,:), & Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(j,k,:), & press(j,k,:)) end do File [modified]: GRHydro_Prim2ConM.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_Prim2ConM.F90 2012-08-27 19:19:32 UTC (rev 414) +++ trunk/src/GRHydro_Prim2ConM.F90 2012-08-27 19:19:35 UTC (rev 415) @@ -77,7 +77,7 @@ if(evolve_temper.ne.1) then !$OMP PARALLEL DO PRIVATE(k,j,i,avg_detl,avg_detr,& - !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,& !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) File [modified]: GRHydro_Reconstruct.F90 Delta lines: +14 -10 =================================================================== --- trunk/src/GRHydro_Reconstruct.F90 2012-08-27 19:19:32 UTC (rev 414) +++ trunk/src/GRHydro_Reconstruct.F90 2012-08-27 19:19:35 UTC (rev 415) @@ -118,17 +118,21 @@ velxminus(i,j,k) = vup(i,j,k,1) velyminus(i,j,k) = vup(i,j,k,2) velzminus(i,j,k) = vup(i,j,k,3) + if(evolve_y_e.ne.0) then + y_e_plus(i,j,k) = y_e(i,j,k) + y_e_minus(i,j,k) = y_e(i,j,k) + endif + if(evolve_tracer.ne.0) then + where(tracerplus(i,j,k,:).le.local_min_tracer .or. & + tracerminus(i,j,k,:).le.local_min_tracer) + tracerplus(i,j,k,:) = tracer(i,j,k,:) + tracerminus(i,j,k,:) = tracer(i,j,k,:) + end where + end if end if - if(evolve_tracer.ne.0) then - where(tracerplus(i,j,k,:).le.local_min_tracer .or. & - tracerminus(i,j,k,:).le.local_min_tracer) - tracerplus(i,j,k,:) = tracer(i,j,k,:) - tracerminus(i,j,k,:) = tracer(i,j,k,:) - end where - end if - ! Riemann problem might not be trivial anymore!!!! - SpaceMask_SetStateBitsF90(space_mask, i-xoffset, j-yoffset, k-zoffset, type_bits, not_trivial) - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, not_trivial) + ! Riemann problem might not be trivial anymore!!!! + SpaceMask_SetStateBitsF90(space_mask, i-xoffset, j-yoffset, k-zoffset, type_bits, not_trivial) + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, not_trivial) enddo enddo enddo File [modified]: GRHydro_ReconstructPoly.F90 Delta lines: +5 -5 =================================================================== --- trunk/src/GRHydro_ReconstructPoly.F90 2012-08-27 19:19:32 UTC (rev 414) +++ trunk/src/GRHydro_ReconstructPoly.F90 2012-08-27 19:19:35 UTC (rev 415) @@ -496,11 +496,11 @@ do j = GRHydro_stencil, nx - GRHydro_stencil + 1 call SimplePPM_1dM(GRHydro_eos_handle,1,nz,CCTK_DELTA_SPACE(3),& rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),& - Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),psidc(:,j,k),eps(j,k,:),press(j,k,:),& + Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),psidc(j,k,:),eps(j,k,:),press(j,k,:),& rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),& - Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),psidcminus(:,j,k),epsminus(j,k,:),& + Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),psidcminus(j,k,:),epsminus(j,k,:),& rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& - Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),psidcplus(:,j,k),epsplus(j,k,:),& + Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),psidcplus(j,k,:),epsplus(j,k,:),& 1,trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& @@ -526,9 +526,9 @@ rho(j,k,:),velz(j,k,:),velx(j,k,:),vely(j,k,:),& Bvecz(j,k,:),Bvecx(j,k,:),Bvecy(j,k,:),dum(:,j,k),eps(j,k,:),press(j,k,:),& rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),& - Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),dumm(:,j,k),epsminus(j,k,:),& + Bveczminus(j,k,:),Bvecxminus(j,k,:),Bvecyminus(j,k,:),dumm(j,k,:),epsminus(j,k,:),& rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& - Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),dump(:,j,k),epsplus(j,k,:),& + Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),dump(j,k,:),epsplus(j,k,:),& 0,trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& File [modified]: GRHydro_UpdateMaskM.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_UpdateMaskM.F90 2012-08-27 19:19:32 UTC (rev 414) +++ trunk/src/GRHydro_UpdateMaskM.F90 2012-08-27 19:19:35 UTC (rev 415) @@ -123,7 +123,7 @@ if(evolve_temper.ne.0) then -! ! set the temperature to be relatively low + ! set the temperature to be relatively low temperature(i,j,k) = grhydro_hot_atmo_temp y_e(i,j,k) = grhydro_hot_atmo_Y_e keytemp = 1

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits