Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/examples/tov/ (Rev. 80)
by knarf＠cct.lsu.edu 28 Apr '11

28 Apr '11

User: knarf Date: 2011/04/28 01:46 PM Modified: /examples/tov/ tov.rpar Log: fix parfile File Changes: Directory: /examples/tov/ ========================= File [modified]: tov.rpar Delta lines: +0 -1 =================================================================== --- examples/tov/tov.rpar 2011-04-28 18:26:13 UTC (rev 79) +++ examples/tov/tov.rpar 2011-04-28 18:45:59 UTC (rev 80) @@ -216,7 +216,6 @@ ML_BSSN::ML_Ham ML_BSSN::ML_mom GRHydro::dens - GRHydro::eps " IOASCII::out1D_every = $out_every

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/examples/tov/ (Rev. 79)
by knarf＠cct.lsu.edu 28 Apr '11

28 Apr '11

User: knarf Date: 2011/04/28 01:26 PM Modified: /examples/tov/ tov.rpar Log: fix naming, add a bit more output File Changes: Directory: /examples/tov/ ========================= File [modified]: tov.rpar Delta lines: +5 -3 =================================================================== --- examples/tov/tov.rpar 2011-04-28 18:22:51 UTC (rev 78) +++ examples/tov/tov.rpar 2011-04-28 18:26:13 UTC (rev 79) @@ -1,13 +1,13 @@ #!/usr/bin/perl -w use strict; -foreach my $type (("", "pert")) { +foreach my $type (("", "_pert")) { foreach my $res ((16, 8, 4, 2)) { my $out_every = 256/$res; my $K_TOV = 100; my $K_change = ""; - if ($type eq "pert") { + if ($type eq "_pert") { $K_TOV = 99; $K_change = "TOVSolver::TOV_Transform_to_K = 100"; } @@ -215,6 +215,8 @@ ADMBase::curv ML_BSSN::ML_Ham ML_BSSN::ML_mom + GRHydro::dens + GRHydro::eps " IOASCII::out1D_every = $out_every @@ -235,7 +237,7 @@ " PAR_EOF - open(OUT, ">tov_${res}_${type}.par"); + open(OUT, ">tov_${res}${type}.par"); print OUT "$par"; close(OUT); }}

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 78)
by schnetter＠cct.lsu.edu 28 Apr '11

28 Apr '11

User: eschnett Date: 2011/04/28 01:22 PM Modified: / ET.tex Log: Typo File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +1 -1 =================================================================== --- ET.tex 2011-04-28 18:22:37 UTC (rev 77) +++ ET.tex 2011-04-28 18:22:51 UTC (rev 78) @@ -925,7 +925,7 @@ the ODE solution nor the resulting density profile. -\subsection{Equations of States}\pages{1 Christian} +\subsection{Equations of State}\pages{1 Christian} \subsection{Spacetime Curvature and Hydrodynamics Evolution} \todo{Christian in charge}

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/local_bibtex/ (Rev. 77)
by schnetter＠cct.lsu.edu 28 Apr '11

28 Apr '11

User: eschnett Date: 2011/04/28 01:22 PM Modified: /local_bibtex/ ott_references.bib, references.bib Log: Expand \prd File Changes: Directory: /local_bibtex/ ========================= File [modified]: ott_references.bib Delta lines: +1 -1 =================================================================== --- local_bibtex/ott_references.bib 2011-04-28 17:47:46 UTC (rev 76) +++ local_bibtex/ott_references.bib 2011-04-28 18:22:37 UTC (rev 77) @@ -169,7 +169,7 @@ @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}, + journal = {Phys. Rev. D}, keywords = {Numerical relativity, Finite-difference methods, Wave generation and sources, Supernovae}, year = 2011, month = mar, File [modified]: references.bib Delta lines: +2 -2 =================================================================== --- local_bibtex/references.bib 2011-04-28 17:47:46 UTC (rev 76) +++ local_bibtex/references.bib 2011-04-28 18:22:37 UTC (rev 77) @@ -5519,7 +5519,7 @@ {Pfeiffer}, H.~P. and {Scheel}, M.~A. and {Cook}, G.~B. and {Teukolsky}, S.~A.}, title = "{High-accuracy comparison of numerical relativity simulations with post-Newtonian expansions}", - journal = {\prd}, + journal = {Phys. Rev. D}, eprint = {arXiv:0710.0158}, year = 2007, month = dec, @@ -29270,7 +29270,7 @@ @ARTICLE{Nollert1993b}, author = {{Nollert}, H.-P.}, title = "{Quasinormal modes of Schwarzschild black holes: The determination of quasinormal frequencies with very large imaginary parts}", - journal = {\prd}, + journal = {Phys. Rev. D}, year = 1993, month = jun, volume = 47,

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 244)
by tanja.bode＠physics.gatech.edu 28 Apr '11

28 Apr '11

User: tbode Date: 2011/04/28 12:55 PM Modified: /trunk/ interface.ccl, schedule.ccl /trunk/src/ GRHydro_Boundaries.F90, GRHydro_CalcUpdate.F90, GRHydro_Con2PrimM.F90, GRHydro_Con2PrimM_polytype_pt.c, GRHydro_Con2PrimM_pt.c, GRHydro_FluxM.F90, GRHydro_HLLEM.F90, GRHydro_InterfacesM.h, GRHydro_Prim2ConM.F90, GRHydro_Reconstruct.F90, GRHydro_ReconstructPoly.F90, GRHydro_RegisterVars.cc, GRHydro_SourceM.F90, GRHydro_UpdateMaskM.F90 Log: MERGE divergence cleaning feature into trunk. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Boundaries.F90 Delta lines: <Unavailable> File [modified]: GRHydro_CalcUpdate.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM_polytype_pt.c Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM_pt.c Delta lines: <Unavailable> File [modified]: GRHydro_FluxM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_HLLEM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_InterfacesM.h 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_SourceM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_UpdateMaskM.F90 Delta lines: <Unavailable> Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: <Unavailable> File [modified]: schedule.ccl Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/examples/tov/ (Rev. 76)
by knarf＠cct.lsu.edu 28 Apr '11

28 Apr '11

User: knarf Date: 2011/04/28 12:47 PM Modified: /examples/tov/ tov.rpar Log: add perturbed tov parameters (not in ET yet File Changes: Directory: /examples/tov/ ========================= File [modified]: tov.rpar Delta lines: +13 -3 =================================================================== --- examples/tov/tov.rpar 2011-04-28 17:23:25 UTC (rev 75) +++ examples/tov/tov.rpar 2011-04-28 17:47:46 UTC (rev 76) @@ -1,9 +1,18 @@ #!/usr/bin/perl -w use strict; +foreach my $type (("", "pert")) { foreach my $res ((16, 8, 4, 2)) { my $out_every = 256/$res; + my $K_TOV = 100; + my $K_change = ""; + if ($type eq "pert") { + $K_TOV = 99; + $K_change = "TOVSolver::TOV_Transform_to_K = 100"; + } + + my $par = <<PAR_EOF; ## mclachlan tov_static ActiveThorns = "Time MoL" @@ -184,7 +193,8 @@ TOVSolver::TOV_Rho_Central[0] = 1.28e-3 TOVSolver::TOV_Gamma[0] = 2.0 -TOVSolver::TOV_K[0] = 100.0 +TOVSolver::TOV_K[0] = $K_TOV +$K_change # I/O @@ -225,8 +235,8 @@ " PAR_EOF - open(OUT, ">tov_${res}.par"); + open(OUT, ">tov_${res}_${type}.par"); print OUT "$par"; close(OUT); -} +}}

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/examples/tov/ (Rev. 75)
by knarf＠cct.lsu.edu 28 Apr '11

28 Apr '11

User: knarf Date: 2011/04/28 12:23 PM Modified: /examples/tov/ tov.rpar Log: output simulation speed File Changes: Directory: /examples/tov/ ========================= File [modified]: tov.rpar Delta lines: +1 -1 =================================================================== --- examples/tov/tov.rpar 2011-04-28 17:22:12 UTC (rev 74) +++ examples/tov/tov.rpar 2011-04-28 17:23:25 UTC (rev 75) @@ -191,7 +191,7 @@ IO::out_dir = \$parfile IOBasic::outInfo_every = $out_every -IOBasic::outInfo_vars = "HydroBase::rho ML_BSSN::ML_Ham" +IOBasic::outInfo_vars = "Carpet::physical_time_per_hour HydroBase::rho ML_BSSN::ML_Ham" IOScalar::outScalar_every = $out_every IOScalar::one_file_per_group = yes

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/examples/tov/ (Rev. 74)
by knarf＠cct.lsu.edu 28 Apr '11

28 Apr '11

User: knarf Date: 2011/04/28 12:22 PM Modified: /examples/tov/ tov.rpar Log: small updates File Changes: Directory: /examples/tov/ ========================= File [modified]: tov.rpar Delta lines: +10 -16 =================================================================== --- examples/tov/tov.rpar 2011-04-28 16:21:09 UTC (rev 73) +++ examples/tov/tov.rpar 2011-04-28 17:22:12 UTC (rev 74) @@ -1,7 +1,7 @@ #!/usr/bin/perl -w use strict; -foreach my $res ((8, 4, 2)) { +foreach my $res ((16, 8, 4, 2)) { my $out_every = 256/$res; my $par = <<PAR_EOF; @@ -15,8 +15,6 @@ ActiveThorns = "Carpet CarpetLib CarpetReduce CarpetRegrid2 CarpetInterp" ActiveThorns = "CarpetIOASCII CarpetIOScalar CarpetIOHDF5 CarpetIOBasic" -ActiveThorns = "ADMConstraints" - # grid parameters CartGrid3D::type = "coordbase" CartGrid3D::domain = "full" @@ -24,9 +22,9 @@ CoordBase::xmin = 0.0 CoordBase::ymin = 0.0 CoordBase::zmin = 0.0 -CoordBase::xmax = 240.0 -CoordBase::ymax = 240.0 -CoordBase::zmax = 240.0 +CoordBase::xmax = 320.0 +CoordBase::ymax = 320.0 +CoordBase::zmax = 320.0 CoordBase::dx = $res CoordBase::dy = $res CoordBase::dz = $res @@ -69,14 +67,10 @@ ADMBase::metric_type = "physical" -ADMConstraints::bound = "static" -ADMConstraints::constraints_timelevels = 3 -ADMConstraints::constraints_persist = yes - SpaceMask::use_mask = "yes" Cactus::terminate = "time" -Cactus::cctk_final_time = 1000 +Cactus::cctk_final_time = 10000 Carpet::domain_from_coordbase = "yes" Carpet::enable_all_storage = no @@ -197,7 +191,7 @@ IO::out_dir = \$parfile IOBasic::outInfo_every = $out_every -IOBasic::outInfo_vars = "HydroBase::rho ADMBase::lapse" +IOBasic::outInfo_vars = "HydroBase::rho ML_BSSN::ML_Ham" IOScalar::outScalar_every = $out_every IOScalar::one_file_per_group = yes @@ -209,8 +203,8 @@ ADMBase::lapse ADMBase::metric ADMBase::curv - ADMConstraints::ham - ADMConstraints::momentum + ML_BSSN::ML_Ham + ML_BSSN::ML_mom " IOASCII::out1D_every = $out_every @@ -226,8 +220,8 @@ ADMBase::lapse ADMBase::metric ADMBase::curv - ADMConstraints::ham - ADMConstraints::momentum + ML_BSSN::ML_Ham + ML_BSSN::ML_mom " PAR_EOF

1 0

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

28 Apr '11

User: tbode Date: 2011/04/28 11:47 AM Modified: /branches/divclean/src/ GRHydro_CalcUpdate.F90, GRHydro_FluxM.F90, GRHydro_HLLEM.F90, GRHydro_Prim2ConM.F90, GRHydro_Reconstruct.F90, GRHydro_SourceM.F90 Log: Josh's fixes. File Changes: Directory: /branches/divclean/src/ ================================== File [modified]: GRHydro_CalcUpdate.F90 Delta lines: +7 -0 =================================================================== --- branches/divclean/src/GRHydro_CalcUpdate.F90 2011-04-28 01:28:55 UTC (rev 242) +++ branches/divclean/src/GRHydro_CalcUpdate.F90 2011-04-28 16:47:14 UTC (rev 243) @@ -102,6 +102,13 @@ Bvec(i+xoffset,j+yoffset,k+zoffset,flux_direction)) divB(i,j,k) = divB(i,j,k) + ( alp_l * Bvec_l - alp_r * Bvec_r ) * idx endif + +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'update:', & +!!$ Bconsrhs(i,j,k,1),Bconsxflux(i,j,k), & +!!$ Bconsxflux(i-xoffset,j-yoffset,k-zoffset), & +!!$ Bconsrhs(i,j,k,3),Bconszflux(i,j,k), & +!!$ Bconszflux(i-xoffset,j-yoffset,k-zoffset) + endif if (evolve_tracer .ne. 0) then File [modified]: GRHydro_FluxM.F90 Delta lines: +7 -10 =================================================================== --- branches/divclean/src/GRHydro_FluxM.F90 2011-04-28 01:28:55 UTC (rev 242) +++ branches/divclean/src/GRHydro_FluxM.F90 2011-04-28 16:47:14 UTC (rev 243) @@ -23,13 +23,10 @@ CCTK_REAL :: vxt,vyt,vzt,bsubx,bsuby,bsubz,ab0,w CCTK_REAL :: det,alp,beta,pressstar CCTK_REAL :: velm - CCTK_REAL :: sdet,psipstar,psiBx,psiBy,psiBz + CCTK_REAL :: sdet,psipstar sdet=sqrt(det) psipstar=pressstar*sdet - psiBx=Bx*sdet - psiBy=By*sdet - psiBz=Bz*sdet !!$ We actually need all three values of vtilde = v^i - beta^i/alp, as well as velm=vxt+beta/alp @@ -38,19 +35,19 @@ !!$ In the notation of Anton et al.: [psi^6 D] vtilde^i densf = dens * vxt - sxf = sx*vxt+psipstar-bsubx*psiBx/w + sxf = sx*vxt+psipstar-bsubx*Bx/w - syf = sy*vxt-bsuby*psiBx/w + syf = sy*vxt-bsuby*Bx/w - szf = sz*vxt-bsubz*psiBx/w + szf = sz*vxt-bsubz*Bx/w !!$ [psi^6 tau] vtilde^i +p* v^i - alp b^0 B^i/w - tauf = tau*vxt + psipstar*velm - ab0*psiBx/w + tauf = tau*vxt + psipstar*velm - ab0*Bx/w !!$ [psi^6 (B^k vtilde^i - B^i vtilde^k)] bxf = 0.0 - byf = psiBy * vxt - psiBx*vyt - bzf = psiBz * vxt - psiBx*vzt + byf = By * vxt - Bx*vyt + bzf = Bz * vxt - Bx*vzt end subroutine num_x_fluxM File [modified]: GRHydro_HLLEM.F90 Delta lines: +88 -78 =================================================================== --- branches/divclean/src/GRHydro_HLLEM.F90 2011-04-28 01:28:55 UTC (rev 242) +++ branches/divclean/src/GRHydro_HLLEM.F90 2011-04-28 16:47:14 UTC (rev 243) @@ -40,10 +40,9 @@ DECLARE_CCTK_PARAMETERS integer :: i, j, k, m - CCTK_REAL, dimension(8) :: fplus,fminus,f1,qdiff - CCTK_REAL, dimension(7) :: prim_p, prim_m - CCTK_REAL, dimension(5) :: lamminus,lamplus,cons_p,cons_m - CCTK_REAL, dimension(3) :: mag_p,mag_m + CCTK_REAL, dimension(8) :: cons_p,cons_m,fplus,fminus,f1,qdiff + CCTK_REAL, dimension(10) :: prim_p, prim_m + CCTK_REAL, dimension(5) :: lamminus,lamplus CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det, sdet CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & @@ -83,8 +82,6 @@ lamplus = 0.d0 cons_p = 0.d0 cons_m = 0.d0 - mag_p=0.d0 - mag_m=0.d0 fplus = 0.d0 fminus = 0.d0 qdiff = 0.d0 @@ -101,20 +98,20 @@ cons_p(3) = syplus(i,j,k) cons_p(4) = szplus(i,j,k) cons_p(5) = tauplus(i,j,k) - - mag_p(1) = Bvecxplus(i,j,k) - mag_p(2) = Bvecyplus(i,j,k) - mag_p(3) = Bveczplus(i,j,k) + cons_p(6) = Bconsxplus(i,j,k) + cons_p(7) = Bconsyplus(i,j,k) + cons_p(8) = Bconszplus(i,j,k) cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset) cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset) cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset) cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset) cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(6) = Bconsxminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(7) = Bconsyminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(8) = Bconszminus(i+xoffset,j+yoffset,k+zoffset) - mag_m(1) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset) - mag_m(2) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset) - mag_m(3) = Bveczminus(i+xoffset,j+yoffset,k+zoffset) +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'HLLEM0:',cons_m(6),cons_p(6),cons_m(7),cons_p(7),cons_m(8),cons_p(8) prim_p(1) = rhoplus(i,j,k) prim_p(2) = velxplus(i,j,k) @@ -123,6 +120,9 @@ prim_p(5) = epsplus(i,j,k) prim_p(6) = pressplus(i,j,k) prim_p(7) = w_lorentzplus(i,j,k) + prim_p(8) = Bvecxplus(i,j,k) + prim_p(9) = Bvecyplus(i,j,k) + prim_p(10) = Bveczplus(i,j,k) prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset) prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset) @@ -131,6 +131,9 @@ prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) prim_m(7) = w_lorentzminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(8) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(9) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(10)= Bveczminus(i+xoffset,j+yoffset,k+zoffset) if(clean_divergence.ne.0) then psidcp = psidcplus(i,j,k) @@ -191,11 +194,11 @@ vztm = prim_m(4)-avg_betaz/avg_alp call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & - prim_p(2),prim_p(3),prim_p(4),mag_p(1),mag_p(2),mag_p(3), & + prim_p(2),prim_p(3),prim_p(4),prim_p(8),prim_p(9),prim_p(10), & velxlowp,velylowp,velzlowp,Bvecxlowp,Bvecylowp,Bveczlowp, & bdotvp,b2p,v2p,wp,bxlowp,bylowp,bzlowp) call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & - prim_m(2),prim_m(3),prim_m(4),mag_m(1),mag_m(2),mag_m(3), & + prim_m(2),prim_m(3),prim_m(4),prim_m(8),prim_m(9),prim_m(10), & velxlowm,velylowm,velzlowm,Bvecxlowm,Bvecylowm,Bveczlowm, & bdotvm,b2m,v2m,wm,bxlowm,bylowm,bzlowm) @@ -224,41 +227,41 @@ if (flux_direction == 1) then call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& - mag_p(1),mag_p(2),mag_p(3),& + cons_p(6),cons_p(7),cons_p(8),& f1(1),f1(2),f1(3),f1(4),f1(5),f1(6),f1(7),f1(8),& vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6) + avg_alp*sdet*uxxh*psidcp - sdet*mag_p(1)*avg_betax - f1(7)=f1(7) + avg_alp*sdet*uxyh*psidcp - sdet*mag_p(1)*avg_betay - f1(8)=f1(8) + avg_alp*sdet*uxzh*psidcp - sdet*mag_p(1)*avg_betaz - psidcf = avg_alp*mag_p(1)-psidcp*avg_betax + f1(6)=f1(6) + avg_alp*sdet*uxxh*psidcp - cons_p(6)*avg_betax + f1(7)=f1(7) + avg_alp*sdet*uxyh*psidcp - cons_p(6)*avg_betay + f1(8)=f1(8) + avg_alp*sdet*uxzh*psidcp - cons_p(6)*avg_betaz + psidcf = avg_alp*cons_p(6)/sdet-psidcp*avg_betax endif else if (flux_direction == 2) then call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& - mag_p(2),mag_p(3),mag_p(1),& + cons_p(7),cons_p(8),cons_p(6),& f1(1),f1(3),f1(4),f1(2),f1(5),f1(7),f1(8),f1(6),& vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6) + avg_alp*sdet*uxyh*psidcp - sdet*mag_p(2)*avg_betax - f1(7)=f1(7) + avg_alp*sdet*uyyh*psidcp - sdet*mag_p(2)*avg_betay - f1(8)=f1(8) + avg_alp*sdet*uyzh*psidcp - sdet*mag_p(2)*avg_betaz - psidcf = avg_alp*mag_p(2)-psidcp*avg_betay + f1(6)=f1(6) + avg_alp*sdet*uxyh*psidcp - cons_p(7)*avg_betax + f1(7)=f1(7) + avg_alp*sdet*uyyh*psidcp - cons_p(7)*avg_betay + f1(8)=f1(8) + avg_alp*sdet*uyzh*psidcp - cons_p(7)*avg_betaz + psidcf = avg_alp*cons_p(7)/sdet-psidcp*avg_betay endif else if (flux_direction == 3) then call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& - mag_p(3),mag_p(1),mag_p(2),& + cons_p(8),cons_p(6),cons_p(7),& f1(1),f1(4),f1(2),f1(3),f1(5),f1(8),f1(6),f1(7), & vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - f1(6)=f1(6) + avg_alp*sdet*uxzh*psidcp - sdet*mag_p(3)*avg_betax - f1(7)=f1(7) + avg_alp*sdet*uyzh*psidcp - sdet*mag_p(3)*avg_betay - f1(8)=f1(8) + avg_alp*sdet*uzzh*psidcp - sdet*mag_p(3)*avg_betaz - psidcf = avg_alp*mag_p(3)-psidcp*avg_betaz + f1(6)=f1(6) + avg_alp*sdet*uxzh*psidcp - cons_p(8)*avg_betax + f1(7)=f1(7) + avg_alp*sdet*uyzh*psidcp - cons_p(8)*avg_betay + f1(8)=f1(8) + avg_alp*sdet*uzzh*psidcp - cons_p(8)*avg_betaz + psidcf = avg_alp*cons_p(8)/sdet-psidcp*avg_betaz endif else @@ -284,11 +287,10 @@ qdiff(3) = cons_m(3) - cons_p(3) qdiff(4) = cons_m(4) - cons_p(4) qdiff(5) = cons_m(5) - cons_p(5) + qdiff(6) = cons_m(6) - cons_p(6) + qdiff(7) = cons_m(7) - cons_p(7) + qdiff(8) = cons_m(8) - cons_p(8) - qdiff(6) = mag_m(1) - mag_p(1) - qdiff(7) = mag_m(2) - mag_p(2) - qdiff(8) = mag_m(3) - mag_p(3) - if (clean_divergence.ne.0) then psidcdiff = psidcm - psidcp endif @@ -298,107 +300,108 @@ if (flux_direction == 1) then call eigenvaluesM(GRHydro_eos_handle,& prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & - mag_m(1),mag_m(2),mag_m(3),& + prim_m(8),prim_m(9),prim_m(10),& lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & - mag_p(1),mag_p(2),mag_p(3),& + prim_p(8),prim_p(9),prim_p(10),& lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& - mag_p(1),mag_p(2),mag_p(3),& + cons_p(6),cons_p(7),cons_p(8),& fplus(1),fplus(2),fplus(3),fplus(4),fplus(5),fplus(6),fplus(7),fplus(8),& vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) call num_x_fluxM(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),& - mag_m(1),mag_m(2),mag_m(3),& + cons_m(6),cons_m(7),cons_m(8),& fminus(1),fminus(2),fminus(3),fminus(4),fminus(5),& fminus(6),fminus(7),fminus(8),& vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + avg_alp*sdet*uxxh*psidcm - sdet*mag_m(1)*avg_betax - fminus(7)=fminus(7) + avg_alp*sdet*uxyh*psidcm - sdet*mag_m(1)*avg_betay - fminus(8)=fminus(8) + avg_alp*sdet*uxzh*psidcm - sdet*mag_m(1)*avg_betaz - fplus(6)=fplus(6) + avg_alp*sdet*uxxh*psidcp - sdet*mag_p(1)*avg_betax - fplus(7)=fplus(7) + avg_alp*sdet*uxyh*psidcp - sdet*mag_p(1)*avg_betay - fplus(8)=fplus(8) + avg_alp*sdet*uxzh*psidcp - sdet*mag_p(1)*avg_betaz - psidcfp = avg_alp*mag_p(1)-avg_betax*psidcp - psidcfm = avg_alp*mag_m(1)-avg_betax*psidcm + fminus(6)=fminus(6) + avg_alp*sdet*uxxh*psidcm - cons_m(6)*avg_betax + fminus(7)=fminus(7) + avg_alp*sdet*uxyh*psidcm - cons_m(6)*avg_betay + fminus(8)=fminus(8) + avg_alp*sdet*uxzh*psidcm - cons_m(6)*avg_betaz + fplus(6)=fplus(6) + avg_alp*sdet*uxxh*psidcp - cons_p(6)*avg_betax + fplus(7)=fplus(7) + avg_alp*sdet*uxyh*psidcp - cons_p(6)*avg_betay + fplus(8)=fplus(8) + avg_alp*sdet*uxzh*psidcp - cons_p(6)*avg_betaz + psidcfp = avg_alp*cons_p(6)/sdet-avg_betax*psidcp + psidcfm = avg_alp*cons_m(6)/sdet-avg_betax*psidcm endif else if (flux_direction == 2) then call eigenvaluesM(GRHydro_eos_handle,& prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & - mag_m(2),mag_m(3),mag_m(1),& + prim_m(9),prim_m(10),prim_m(8),& lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle, & prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & - mag_p(2),mag_p(3),mag_p(1),& + prim_p(9),prim_p(10),prim_p(8),& lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& - mag_p(2),mag_p(3),mag_p(1),& + cons_p(7),cons_p(8),cons_p(6),& fplus(1),fplus(3),fplus(4),fplus(2),fplus(5),fplus(7),fplus(8),fplus(6),& vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) call num_x_fluxM(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),& - mag_m(2),mag_m(3),mag_m(1),& + cons_m(7),cons_m(8),cons_m(6),& fminus(1),fminus(3),fminus(4),fminus(2),fminus(5),& fminus(7),fminus(8),fminus(6),& vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + avg_alp*sdet*uxyh*psidcm - sdet*mag_m(2)*avg_betax - fminus(7)=fminus(7) + avg_alp*sdet*uyyh*psidcm - sdet*mag_m(2)*avg_betay - fminus(8)=fminus(8) + avg_alp*sdet*uyzh*psidcm - sdet*mag_m(2)*avg_betaz - fplus(6)=fplus(6) + avg_alp*sdet*uxyh*psidcp - sdet*mag_p(2)*avg_betax - fplus(7)=fplus(7) + avg_alp*sdet*uyyh*psidcp - sdet*mag_p(2)*avg_betay - fplus(8)=fplus(8) + avg_alp*sdet*uyzh*psidcp - sdet*mag_p(2)*avg_betaz - psidcfp = avg_alp*mag_p(2)-avg_betay*psidcp - psidcfm = avg_alp*mag_m(2)-avg_betay*psidcm + fminus(6)=fminus(6) + avg_alp*sdet*uxyh*psidcm - cons_m(7)*avg_betax + fminus(7)=fminus(7) + avg_alp*sdet*uyyh*psidcm - cons_m(7)*avg_betay + fminus(8)=fminus(8) + avg_alp*sdet*uyzh*psidcm - cons_m(7)*avg_betaz + fplus(6)=fplus(6) + avg_alp*sdet*uxyh*psidcp - cons_p(7)*avg_betax + fplus(7)=fplus(7) + avg_alp*sdet*uyyh*psidcp - cons_p(7)*avg_betay + fplus(8)=fplus(8) + avg_alp*sdet*uyzh*psidcp - cons_p(7)*avg_betaz + psidcfp = avg_alp*cons_p(7)/sdet-avg_betay*psidcp + psidcfm = avg_alp*cons_m(7)/sdet-avg_betay*psidcm endif else if (flux_direction == 3) then call eigenvaluesM(GRHydro_eos_handle,& prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & - mag_m(3),mag_m(1),mag_m(2),& + prim_m(10),prim_m(8),prim_m(9),& lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call eigenvaluesM(GRHydro_eos_handle,& prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & - mag_p(3),mag_p(1),mag_p(2),& + prim_p(10),prim_p(8),prim_p(9),& lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& - mag_p(3),mag_p(1),mag_p(2),& + cons_p(8),cons_p(6),cons_p(7),& fplus(1),fplus(4),fplus(2),fplus(3),fplus(5),fplus(8),fplus(6),fplus(7), & vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, & avg_det,avg_alp,avg_beta) call num_x_fluxM(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),& - mag_m(3),mag_m(1),mag_m(2),& + cons_m(8),cons_m(6),cons_m(7),& fminus(1),fminus(4),fminus(2),fminus(3),fminus(5), & fminus(8),fminus(6),fminus(7), & vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, & avg_det,avg_alp,avg_beta) if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + avg_alp*sdet*uxzh*psidcm - sdet*mag_m(3)*avg_betax - fminus(7)=fminus(7) + avg_alp*sdet*uyzh*psidcm - sdet*mag_m(3)*avg_betay - fminus(8)=fminus(8) + avg_alp*sdet*uzzh*psidcm - sdet*mag_m(3)*avg_betaz - fplus(6)=fplus(6) + avg_alp*sdet*uxzh*psidcp - sdet*mag_p(3)*avg_betax - fplus(7)=fplus(7) + avg_alp*sdet*uyzh*psidcp - sdet*mag_p(3)*avg_betay - fplus(8)=fplus(8) + avg_alp*sdet*uzzh*psidcp - sdet*mag_p(3)*avg_betaz - psidcfp = avg_alp*mag_p(3)-avg_betaz*psidcp - psidcfm = avg_alp*mag_m(3)-avg_betaz*psidcm + fminus(6)=fminus(6) + avg_alp*sdet*uxzh*psidcm - cons_m(8)*avg_betax + fminus(7)=fminus(7) + avg_alp*sdet*uyzh*psidcm - cons_m(8)*avg_betay + fminus(8)=fminus(8) + avg_alp*sdet*uzzh*psidcm - cons_m(8)*avg_betaz + fplus(6)=fplus(6) + avg_alp*sdet*uxzh*psidcp - cons_p(8)*avg_betax + fplus(7)=fplus(7) + avg_alp*sdet*uyzh*psidcp - cons_p(8)*avg_betay + fplus(8)=fplus(8) + avg_alp*sdet*uzzh*psidcp - cons_p(8)*avg_betaz + psidcfp = avg_alp*cons_p(8)/sdet-avg_betaz*psidcp + psidcfm = avg_alp*cons_m(8)/sdet-avg_betaz*psidcm endif else call CCTK_WARN(0, "Flux direction not x,y,z") end if - + + !!$ Find minimum and maximum wavespeeds charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), & @@ -410,18 +413,23 @@ lamminus(4),lamminus(5)) charpm = charmax - charmin + +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'HLLEM:',fplus(6),fminus(6),fplus(8),fminus(8),charmax,charmin,cons_p(6),cons_m(6) !!$ Calculate flux by standard formula - do m = 1,8 - - qdiff(m) = cons_m(m) - cons_p(m) - - f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & - charmax * charmin * qdiff(m)) / charpm - + do m = 1,8 + + qdiff(m) = cons_m(m) - cons_p(m) + + f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & + charmax * charmin * qdiff(m)) / charpm + end do +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'HLLEM2:',fplus(6),fminus(6),fplus(8),fminus(8),charmax,charmin,qdiff(6),qdiff(8),f1(6),f1(8) + + if(clean_divergence.ne.0) then psidcdiff = psidcm - psidcp @@ -448,6 +456,8 @@ Bconsyflux(i, j, k) = f1(7) Bconszflux(i, j, k) = f1(8) +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'HLLEM2:',fplus(6),fminus(6),fplus(8),fminus(8),f1(6),f1(8) + if(clean_divergence.ne.0) then psidcflux(i,j,k) = psidcf endif File [modified]: GRHydro_Prim2ConM.F90 Delta lines: +2 -2 =================================================================== --- branches/divclean/src/GRHydro_Prim2ConM.F90 2011-04-28 01:28:55 UTC (rev 242) +++ branches/divclean/src/GRHydro_Prim2ConM.F90 2011-04-28 16:47:14 UTC (rev 243) @@ -88,10 +88,10 @@ call prim2conM(GRHydro_eos_handle, gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & avg_detr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& - Bvecxplus(i,j,k),Bvecyplus(i,j,k),Bveczplus(i,j,k), & + Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k), & rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& - Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(i,j,k), & + Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k), & w_lorentzplus(i,j,k)) end do File [modified]: GRHydro_Reconstruct.F90 Delta lines: +6 -0 =================================================================== --- branches/divclean/src/GRHydro_Reconstruct.F90 2011-04-28 01:28:55 UTC (rev 242) +++ branches/divclean/src/GRHydro_Reconstruct.F90 2011-04-28 16:47:14 UTC (rev 243) @@ -219,6 +219,9 @@ call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask) endif + +!!$ write(6,*)'recon0:',Bvecxplus(75,5,5),Bvecxminus(75,5,5) + else if (CCTK_EQUALS(recon_vars,"conservative")) then call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & dens, densplus, densminus, trivial_rp, hydro_excision_mask) @@ -920,6 +923,9 @@ CCTK_EQUALS(recon_method,"ppm")) then if(evolve_mhd.ne.0) then call primitive2conservativeM(CCTK_PASS_FTOF) + +!!$ write(6,*)'recon1:',Bvecxplus(75,5,5),Bvecxminus(75,5,5),Bconsxplus(75,5,5),Bconsxminus(75,5,5) + else call primitive2conservative(CCTK_PASS_FTOF) endif File [modified]: GRHydro_SourceM.F90 Delta lines: +7 -0 =================================================================== --- branches/divclean/src/GRHydro_SourceM.F90 2011-04-28 01:28:55 UTC (rev 242) +++ branches/divclean/src/GRHydro_SourceM.F90 2011-04-28 16:47:14 UTC (rev 243) @@ -442,6 +442,10 @@ srhs(i,j,k,2) = alp(i,j,k)*sqrtdet*sy_source srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sz_source taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source + +!!$ if(i.eq.75.and.j.eq.5.and.k.eq.5)write(6,*)'source:',i,j,k, & +!!$ Bconsrhs(i,j,k,1),Bconsrhs(i,j,k,3) + if(clean_divergence.ne.0) then psidcrhs(i,j,k) = -1.d0 * ( ch_dc*ch_dc/cp_dc/cp_dc*alp(i,j,k) + & dx_betax + dy_betay + dz_betaz ) * psidc(i,j,k) + & @@ -500,6 +504,9 @@ psidc(i,j,k)*alp(i,j,k)*sqrtdet*( uxz*gdg_x + uyz*gdg_y + & uzz*gdg_z ) +!!$ if(i.eq.5.and.j.eq.5.and.k.eq.5)write(6,*)'source:',i,j,k, & +!!$ bvcx_source,bvcz_source + Bconsrhsx(i,j,k) = bvcx_source Bconsrhsy(i,j,k) = bvcy_source Bconsrhsz(i,j,k) = bvcz_source

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/examples/tov/ (Rev. 73)
by knarf＠cct.lsu.edu 28 Apr '11

28 Apr '11

User: knarf Date: 2011/04/28 11:21 AM Modified: /examples/tov/ tov.rpar Log: remove parameter which isn't present anymore in the hg-version of Carpet File Changes: Directory: /examples/tov/ ========================= File [modified]: tov.rpar Delta lines: +0 -1 =================================================================== --- examples/tov/tov.rpar 2011-04-27 22:55:25 UTC (rev 72) +++ examples/tov/tov.rpar 2011-04-28 16:21:09 UTC (rev 73) @@ -84,7 +84,6 @@ Carpet::poison_new_timelevels = "yes" Carpet::check_for_poison = "no" -Carpet::poison_value = 113 Carpet::init_3_timelevels = no Carpet::init_fill_timelevels = "yes"

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits