Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 72)
by knarf＠cct.lsu.edu 27 Apr '11

27 Apr '11

User: knarf Date: 2011/04/27 05:55 PM Modified: / ET.tex Log: make clear that we should put all references in one file, using one convention File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +1 -0 =================================================================== --- ET.tex 2011-04-27 18:57:16 UTC (rev 71) +++ ET.tex 2011-04-27 22:55:25 UTC (rev 72) @@ -1889,6 +1889,7 @@ project cfz-411-aa}. \bibliographystyle{amsplain-url} +% TODO: all references should be in manifest/einsteintoolkit \bibliography{manifest/einsteintoolkit,local_bibtex/ott_references,local_bibtex/references} \end{document}

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 125)
by bcmsma＠astro.rit.edu 27 Apr '11

27 Apr '11

User: bmundim Date: 2011/04/27 05:37 PM Modified: /trunk/ interface.ccl Log: Change Prim2Con functions interface. File Changes: Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +17 -14 =================================================================== --- trunk/interface.ccl 2011-04-27 22:26:38 UTC (rev 124) +++ trunk/interface.ccl 2011-04-27 22:37:50 UTC (rev 125) @@ -40,18 +40,19 @@ CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz) SUBROUTINE Prim2ConPolyM(CCTK_INT IN handle, \ - CCTK_REAL IN gxx, CCTK_REAL IN gxy, CCTK_REAL IN gxz, \ - CCTK_REAL IN gyy, CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ - CCTK_REAL IN det, CCTK_REAL OUT dens, \ - CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ - CCTK_REAL OUT sz, CCTK_REAL OUT tau, \ - CCTK_REAL IN Bvecx, CCTK_REAL IN Bvecy, \ - CCTK_REAL IN Bvecz, CCTK_REAL IN rho, CCTK_REAL IN velx, \ - CCTK_REAL IN vely, \ - CCTK_REAL IN velz, CCTK_REAL OUT epsilon, \ - CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz) + CCTK_REAL IN gxx, CCTK_REAL IN gxy, \ + CCTK_REAL IN gxz, CCTK_REAL IN gyy, \ + CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ + CCTK_REAL IN det, CCTK_REAL OUT dens, \ + CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ + CCTK_REAL OUT sz, CCTK_REAL OUT tau, \ + CCTK_REAL OUT Bconsx, CCTK_REAL OUT Bconsy, \ + CCTK_REAL OUT Bconsz, CCTK_REAL IN rho, CCTK_REAL IN velx, \ + CCTK_REAL IN vely, \ + CCTK_REAL IN velz, CCTK_REAL OUT epsilon, \ + CCTK_REAL OUT press, CCTK_REAL IN Bvecx, CCTK_REAL IN Bvecy, \ + CCTK_REAL IN Bvecz, CCTK_REAL OUT w_lorentz) - SUBROUTINE Prim2ConGenM(CCTK_INT IN handle, \ CCTK_REAL IN gxx, CCTK_REAL IN gxy, \ CCTK_REAL IN gxz, CCTK_REAL IN gyy, \ @@ -59,12 +60,14 @@ CCTK_REAL IN det, CCTK_REAL OUT dens, \ CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ CCTK_REAL OUT sz, CCTK_REAL OUT tau, \ - CCTK_REAL IN Bvecx, CCTK_REAL IN Bvecy, \ - CCTK_REAL IN Bvecz, CCTK_REAL IN rho, CCTK_REAL IN velx, \ + CCTK_REAL OUT Bconsx, CCTK_REAL OUT Bconsy, \ + CCTK_REAL OUT Bconsz, CCTK_REAL IN rho, CCTK_REAL IN velx, \ CCTK_REAL IN vely, \ CCTK_REAL IN velz, CCTK_REAL IN epsilon, \ - CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz) + CCTK_REAL OUT press, CCTK_REAL IN Bvecx, CCTK_REAL IN Bvecy, \ + CCTK_REAL IN Bvecz, CCTK_REAL OUT w_lorentz) + SUBROUTINE Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL OUT dens, \ CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ CCTK_REAL OUT sz, CCTK_REAL OUT tau, \

1 0

[svn:einsteintoolkit] GRHydro/branches/divclean/ (Rev. 234)
by tanja.bode＠physics.gatech.edu 27 Apr '11

27 Apr '11

User: tbode Date: 2011/04/27 05:37 PM Modified: /branches/divclean/ interface.ccl Log: fix typo. File Changes: Directory: /branches/divclean/ ============================== File [modified]: interface.ccl Delta lines: +1 -1 =================================================================== --- branches/divclean/interface.ccl 2011-04-27 22:20:36 UTC (rev 233) +++ branches/divclean/interface.ccl 2011-04-27 22:37:19 UTC (rev 234) @@ -424,7 +424,7 @@ Bconsxplus,Bconsyplus,Bconszplus,Bconsxminus,Bconsyminus,Bconszminus } "Conservative variables extended to the cell boundaries" -real GRHydro_MHD_prin_bext type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' +real GRHydro_MHD_prim_bext type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' { Bvecxplus,Bvecyplus,Bveczplus,Bvecxminus,Bvecyminus,Bveczminus } "Primitive mhd variables extended to the cell boundaries"

1 0

[svn:einsteintoolkit] manifest/trunk/ (Rev. 95)
by knarf＠cct.lsu.edu 27 Apr '11

27 Apr '11

User: knarf Date: 2011/04/27 05:33 PM Modified: /trunk/ einsteintoolkit.th Log: make the hg-version of Carpet standard in the dev-version of the ET toolkit File Changes: Directory: /trunk/ ================== File [modified]: einsteintoolkit.th Delta lines: +8 -7 =================================================================== --- trunk/einsteintoolkit.th 2011-04-21 19:59:51 UTC (rev 94) +++ trunk/einsteintoolkit.th 2011-04-27 22:33:23 UTC (rev 95) @@ -288,13 +288,14 @@ # Carpet, the AMR driver !TARGET = $ARR -!TYPE = git -!AUTH_URL = carpetgit@carpetcode.org:carpet -!URL = git://carpetcode.org/carpet -# For experimental version, comment previous three lines and uncomment next two -# Note that you will need to get access first (ask Erik Schnetter) -#!TYPE = hg -#!URL = ssh://carpetmercurial@carpetcode.org/carpet +!TYPE = hg +!AUTH_URL = ssh://carpetmercurial@carpetcode.org/carpet +!URL = http://www.carpetcode.org/hg/carpet +# For the git version of Carpet, comment the previous three lines and uncomment +# the next three +#!TYPE = git +#!AUTH_URL = carpetgit@carpetcode.org:carpet +#!URL = git://carpetcode.org/carpet !CHECKOUT = Carpet/doc Carpet/Carpet Carpet/CarpetEvolutionMask

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 124)
by bcmsma＠astro.rit.edu 27 Apr '11

27 Apr '11

User: bmundim Date: 2011/04/27 05:26 PM Added: /trunk/src/ GRHydro_CylindricalExplosionM.F90, GRHydro_MonopoleM.F90 Log: RIT MHD dev: 1) Komissarov's cylindrical explosion test 2) Magnetic monopole test File Changes: Directory: /trunk/src/ ====================== File [added]: GRHydro_CylindricalExplosionM.F90 Delta lines: +164 -0 =================================================================== --- trunk/src/GRHydro_CylindricalExplosionM.F90 (rev 0) +++ trunk/src/GRHydro_CylindricalExplosionM.F90 2011-04-27 22:26:38 UTC (rev 124) @@ -0,0 +1,164 @@ + /*@@ + @file GRHydro_CylindricalExplosionM.F90 + @date Apr 22, 2011 + @author Scott Noble, Joshua Faber, Bruno Mundim + @desc + Cylindrical magnetized shocks. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "GRHydro_Macros.h" + +#define velx(i,j,k) vel(i,j,k,1) +#define vely(i,j,k) vel(i,j,k,2) +#define velz(i,j,k) vel(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) + + + /*@@ + @routine GRHydro_cylindricalexplosionM + @date Fri Apr 22 14:51:37 EDT 2011 + @author Scott Noble, Joshua Faber, Bruno Mundim + @desc + Initial data for cylidrically symmetric shocks with + magnetic fields. Meant to recreate tests demonstrated + in Komissarov (1999). + @enddesc + @calls + @calledby + @history + Using GRHydro_ShockTubeM.F90 as a template. + @endhistory + +@@*/ + +subroutine GRHydro_cylindricalexplosionM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, nx, ny, nz + CCTK_REAL :: direction, det + CCTK_REAL :: rhol, rhor, pressl, pressr + CCTK_REAL :: bvcxl,bvcyl,bvczl + CCTK_REAL :: tmp,tmp2,gam,r_inner,r_outer + + CCTK_REAL :: cyl_fr + +!!$Uses Bx_init, By_init and Bz_init to set magnetic field strength +!!$Original tests had Bx_init = 0.1, 1.0 with By_init=Bz_init = 0 + + bvcxl = Bx_init + bvcyl = By_init + bvczl = Bz_init + +!!$Adiabatic index for test + gam = (4.d0/3.d0) + +!!$Inner radius and outer radius + r_inner = 8.d-1 + r_outer = 1.d0 + +!!$Inner values + rhol = 1.d-4 + pressl = 3.d-5 + +!!$Outer values + rhor = 1.d-2 + pressr = 1.d0 + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + do i=1,nx + do j=1,ny + do k=1,nz + +!!$direction represents the cylindrical radius here. + if (CCTK_EQUALS(shocktube_type,"xshock")) then + direction = sqrt((x(i,j,k)-shock_xpos)**2+& + (y(i,j,k)-shock_ypos)**2) + else if (CCTK_EQUALS(shocktube_type,"yshock")) then + direction = sqrt((y(i,j,k)-shock_ypos)**2+& + (z(i,j,k)-shock_zpos)**2) + else if (CCTK_EQUALS(shocktube_type,"zshock")) then + direction = sqrt((x(i,j,k)-shock_xpos)**2+& + (z(i,j,k)-shock_zpos)**2) + end if + + tmp = cyl_fr(direction,r_inner,r_outer,rhol,rhor) + tmp2 = cyl_fr(direction,r_inner,r_outer,pressl,pressr)/( (gam - 1.d0) * tmp ) + rho(i,j,k) = tmp + eps(i,j,k) = tmp2 + + velx(i,j,k) = 0.d0 + vely(i,j,k) = 0.d0 + velz(i,j,k) = 0.d0 + Bvecx(i,j,k)=bvcxl + Bvecy(i,j,k)=bvcyl + Bvecz(i,j,k)=bvczl + + det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + + if (CCTK_EQUALS(GRHydro_eos_type,"Polytype")) then + call Prim2ConPolyM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + else + call Prim2ConGenM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + end if + enddo + enddo + enddo + + densrhs = 0.d0 + srhs = 0.d0 + taurhs = 0.d0 + Bvecrhs = 0.d0 + + + return + +end subroutine GRHydro_CylindricalExplosionM + + + +function cyl_fr(R,r_inner,r_outer,min_f,max_f) + + implicit none + + CCTK_REAL :: cyl_fr + CCTK_REAL :: R,r_inner,r_outer,min_f,max_f + + if(R .gt. r_outer) then + cyl_fr = min_f + else if( R .lt. r_inner) then + cyl_fr = max_f + else + cyl_fr = exp( ( log(max_f)*(r_outer - R) + log(min_f)*(R - r_inner) ) / (r_outer - r_inner) ) + end if + + return + +end function cyl_fr File [added]: GRHydro_MonopoleM.F90 Delta lines: +109 -0 =================================================================== --- trunk/src/GRHydro_MonopoleM.F90 (rev 0) +++ trunk/src/GRHydro_MonopoleM.F90 2011-04-27 22:26:38 UTC (rev 124) @@ -0,0 +1,109 @@ + /*@@ + @file GRHydro_MonopoleM.F90 + @date Sep 23, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke + @desc + Initial data of the shock tube type - MHD version. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "GRHydro_Macros.h" + +#define velx(i,j,k) vel(i,j,k,1) +#define vely(i,j,k) vel(i,j,k,2) +#define velz(i,j,k) vel(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) + + /*@@ + @routine GRHydro_MonopoleM + @date Sat Jan 26 02:53:49 2002 + @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke + @desc + A monopole in space + @enddesc + @calls + @calledby + @history + Expansion and alteration of the test code from GRAstro_Hydro, + written by Mark Miller. + @endhistory + +@@*/ + +subroutine GRHydro_MonopoleM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, nx, ny, nz + CCTK_REAL :: direction, det + CCTK_REAL :: rhol, rhor, velxl, velxr, velyl, velyr, & + velzl, velzr, epsl, epsr + CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr + CCTK_REAL :: ux,uy,uz,ut,tmp,tmp2,tmp3 + + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + do i=1,nx + do j=1,ny + do k=1,nz + + rho(i,j,k) = 1.0 + velx(i,j,k) = 0.0 + vely(i,j,k) = 0.0 + velz(i,j,k) = 0.0 + eps(i,j,k) = 0.1 + if(i.eq.nx/2+1.and.j.eq.ny/2+1.and.k.eq.nz/2+1) then + Bvecx(i,j,k)=0.1 + else + Bvecx(i,j,k)=0.0 + endif + Bvecy(i,j,k)=0.0 + Bvecz(i,j,k)=0.0 + + det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + + if (CCTK_EQUALS(GRHydro_eos_type,"Polytype")) then + call Prim2ConPolyM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + else + call Prim2ConGenM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + end if + + enddo + enddo + enddo + + densrhs = 0.d0 + srhs = 0.d0 + taurhs = 0.d0 + Bvecrhs = 0.d0 + + + return + +end subroutine GRHydro_MonopoleM

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 123)
by bcmsma＠astro.rit.edu 27 Apr '11

27 Apr '11

User: bmundim Date: 2011/04/27 05:25 PM Modified: /trunk/ param.ccl, schedule.ccl /trunk/par/ balsara1a.par /trunk/src/ CheckParam.c, GRHydro_ShockTube.F90, GRHydro_ShockTubeM.F90, make.code.defn Log: RIT MHD dev: Balsara's and Komissarov's shock tube tests. File Changes: Directory: /trunk/par/ ====================== File [modified]: balsara1a.par Delta lines: +12 -6 =================================================================== --- trunk/par/balsara1a.par 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/par/balsara1a.par 2011-04-27 22:25:23 UTC (rev 123) @@ -5,7 +5,8 @@ driver::ghost_size=3 grhydro::grhydro_stencil=3 -time::dtfac = 0.25 +#time::dtfac = 0.25 +time::dtfac = 0.8 methodoflines::ODE_Method = "rk2" methodoflines::MoL_Intermediate_Steps=2 @@ -37,25 +38,30 @@ grid::type = "BySpacing" grid::domain = "full" -grid::dxyz = 0.01 +#grid::dxyz = 0.01 +grid::dxyz = 0.000625 -driver::global_nx = 150 +driver::global_nx = 1600 driver::global_ny = 10 driver::global_nz = 10 +Cactus::terminate="time" +Cactus::cctk_final_time = 0.4 #cactus::cctk_itlast = 40 -cactus::cctk_itlast = 200 +#cactus::cctk_itlast = 200 +#cactus::cctk_itlast = 2 -IO::out_dir = "balsara1a" +IO::out_dir = $parfile #IOBasic::outInfo_every = 1 #IOBasic::outInfo_vars = "HydroBase::rho" CarpetIOBasic::outInfo_vars="hydrobase::rho" CarpetIOBasic::outInfo_every=1 CarpetIOASCII::out1D_criterion = "time" CarpetIOASCII::out1D_dt = 0.01 +CarpetIOASCII::out1D_d=no CarpetIOASCII::out1D_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec" CarpetIOHDF5::out_every = 1 -CarpetIOHDF5::out_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec" +CarpetIOHDF5::out_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel HydroBase::w_lorentz grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec" Directory: /trunk/src/ ====================== File [modified]: CheckParam.c Delta lines: +18 -7 =================================================================== --- trunk/src/CheckParam.c 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/src/CheckParam.c 2011-04-27 22:25:23 UTC (rev 123) @@ -18,20 +18,31 @@ (CCTK_Equals(initial_hydro,"simple_wave")) || (CCTK_Equals(initial_data,"con2primtest")) || (CCTK_Equals(initial_data,"reconstruction_test")) || - (CCTK_Equals(shocktube_type,"diagshock")) || (CCTK_Equals(shocktube_type,"sphere")))) { CCTK_PARAMWARN("That test not yet implemented in MHD!"); } if(!CCTK_Equals(Bvec_evolution_method,"GRHydro") && - ((CCTK_Equals(shock_case,"Balsara1")) || - (CCTK_Equals(shock_case,"Balsara2")) || - (CCTK_Equals(shock_case,"Balsara3")) || - (CCTK_Equals(shock_case,"Balsara4")) || - (CCTK_Equals(shock_case,"Balsara5")))) + ((CCTK_Equals(shock_case,"Balsara0" )) || + (CCTK_Equals(shock_case,"Balsara1" )) || + (CCTK_Equals(shock_case,"Balsara2" )) || + (CCTK_Equals(shock_case,"Balsara3" )) || + (CCTK_Equals(shock_case,"Balsara4" )) || + (CCTK_Equals(shock_case,"Balsara5" )) || + (CCTK_Equals(shock_case,"Alfven" )) || + (CCTK_Equals(shock_case,"Komissarov1")) || + (CCTK_Equals(shock_case,"Komissarov2")) || + (CCTK_Equals(shock_case,"Komissarov3")) || + (CCTK_Equals(shock_case,"Komissarov4")) || + (CCTK_Equals(shock_case,"Komissarov5")) || + (CCTK_Equals(shock_case,"Komissarov6")) || + (CCTK_Equals(shock_case,"Komissarov7")) || + (CCTK_Equals(shock_case,"Komissarov8")) || + (CCTK_Equals(shock_case,"Komissarov9")) || + (CCTK_Equals(initial_hydro,"cylexp")) + )) { CCTK_PARAMWARN("That test requires MHD! Set Bvec_evolution_method=GRHYDRO!"); } } - File [modified]: GRHydro_ShockTube.F90 Delta lines: +11 -0 =================================================================== --- trunk/src/GRHydro_ShockTube.F90 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/src/GRHydro_ShockTube.F90 2011-04-27 22:25:23 UTC (rev 123) @@ -93,6 +93,17 @@ velzr = 0.d0 epsl = 1.5d0 epsr = 1.2d0 + else if (CCTK_EQUALS(shock_case,"Balsaralike1")) then + rhol = 1.0d0 + rhor = 0.125d0 + velxl = 0.0d0 + velxr = 0.0d0 + velyl = 0.0d0 + velyr = 0.0d0 + velzl = 0.0d0 + velzr = 0.0d0 + epsl = 1.0d0/rhol + epsr = 0.1d0/rhor !!$This line only for polytrope, k=1 !!$ epsr = 0.375d0 else if (CCTK_EQUALS(shock_case,"Blast")) then File [modified]: GRHydro_ShockTubeM.F90 Delta lines: +565 -8 =================================================================== --- trunk/src/GRHydro_ShockTubeM.F90 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/src/GRHydro_ShockTubeM.F90 2011-04-27 22:25:23 UTC (rev 123) @@ -23,6 +23,11 @@ #define Bvecy(i,j,k) Bvec(i,j,k,2) #define Bvecz(i,j,k) Bvec(i,j,k,3) +#define OOSQRT2 (0.7071067811865475244008442) +#define OOSQRT3 (0.5773502691896257645091489) +#define OOSQRT6 (0.4082482904638630163662140) + + /*@@ @routine GRHydro_shocktubeM @date Sat Jan 26 02:53:49 2002 @@ -53,7 +58,9 @@ CCTK_REAL :: rhol, rhor, velxl, velxr, velyl, velyr, & velzl, velzr, epsl, epsr CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr + CCTK_REAL :: ux,uy,uz,ut,tmp,tmp2,tmp3 + bvcxl = Bx_init bvcyl = By_init bvczl = Bz_init @@ -69,18 +76,27 @@ do j=1,ny do k=1,nz if (CCTK_EQUALS(shocktube_type,"diagshock")) then - direction = x(i,j,k) - shock_xpos + & - y(i,j,k) - shock_ypos + z(i,j,k) - shock_zpos +!!$ The diagshock choice yields a shock plane perpendicular to the fixed vector (1,1,1) +!!$ This could be changed, but would require 3 new params containing the new shock direction + direction = x(i,j,k) - shock_xpos + & + y(i,j,k) - shock_ypos + z(i,j,k) - shock_zpos + + else if (CCTK_EQUALS(shocktube_type,"diagshock2d")) then +!!$ The diagshock choice yields a shock plane perpendicular to the fixed vector (1,1,0), with similarity in the z-dir. +!!$ This could be changed, but would require 2 new params containing the new shock direction + direction = x(i,j,k) - shock_xpos + & + y(i,j,k) - shock_ypos + else if (CCTK_EQUALS(shocktube_type,"xshock")) then - direction = x(i,j,k) - shock_xpos + direction = x(i,j,k) - shock_xpos else if (CCTK_EQUALS(shocktube_type,"yshock")) then - direction = y(i,j,k) - shock_ypos + direction = y(i,j,k) - shock_ypos else if (CCTK_EQUALS(shocktube_type,"zshock")) then - direction = z(i,j,k) - shock_zpos + direction = z(i,j,k) - shock_zpos else if (CCTK_EQUALS(shocktube_type,"sphere")) then - direction = sqrt((x(i,j,k)-shock_xpos)**2+& - (y(i,j,k)-shock_ypos)**2+& - (z(i,j,k)-shock_zpos)**2)-shock_radius + direction = sqrt((x(i,j,k)-shock_xpos)**2+& + (y(i,j,k)-shock_ypos)**2+& + (z(i,j,k)-shock_zpos)**2)-shock_radius end if if (CCTK_EQUALS(shock_case,"Simple")) then rhol = 10.d0 @@ -117,6 +133,31 @@ velzr = 0.d0 epsl = 1500.d0 epsr = 1.5d-2 + +!!$ The following shocktubes are from Balsara 2001 . +!!$ All use n=1600 cells, over domain x=[-0.5,0.5] +!!$ All assume ideal-gas or gamma-law EOS, the first test uses GAMMA=2. while the rest GAMMA=5./3. + +!!$ Unmagnetized Test 1 (rel. Brio & Wu 1988 by van Putten 1993) of Balsara 2001 -- compare at t=0.4 + else if (CCTK_EQUALS(shock_case,"Balsara0")) then + rhol = 1.0d0 + rhor = 0.125d0 + velxl = 0.0d0 + velxr = 0.0d0 + velyl = 0.0d0 + velyr = 0.0d0 + velzl = 0.0d0 + velzr = 0.0d0 + bvcxl=0.0d0 + bvcxr=0.0d0 + bvcyl=0.0d0 + bvcyr=0.0d0 + bvczl=0.0d0 + bvczr=0.0d0 + epsl = 1.0d0/rhol + epsr = 0.1d0/rhor + +!!$ Test 1 (rel. Brio & Wu 1988 by van Putten 1993) of Balsara 2001 -- compare at t=0.4 else if (CCTK_EQUALS(shock_case,"Balsara1")) then rhol = 1.0d0 rhor = 0.125d0 @@ -134,6 +175,8 @@ bvczr=0.0d0 epsl = 1.0d0/rhol epsr = 0.1d0/rhor + +!!$ Test 2 (blast wave) of Balsara 2001 -- compare at t=0.4 else if (CCTK_EQUALS(shock_case,"Balsara2")) then rhol = 1.d0 rhor = 1.d0 @@ -151,6 +194,8 @@ bvczr=0.7d0 epsl = 1.5d0*30.0d0/rhol epsr = 1.5d0*1.0d0/rhor + +!!$ Test 3 (blast wave) of Balsara 2001 -- compare at t=0.4 else if (CCTK_EQUALS(shock_case,"Balsara3")) then rhol = 1.d0 rhor = 1.d0 @@ -168,6 +213,8 @@ bvczr=0.7d0 epsl = 1.5d0*1000.0d0/rhol epsr = 1.5d0*0.1d0/rhor + +!!$ Test 4 (rel. version of Noh 1987) of Balsara 2001 -- compare at t=0.4 else if (CCTK_EQUALS(shock_case,"Balsara4")) then rhol = 1.d0 rhor = 1.d0 @@ -185,6 +232,8 @@ bvczr=-7.d0 epsl = 1.5d0*0.1d0/rhol epsr = 1.5d0*0.1d0/rhor + +!!$ Test 5 (non-coplanar set of waves) of Balsara 2001 -- compare at t=0.55 else if (CCTK_EQUALS(shock_case,"Balsara5")) then rhol = 1.08d0 rhor = 1.d0 @@ -202,12 +251,286 @@ bvczr=0.5d0 epsl = 1.5d0*0.95d0/rhol epsr = 1.5d0*1.0d0/rhor + +!!$ "Generic Alfven Test of Giacomazzo and Rezzolla J.Comp.Phys (2006) + else if (CCTK_EQUALS(shock_case,"Alfven")) then + rhol = 1.d0 + rhor = 0.9d0 + velxl = 0.d0 + velxr = 0.d0 + velyl = 0.3d0 + velyr = 0.d0 + velzl = 0.4d0 + velzr = 0.d0 + bvcxl=1.0d0 + bvcxr=1.0d0 + bvcyl=6.d0 + bvcyr=5.d0 + bvczl=2.d0 + bvczr=2.d0 + epsl = 1.5d0*5.d0/rhol + epsr = 1.5d0*5.3d0/rhor + +!!$ The following 9 tests are from Komissarov 1999 . +!!$ Note that the data is specified in terms of the 4-velocity, so some conversion is necessary. +!!$ All assume ideal-gas or gamma-law EOS with GAMMA=4./3. + +!!$ Fast Shock test of Komissarov 1999 -- compare at t=2.5 , n=40, x=[-1,1] + else if (CCTK_EQUALS(shock_case,"Komissarov1")) then + rhol = 1.d0 + rhor = 25.48d0 + bvcxl=20.d0 + bvcxr=20.d0 + bvcyl=25.02d0 + bvcyr=49.d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d0 /rhol + epsr = 3.d0 * 367.5d0 /rhor + + ux=25.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=1.091d0 + uy=0.3923d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + +!!$ Slow Shock test of Komissarov 1999 -- compare at t=2. , n=200, x=[-0.5,1.5] + else if (CCTK_EQUALS(shock_case,"Komissarov2")) then + rhol = 1.d0 + rhor = 3.323d0 + bvcxl=10.d0 + bvcxr=10.d0 + bvcyl=18.28d0 + bvcyr=14.49d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d0 /rhol + epsr = 3.d0 * 55.36d0 /rhor + + ux=1.53d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=0.9571d0 + uy=-0.6822d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + +!!$ Switch-off Fast test of Komissarov 1999 -- compare at t=1. , n=150, x=[-1,1] + else if (CCTK_EQUALS(shock_case,"Komissarov3")) then + rhol = 0.1d0 + rhor = 0.562d0 + bvcxl=2.d0 + bvcxr=2.d0 + bvcyl=0.d0 + bvcyr=4.710d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d0 /rhol + epsr = 3.d0 * 10.d0 /rhor + + ux=-2.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=-0.212d0 + uy=-0.590d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + +!!$ Switch-on Slow test of Komissarov 1999 -- compare at t=2. , n=150, x=[-1,1.5] + else if (CCTK_EQUALS(shock_case,"Komissarov4")) then + rhol = 1.78d-2 + rhor = 1.d-2 + bvcxl=1.d0 + bvcxr=1.d0 + bvcyl=1.022d0 + bvcyr=0.d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 0.1d0 /rhol + epsr = 3.d0 * 1.d0 /rhor + + ux=-0.765d0 + uy=-1.386d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=0.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + +!!$ Alfven wave test of Komissarov 1999 -- compare at t=2. , n=200, x=[-1,1.5] +!!$ Needs special setup -- FIX + else if (CCTK_EQUALS(shock_case,"Komissarov5")) then + rhol = 1.d0 + rhor = 1.d0 + bvcxl=3.d0 + bvcxr=3.d0 + bvcyl=3.d0 + bvcyr=-6.857d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d0 /rhol + epsr = 3.d0 * 1.d0 /rhor + + ux=0.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=3.70d0 + uy=5.76d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + +!!$ Compound wave test of Komissarov 1999 -- compare at t=0.1,0.75,1.5 , n=200, x=[-0.5,1.5] +!!$ Needs special setup -- FIX + else if (CCTK_EQUALS(shock_case,"Komissarov6")) then + rhol = 1.d0 + rhor = 1.d0 + bvcxl=3.d0 + bvcxr=3.d0 + bvcyl=3.d0 + bvcyr=-6.857d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d0 /rhol + epsr = 3.d0 * 1.d0 /rhor + + ux=0.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=3.70d0 + uy=5.76d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + +!!$ Shock Tube 1 test of Komissarov 1999 -- compare at t=1. , n=400, x=[-1,1.5] + else if (CCTK_EQUALS(shock_case,"Komissarov7")) then + rhol = 1.d0 + rhor = 1.d0 + bvcxl=1.d0 + bvcxr=1.d0 + bvcyl=0.d0 + bvcyr=0.d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d3 /rhol + epsr = 3.d0 * 1.d0 /rhor + + velxl = 0. + velyl = 0. + velzl = 0. + + velxr = 0. + velyr = 0. + velzr = 0. + +!!$ Shock Tube 2 test of Komissarov 1999 -- compare at t=1. , n=500, x=[-1.25,1.25] + else if (CCTK_EQUALS(shock_case,"Komissarov8")) then + rhol = 1.d0 + rhor = 1.d0 + bvcxl=0.d0 + bvcxr=0.d0 + bvcyl=20.d0 + bvcyr=0.d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 30.d0 /rhol + epsr = 3.d0 * 1.d0 /rhor + + velxl = 0. + velyl = 0. + velzl = 0. + + velxr = 0. + velyr = 0. + velzr = 0. + +!!$ Collision test of Komissarov 1999 -- compare at t=1.22 , n=200, x=[-1,1] + else if (CCTK_EQUALS(shock_case,"Komissarov9")) then + rhol = 1.d0 + rhor = 1.d0 + bvcxl=10.d0 + bvcxr=10.d0 + bvcyl=10.d0 + bvcyr=-10.d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d0 /rhol + epsr = 3.d0 * 1.d0 /rhor + + ux=5.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=-5.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + else call CCTK_WARN(0,"Shock case not recognized") end if if ( ((change_shock_direction==0).and.(direction .lt. 0.0d0)).or.& ((change_shock_direction==1).and.(direction .gt. 0.0d0)) ) then + +!!$ Left state + rho(i,j,k) = rhol velx(i,j,k) = velxl vely(i,j,k) = velyl @@ -217,6 +540,9 @@ Bvecy(i,j,k)=bvcyl Bvecz(i,j,k)=bvczl else + +!!$ Right state + rho(i,j,k) = rhor velx(i,j,k) = velxr vely(i,j,k) = velyr @@ -227,6 +553,58 @@ Bvecz(i,j,k)=bvczr end if + + if (CCTK_EQUALS(shocktube_type,"yshock")) then +!!$ Cycle x,y,z forward + tmp=velx(i,j,k) + velx(i,j,k)=velz(i,j,k) + velz(i,j,k)=vely(i,j,k) + vely(i,j,k)=tmp + tmp=Bvecx(i,j,k) + Bvecx(i,j,k)=Bvecz(i,j,k) + Bvecz(i,j,k)=Bvecy(i,j,k) + Bvecy(i,j,k)=tmp + else if (CCTK_EQUALS(shocktube_type,"zshock")) then +!!$ Cycle x,y,z backward + tmp=velx(i,j,k) + velx(i,j,k)=vely(i,j,k) + vely(i,j,k)=velz(i,j,k) + velz(i,j,k)=tmp + tmp=Bvecx(i,j,k) + Bvecx(i,j,k)=Bvecy(i,j,k) + Bvecy(i,j,k)=Bvecz(i,j,k) + Bvecz(i,j,k)=tmp + else if (CCTK_EQUALS(shocktube_type,"diagshock")) then +!!$ Rotated basis vectors necessary to evaluate the orthogonal matrix elements: +!!$ xhat = 1/sqrt(3)[1,1,1], yhat = 1/sqrt(2)[-1,1,0]; zhat = 1/sqrt(6)[-1,-1,2] +!!$ Orthogonal matrix constructed from the tensor product between the original +!!$ cartesian basis x's and new basis vectors xhat's, rotated towards the diagonal +!!$ shock normal and tangent directions. + tmp = OOSQRT3*velx(i,j,k) - OOSQRT2*vely(i,j,k) - OOSQRT6*velz(i,j,k) + tmp2 = OOSQRT3*velx(i,j,k) + OOSQRT2*vely(i,j,k) - OOSQRT6*velz(i,j,k) + tmp3 = OOSQRT3*velx(i,j,k) + 2.d0*OOSQRT6*velz(i,j,k) + velx(i,j,k)=tmp + vely(i,j,k)=tmp2 + velz(i,j,k)=tmp3 + tmp = OOSQRT3*Bvecx(i,j,k) - OOSQRT2*Bvecy(i,j,k) - OOSQRT6*Bvecz(i,j,k) + tmp2 = OOSQRT3*Bvecx(i,j,k) + OOSQRT2*Bvecy(i,j,k) - OOSQRT6*Bvecz(i,j,k) + tmp3 = OOSQRT3*Bvecx(i,j,k) + 2.d0*OOSQRT6*Bvecz(i,j,k) + Bvecx(i,j,k)=tmp + Bvecy(i,j,k)=tmp2 + Bvecz(i,j,k)=tmp3 + else if (CCTK_EQUALS(shocktube_type,"diagshock2d")) then +!!$ New basis: +!!$ xhat = 1/sqrt(2)[1,1,0], yhat = 1/sqrt(2)[-1,1,0]; zhat = [0,0,1] + tmp = OOSQRT2*velx(i,j,k) - OOSQRT2*vely(i,j,k) + tmp2 = OOSQRT2*velx(i,j,k) + OOSQRT2*vely(i,j,k) + velx(i,j,k)=tmp + vely(i,j,k)=tmp2 + tmp = OOSQRT2*Bvecx(i,j,k) - OOSQRT2*Bvecy(i,j,k) + tmp2 = OOSQRT2*Bvecx(i,j,k) + OOSQRT2*Bvecy(i,j,k) + Bvecx(i,j,k)=tmp + Bvecy(i,j,k)=tmp2 + endif + det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) if (CCTK_EQUALS(GRHydro_eos_type,"Polytype")) then @@ -257,3 +635,182 @@ return end subroutine GRHydro_shocktubeM + +subroutine GRHydro_Diagshock_BoundaryM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, nx, ny, nz, stenp1, minsum, maxsum, inew, jnew, knew + CCTK_INT :: xoff,yoff,zoff,indsum + CCTK_REAL :: det + + stenp1=GRHydro_stencil + 1 + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + minsum = 3*stenp1 + maxsum = nx+ny+nz-3*stenp1 + 3 + + xoff=0 + yoff=0 + zoff=0 + + do k=1,nz + if(k.lt.stenp1)then + zoff=k-stenp1 + else if(k.gt.nz-stenp1+1) then + zoff=k-(nz-stenp1+1) + else + zoff=0 + endif + + do j=1,ny + if(j.lt.stenp1)then + yoff=j-stenp1 + else if(j.gt.ny-stenp1+1) then + yoff=j-(ny-stenp1+1) + else + yoff=0 + endif + + do i=1,nx + if(i.lt.stenp1) then + xoff=i-stenp1 + else if(i.gt.nx-stenp1+1) then + xoff=i-(nx-stenp1+1) + else + xoff=0 + endif + + indsum = i+j+k + + if( (xoff.ne.0.or.yoff.ne.0.or.zoff.ne.0) .and. & + indsum.ge.minsum.and.indsum.le.maxsum) then +!!$ We can map the point to the interior diagonal, orthogonal to the shock. + + inew=indsum/3 + jnew=(indsum-inew)/2 + knew=indsum-inew-jnew + + dens(i,j,k) = dens(inew,jnew,knew) + sx(i,j,k) = sx(inew,jnew,knew) + sy(i,j,k) = sy(inew,jnew,knew) + sz(i,j,k) = sz(inew,jnew,knew) + tau(i,j,k) = tau(inew,jnew,knew) + Bvecx(i,j,k)=Bvecx(inew,jnew,knew) + Bvecy(i,j,k)=Bvecy(inew,jnew,knew) + Bvecz(i,j,k)=Bvecz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + endif + + enddo + enddo + enddo + +end subroutine GRHydro_Diagshock_BoundaryM + +subroutine GRHydro_Diagshock2D_BoundaryM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, kc, nx, ny, nz, stenp1, minsum, maxsum, inew, jnew, knew + CCTK_INT :: xoff,yoff,zoff,indsum + CCTK_REAL :: det + + stenp1=GRHydro_stencil + 1 + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + minsum = 2*stenp1 + maxsum = nx+ny-2*stenp1 + 2 + + xoff=0 + yoff=0 + zoff=0 + + do k=1,nz + if(k.lt.stenp1)then + zoff=k-stenp1 + else if(k.gt.nz-stenp1+1) then + zoff=k-(nz-stenp1+1) + else + zoff=0 + endif + + do j=1,ny + if(j.lt.stenp1)then + yoff=j-stenp1 + else if(j.gt.ny-stenp1+1) then + yoff=j-(ny-stenp1+1) + else + yoff=0 + endif + + do i=1,nx + if(i.lt.stenp1) then + xoff=i-stenp1 + else if(i.gt.nx-stenp1+1) then + xoff=i-(nx-stenp1+1) + else + xoff=0 + endif + + indsum = i+j + + if( (xoff.ne.0.or.yoff.ne.0) .and. & + indsum.ge.minsum.and.indsum.le.maxsum) then +!!$ We can map the point to the interior diagonal, orthogonal to the shock. + + inew=indsum/2 + jnew=indsum-inew + + dens(i,j,k) = dens(inew,jnew,k) + sx(i,j,k) = sx(inew,jnew,k) + sy(i,j,k) = sy(inew,jnew,k) + sz(i,j,k) = sz(inew,jnew,k) + tau(i,j,k) = tau(inew,jnew,k) + Bvecx(i,j,k)=Bvecx(inew,jnew,k) + Bvecy(i,j,k)=Bvecy(inew,jnew,k) + Bvecz(i,j,k)=Bvecz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + + else if( zoff.ne.0) then + + kc = nz/2 + + dens(i,j,k) = dens(i,j,kc) + sx(i,j,k) = sx(i,j,kc) + sy(i,j,k) = sy(i,j,kc) + sz(i,j,k) = sz(i,j,kc) + tau(i,j,k) = tau(i,j,kc) + Bvecx(i,j,k)=Bvecx(i,j,kc) + Bvecy(i,j,k)=Bvecy(i,j,kc) + Bvecz(i,j,k)=Bvecz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + endif + + enddo + enddo + enddo + +end subroutine GRHydro_Diagshock2D_BoundaryM File [modified]: make.code.defn Delta lines: +3 -1 =================================================================== --- trunk/src/make.code.defn 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/src/make.code.defn 2011-04-27 22:25:23 UTC (rev 123) @@ -14,7 +14,9 @@ CheckParam.c\ GRHydro_P2C2P.F90 \ GRHydro_P2C2PM.F90 \ - GRHydro_P2C2PM_polytype.F90 + GRHydro_P2C2PM_polytype.F90 \ + GRHydro_MonopoleM.F90 \ + GRHydro_CylindricalExplosionM.F90 # Subdirectories containing source files SUBDIRS = Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +18 -1 =================================================================== --- trunk/param.ccl 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/param.ccl 2011-04-27 22:25:23 UTC (rev 123) @@ -12,6 +12,8 @@ "only_atmo" :: "Set only a low atmosphere" "read_conformal":: "After reading in initial alp, rho and gxx from h5 files, sets the other quantities" "simple_wave" :: "Set initial data from Anile Miller Motta, Phys.Fluids. 26, 1450 (1983)" + "Monopole" :: "Monopole at the center" + "cylexp" :: "Cylindrical Explosion" } shares:ADMBase @@ -31,22 +33,35 @@ KEYWORD shocktube_type "Diagonal or parallel shock?" { "diagshock" :: "Diagonal across all axes" + "diagshock2d" :: "Diagonal across x-y axes" "xshock" :: "Parallel to x axis" "yshock" :: "Parallel to y axis" "zshock" :: "Parallel to z axis" "sphere" :: "spherically symmetric shock" -} "diagshock" +} "xshock" KEYWORD shock_case "Simple, Sod's problem or other?" { "Simple" :: "GRAstro_Hydro test case" "Sod" :: "Sod's problem" "Blast" :: "Strong blast wave" + "Balsaralike1" :: "Hydro version of Balsara Test #1" + "Balsara0" :: "Balsara Test #1, but unmagnetized" "Balsara1" :: "Balsara Test #1" "Balsara2" :: "Balsara Test #2" "Balsara3" :: "Balsara Test #3" "Balsara4" :: "Balsara Test #4" "Balsara5" :: "Balsara Test #5" + "Alfven" :: "Generical Alfven Test" + "Komissarov1" :: "Komissarov Test #1" + "Komissarov2" :: "Komissarov Test #2" + "Komissarov3" :: "Komissarov Test #3" + "Komissarov4" :: "Komissarov Test #4" + "Komissarov5" :: "Komissarov Test #5" + "Komissarov6" :: "Komissarov Test #6" + "Komissarov7" :: "Komissarov Test #7" + "Komissarov8" :: "Komissarov Test #8" + "Komissarov9" :: "Komissarov Test #9" } "Sod" REAL shock_xpos "Position of shock plane: x" @@ -124,3 +139,5 @@ USES REAL initial_rho_abs_min USES REAL initial_rho_rel_min USES REAL initial_atmosphere_factor +USES int GRHydro_stencil +USES BOOLEAN clean_divergence File [modified]: schedule.ccl Delta lines: +30 -0 =================================================================== --- trunk/schedule.ccl 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/schedule.ccl 2011-04-27 22:25:23 UTC (rev 123) @@ -5,6 +5,13 @@ LANG: C } "Check parameters" +if (CCTK_Equals(initial_hydro,"Monopole")) { + schedule GRHydro_MonopoleM in HydroBase_Initial + { + LANG: Fortran + } "Monopole initial data - MHD version" +} + if (CCTK_Equals(initial_hydro,"shocktube")) { if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { @@ -13,6 +20,22 @@ LANG: Fortran } "Shocktube initial data - MHD version" + if(CCTK_Equals(shocktube_type,"diagshock")) + { + schedule GRHydro_Diagshock_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection + { + LANG: Fortran + } "Diagonal shock boundary conditions" + } + + if(CCTK_Equals(shocktube_type,"diagshock2d")) + { + schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection + { + LANG: Fortran + } "2-D Diagonal shock boundary conditions" + } + } else { schedule GRHydro_shocktube in HydroBase_Initial { @@ -22,6 +45,13 @@ } } +if (CCTK_Equals(initial_hydro,"cylexp")) { + schedule GRHydro_CylindricalExplosionM in HydroBase_Initial + { + LANG: Fortran + } "Cylindrical Explosion initial data - MHD-only" +} + if (CCTK_Equals(initial_data,"con2primtest")) { STORAGE:GRHydro_init_data_reflevel schedule GRHydro_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE GRHydro_con2primtest

1 0

[svn:einsteintoolkit] GRHydro/branches/divclean/src/ (Rev. 233)
by tanja.bode＠physics.gatech.edu 27 Apr '11

27 Apr '11

User: tbode Date: 2011/04/27 05:20 PM Modified: /branches/divclean/src/ GRHydro_Con2Prim.F90, GRHydro_Con2PrimM.F90 Log: Fix wrongly overwritten Con2Prim. File Changes: Directory: /branches/divclean/src/ ================================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM.F90 Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] GRHydro/branches/divclean/ (Rev. 232)
by tanja.bode＠physics.gatech.edu 27 Apr '11

27 Apr '11

User: tbode Date: 2011/04/27 04:02 PM Modified: /branches/divclean/ interface.ccl, schedule.ccl /branches/divclean/src/ GRHydro_Boundaries.F90, GRHydro_CalcUpdate.F90, GRHydro_Con2Prim.F90, GRHydro_HLLEM.F90, GRHydro_Prim2ConM.F90, GRHydro_Reconstruct.F90, GRHydro_ReconstructPoly.F90, GRHydro_RegisterVars.cc, GRHydro_UpdateMaskM.F90 Log: Add Josh's changes. File Changes: Directory: /branches/divclean/ ============================== File [modified]: interface.ccl Delta lines: <Unavailable> File [modified]: schedule.ccl Delta lines: <Unavailable> Directory: /branches/divclean/src/ ================================== File [modified]: GRHydro_Boundaries.F90 Delta lines: <Unavailable> File [modified]: GRHydro_CalcUpdate.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2Prim.F90 Delta lines: <Unavailable> File [modified]: GRHydro_HLLEM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Prim2ConM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Reconstruct.F90 Delta lines: <Unavailable> File [modified]: GRHydro_ReconstructPoly.F90 Delta lines: <Unavailable> File [modified]: GRHydro_RegisterVars.cc Delta lines: <Unavailable> File [modified]: GRHydro_UpdateMaskM.F90 Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] GRHydro/branches/divclean/src/ (Rev. 231)
by tanja.bode＠physics.gatech.edu 27 Apr '11

27 Apr '11

User: tbode Date: 2011/04/27 03:30 PM Modified: /branches/divclean/src/ GRHydro_SourceM.F90 Log: Modify rhs for divergence cleaning. File Changes: Directory: /branches/divclean/src/ ================================== File [modified]: GRHydro_SourceM.F90 Delta lines: +50 -3 =================================================================== --- branches/divclean/src/GRHydro_SourceM.F90 2011-04-27 17:51:47 UTC (rev 230) +++ branches/divclean/src/GRHydro_SourceM.F90 2011-04-27 20:30:04 UTC (rev 231) @@ -33,6 +33,12 @@ #define Bvecx(i,j,k) Bvec(i,j,k,1) #define Bvecy(i,j,k) Bvec(i,j,k,2) #define Bvecz(i,j,k) Bvec(i,j,k,3) +#define Bconsx(i,j,k) Bcons(i,j,k,1) +#define Bconsy(i,j,k) Bcons(i,j,k,2) +#define Bconsz(i,j,k) Bcons(i,j,k,3) +#define Bconsrhsx(i,j,k) Bconsrhs(i,j,k,1) +#define Bconsrhsy(i,j,k) Bconsrhs(i,j,k,2) +#define Bconsrhsz(i,j,k) Bconsrhs(i,j,k,3) /*@@ @routine SourceTermsM @@ -76,8 +82,9 @@ CCTK_REAL :: halfshiftdgx, halfshiftdgy, halfshiftdgz CCTK_REAL :: halfTdgx, halfTdgy, halfTdgz CCTK_REAL :: invalp, invalp2 - CCTK_REAL :: dx_psidc, dy_psidc, dz_psidc + !!CCTK_REAL :: dx_psidc, dy_psidc, dz_psidc CCTK_REAL :: bvcx_source, bvcy_source, bvcz_source + CCTK_REAL :: dx_det_bydet, dy_det_bydet, dz_det_bydet CCTK_REAL :: Bvecxlow,Bvecylow,Bveczlow,bdotv,b2,dum,bxlow,bylow,bzlow CCTK_REAL :: bt,bx,by,bz,rhohstarW2,pstar @@ -103,7 +110,7 @@ densrhs = 0.d0 srhs = 0.d0 taurhs = 0.d0 - Bvecrhs=0.d0 + Bconsrhs=0.d0 if (evolve_tracer .ne. 0) then cons_tracerrhs = 0.d0 @@ -436,7 +443,47 @@ srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sz_source taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source if(clean_divergence.ne.0) then - psidcrhs(i,j,k) = -1.d0*ch_dc*ch_dc/cp_dc/cp_dc*psidc(i,j,k) + psidcrhs(i,j,k) = -1.d0 * ( ch_dc*ch_dc/cp_dc/cp_dc*alp(i,j,k) + & + dx_shiftx + dy_shifty + dz_shiftz ) * psidc(i,j,k) + & + Bconsx(i,j,k) * ( dx_alp - half*alp(i,j,k) * & + ( uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + 2.d0*uxy*dx_gxy + & + 2.d0*uxz*dx_gxz + 2.d0*uyz*dx_gyz ) )/ sqrtdet + & + Bconsy(i,j,k) * (dy_alp - half*alp(i,j,k) * & + ( uxx*dy_gxx + uyy*dy_gyy + uzz*dy_gzz + 2.d0*uxy*dy_gxy + & + 2.d0*uxz*dy_gxz + 2.d0*uyz*dy_gyz ) )/ sqrtdet + & + Bconsz(i,j,k) * (dz_alp - half*alp(i,j,k) * & + ( uxx*dz_gxx + uyy*dy_gyy + uzz*dy_gzz + 2.d0*uxy*dy_gxy + & + 2.d0*uxz*dz_gxz + 2.d0*uyz*dy_gyz ) )/ sqrtdet + + !!$ g^{jk} d_i g_{kj} = d_i (g) / det + dx_det_bydet = uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + & + 2.d0*(uxy*dx_gxy+uxz*dx_gxz+uyz*dx_gyz) + dy_det_bydet = uxx*dy_gxx + uyy*dy_gyy + uzz*dy_gzz + & + 2.d0*(uxy*dy_gxy+uxz*dy_gxz+uyz*dy_gyz) + dz_det_bydet = uxx*dz_gxx + uyy*dz_gyy + uzz*dz_gzz + & + 2.d0*(uxy*dz_gxy+uxz*dz_gxz+uyz*dz_gyz) + + bvcx_source = -1.d0 * ( Bconsx(i,j,k)*dx_shiftx + & + Bconsy(i,j,k)*dy_shiftx + Bconsz(i,j,k)*dz_shiftx ) + & + psidc(i,j,k)*sqrtdet*( uxx*dx_alp+uxy*dy_alp+uxz*dz_alp ) + & + psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxx*dx_det_bydet + & + uxy*dy_det_bydet + uxz*dz_det_bydet ) + + bvcy_source = -1.d0 * ( Bconsx(i,j,k)*dx_shifty + & + Bconsy(i,j,k)*dy_shifty + Bconsz(i,j,k)*dz_shifty ) + & + psidc(i,j,k)*sqrtdet*( uxy*dx_alp+uyy*dy_alp+uyz*dz_alp ) + & + psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxy*dx_det_bydet + & + uyy*dy_det_bydet + uyz*dz_det_bydet ) + + bvcz_source = -1.d0 * ( Bconsx(i,j,k)*dx_shiftz + & + Bconsy(i,j,k)*dy_shiftz + Bconsz(i,j,k)*dz_shiftz ) + & + psidc(i,j,k)*sqrtdet*( uxz*dx_alp+uyz*dy_alp+uzz*dz_alp ) + & + psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxz*dx_det_bydet + & + uyz*dy_det_bydet + uzz*dz_det_bydet ) + + Bconsrhsx(i,j,k) = bvcx_source + Bconsrhsy(i,j,k) = bvcy_source + Bconsrhsz(i,j,k) = bvcz_source endif enddo

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 71)
by cott＠tapir.caltech.edu 27 Apr '11

27 Apr '11

User: cott Date: 2011/04/27 01:57 PM Modified: / ET.tex /local_bibtex/ ott_references.bib Log: * add a ref File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +1 -1 =================================================================== --- ET.tex 2011-04-27 18:54:54 UTC (rev 70) +++ ET.tex 2011-04-27 18:57:16 UTC (rev 71) @@ -965,7 +965,7 @@ specific internal energy and the pressure, respectively. \subsubsection{Spacetime Curvature Evolution} The Einstein Toolkit -curvature evolution code \codename{McLachlan}~\cite{brown:09} is +curvature evolution code \codename{McLachlan}~\cite{brown:09,reisswig:11ccwave} is auto-generated from tensor equations via \codename{Kranc} (Sec.~\ref{sec:kranc}). It implements the Einstein equations in a $3+1$ split as a Cauchy initial boundary value problem Directory: /local_bibtex/ ========================= File [modified]: ott_references.bib Delta lines: +12 -0 =================================================================== --- local_bibtex/ott_references.bib 2011-04-27 18:54:54 UTC (rev 70) +++ local_bibtex/ott_references.bib 2011-04-27 18:57:16 UTC (rev 71) @@ -165,3 +165,15 @@ fulltexturl = {http://link.aps.org/abstract/PRD/v79/e044023}, doi = {10.1103/PhysRevD.79.044023}, } + +@ARTICLE{reisswig:11ccwave, + author = {{Reisswig}, C. and {Ott}, C.~D. and {Sperhake}, U. and {Schnetter}, E.}, + title = "{Gravitational wave extraction in simulations of rotating stellar core collapse}", + journal = {\prd}, + keywords = {Numerical relativity, Finite-difference methods, Wave generation and sources, Supernovae}, + year = 2011, + month = mar, + volume = 83, + number = 6, + pages = {064008}, +}

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits