Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 296)
by roland.haas＠physics.gatech.edu 02 Nov '11

02 Nov '11

User: rhaas Date: 2011/11/02 10:04 AM Modified: /trunk/ interface.ccl, param.ccl, schedule.ccl /trunk/src/ GRHydro_CalcUpdate.F90, GRHydro_Source.F90, GRHydro_SourceM.F90 Log: First try of a constraint transport scheme for unigrid File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_CalcUpdate.F90 Delta lines: +48 -3 =================================================================== --- trunk/src/GRHydro_CalcUpdate.F90 2011-11-02 13:44:20 UTC (rev 295) +++ trunk/src/GRHydro_CalcUpdate.F90 2011-11-02 15:04:18 UTC (rev 296) @@ -48,9 +48,9 @@ if (use_weighted_fluxes == 0) then !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,alp_l,alp_r,Bvec_l,Bvec_r) - do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil - do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil - do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil + do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! mean we need to compute fluxes for one more point the GRHydro_stencil would indicate + do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil alp_l = 0.5d0 * (alp(i,j,k) + & alp(i-xoffset,j-yoffset,k-zoffset)) @@ -87,6 +87,19 @@ (alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * psidcflux(i,j,k)) * idx endif + if(transport_constraints.ne.0) then + ! Evec lives on edges of cell: Evec(i,j,k,1) is at edge i,j+1/2,k+1/2 ie. the lower-front edge of cell (i,j,k) + if(flux_direction.eq.1) then + Evec(i,j,k,2) = Evec(i,j,k,2) + 0.25d0 * alp_r*(Bconszflux(i,j,k) + Bconszflux(i,j ,k+1)) + Evec(i,j,k,3) = Evec(i,j,k,3) - 0.25d0 * alp_r*(Bconsxflux(i,j,k) + Bconsyflux(i,j+1,k )) + elseif(flux_direction.eq.2) then + Evec(i,j,k,1) = Evec(i,j,k,1) - 0.25d0 * alp_r*(Bconszflux(i,j,k) + Bconszflux(i ,j,k+1)) + Evec(i,j,k,3) = Evec(i,j,k,3) + 0.25d0 * alp_r*(Bconsxflux(i,j,k) + Bconsxflux(i+1,j,k )) + elseif(flux_direction.eq.3) then + Evec(i,j,k,1) = Evec(i,j,k,1) + 0.25d0 * alp_r*(Bconsyflux(i,j,k) + Bconsyflux(i ,j+1,k)) + Evec(i,j,k,2) = Evec(i,j,k,2) - 0.25d0 * alp_r*(Bconsxflux(i,j,k) + Bconsxflux(i+1,j ,k)) + end if + end if if(track_divB.ne.0) then Bvec_l = 0.5d0 * (Bvec(i,j,k,flux_direction) + & Bvec(i-xoffset,j-yoffset,k-zoffset,flux_direction)) @@ -205,6 +218,10 @@ else if (CCTK_EQUALS(method_type, "Flux split FD")) then + if (transport_constraints .ne. 0) then + call CCTK_WARN(0, "Not supported") + end if + do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil @@ -253,6 +270,34 @@ enddo end if + + if (transport_constraints.ne.0 .and. flux_direction.eq.1) then ! HACK: x direction is last + ! RH: I think one could wrap all of this into a single do loop and remove the + ! Evec storage + do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil + do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil + Bconsrhs(i,j,k,1) = & !Bconsrhs(i,j,k,1) & + - 0.5d0 * (Evec(i-1,j-1,k-1,2)-Evec(i-1,j-1,k ,2)) / CCTK_DELTA_SPACE(3) & + - 0.5d0 * (Evec(i-1,j ,k-1,3)-Evec(i-1,j-1,k-1,3)) / CCTK_DELTA_SPACE(2) & + - 0.5d0 * (Evec(i ,j-1,k-1,2)-Evec(i ,j-1,k ,2)) / CCTK_DELTA_SPACE(3) & + - 0.5d0 * (Evec(i ,j ,k-1,3)-Evec(i ,j-1,k-1,3)) / CCTK_DELTA_SPACE(2) + Bconsrhs(i,j,k,2) = & !Bconsrhs(i,j,k,2) & + - 0.5d0 * (Evec(i-1,j-1,k-1,3)-Evec(i-1,j ,k-1,3)) / CCTK_DELTA_SPACE(1) & + - 0.5d0 * (Evec(i ,j-1,k-1,1)-Evec(i-1,j-1,k-1,1)) / CCTK_DELTA_SPACE(3) & + - 0.5d0 * (Evec(i-1,j-1,k ,3)-Evec(i-1,j ,k ,3)) / CCTK_DELTA_SPACE(1) & + - 0.5d0 * (Evec(i ,j-1,k ,1)-Evec(i-1,j-1,k ,1)) / CCTK_DELTA_SPACE(3) + Bconsrhs(i,j,k,3) = & !Bconsrhs(i,j,k,3) & + - 0.5d0 * (Evec(i-1,j-1,k-1,1)-Evec(i ,j-1,k-1,1)) / CCTK_DELTA_SPACE(1) & + - 0.5d0 * (Evec(i-1,j-1,k ,2)-Evec(i-1,j-1,k-1,2)) / CCTK_DELTA_SPACE(2) & + - 0.5d0 * (Evec(i-1,j ,k-1,1)-Evec(i ,j ,k-1,1)) / CCTK_DELTA_SPACE(1) & + - 0.5d0 * (Evec(i-1,j ,k ,2)-Evec(i-1,j ,k-1,2)) / CCTK_DELTA_SPACE(2) + enddo + enddo + enddo + end if + + return File [modified]: GRHydro_Source.F90 Delta lines: +100 -0 =================================================================== --- trunk/src/GRHydro_Source.F90 2011-11-02 13:44:20 UTC (rev 295) +++ trunk/src/GRHydro_Source.F90 2011-11-02 15:04:18 UTC (rev 296) @@ -93,6 +93,11 @@ CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3 CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + ! poison hack + CCTK_INT, SAVE :: reflevel = -1 + CCTK_INT, SAVE :: last_iteration_seen = -1 + CCTK_INT, SAVE :: mol_substep = -1 + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then g11 => gaa g12 => gab @@ -449,6 +454,101 @@ deallocate(force_spatial_second_order) +#if(0) /* poison edges of domain */ + if(last_iteration_seen .ne. cctk_iteration .or. reflevel .ne. grhydro_reflevel) then + last_iteration_seen = cctk_iteration + reflevel = grhydro_reflevel + mol_substep = 0 + else + mol_substep = mol_substep + 1 + end if + do k = 1, GRHydro_stencil*mol_substep + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do k = cctk_lsh(3)-GRHydro_stencil*mol_substep+1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do i = 1, GRHydro_stencil*mol_substep + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do i = cctk_lsh(1)-GRHydro_stencil*mol_substep+1, cctk_lsh(1) + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do j = 1, GRHydro_stencil*mol_substep + do i = 1, cctk_lsh(1) + do k = 1, cctk_lsh(3) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do j = cctk_lsh(2)-GRHydro_stencil*mol_substep+1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + do k = 1, cctk_lsh(3) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do +#endif + + end subroutine SourceTerms File [modified]: GRHydro_SourceM.F90 Delta lines: +7 -0 =================================================================== --- trunk/src/GRHydro_SourceM.F90 2011-11-02 13:44:20 UTC (rev 295) +++ trunk/src/GRHydro_SourceM.F90 2011-11-02 15:04:18 UTC (rev 296) @@ -127,6 +127,13 @@ if (track_divB .ne. 0) then divB=0.d0 endif + + if (transport_constraints .ne. 0) then + Evec = 0.d0 + Bconsrhsx(i,j,k) = 0.d0 + Bconsrhsy(i,j,k) = 0.d0 + Bconsrhsz(i,j,k) = 0.d0 + endif !!$ Set up the array for checking the order. We switch to second order Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +2 -0 =================================================================== --- trunk/interface.ccl 2011-11-02 13:44:20 UTC (rev 295) +++ trunk/interface.ccl 2011-11-02 15:04:18 UTC (rev 296) @@ -347,6 +347,8 @@ real Bcons[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorparity=-1 tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "B-field conservative variable" +real Evec[3] type = GF Timelevels = 1 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "Electric field at edges" + real Y_e_con type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "Conserved electron fraction" real GRHydro_tracers[number_of_tracers] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar"' File [modified]: param.ccl Delta lines: +4 -0 =================================================================== --- trunk/param.ccl 2011-11-02 13:44:20 UTC (rev 295) +++ trunk/param.ccl 2011-11-02 15:04:18 UTC (rev 296) @@ -537,6 +537,10 @@ { } "no" +boolean transport_constraints "Use constraint transport for magnetic field" +{ +} "no" + boolean calculate_bcom "Calculate the comoving contravariant magnetic 4-vector b^a?" { } "no" File [modified]: schedule.ccl Delta lines: +5 -0 =================================================================== --- trunk/schedule.ccl 2011-11-02 13:44:20 UTC (rev 295) +++ trunk/schedule.ccl 2011-11-02 15:04:18 UTC (rev 296) @@ -107,6 +107,11 @@ STORAGE:GRHydro_tracer_rhs } +if (transport_constraints) +{ + STORAGE:Evec +} + # leave this here for compatibility with older code, should not be used anymore # and go out at some point schedule group GRHydro_Initial IN HydroBase_Initial BEFORE SetTmunu

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 295)
by bcmsma＠astro.rit.edu 02 Nov '11

02 Nov '11

User: bmundim Date: 2011/11/02 08:44 AM Modified: /trunk/ interface.ccl, schedule.ccl Log: RIT dev: Fix tensor parity tags for bcom and Bcons. Add prolongation tag for bcom File Changes: Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +4 -4 =================================================================== --- trunk/interface.ccl 2011-11-01 23:25:12 UTC (rev 294) +++ trunk/interface.ccl 2011-11-02 13:44:20 UTC (rev 295) @@ -345,7 +345,7 @@ real scon[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "generalized momenta" -real Bcons[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "B-field conservative variable" +real Bcons[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorparity=-1 tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "B-field conservative variable" real Y_e_con type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "Conserved electron fraction" @@ -367,9 +367,9 @@ real divB type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Magnetic field constraint" -real bcom[3] type = GF Timelevels = 1 tags='Prolongation="none" tensortypealias="U"' "b^i: comoving contravariant magnetic field 4-vector spatial components" -real bcom0 type = GF Timelevels = 1 tags='Prolongation="none" tensortypealias="Scalar"' "b^0 component of the comoving contravariant magnetic field 4-vector" -real bcom_sq type = GF Timelevels = 1 tags='Prolongation="none" tensortypealias="Scalar"' "half of magnectic pressure: contraction of b_a b^a " +real bcom[3] type = GF Timelevels = 1 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="U" tensorparity=-1 interpolator="matter"' "b^i: comoving contravariant magnetic field 4-vector spatial components" +real bcom0 type = GF Timelevels = 1 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" interpolator="matter"' "b^0 component of the comoving contravariant magnetic field 4-vector" +real bcom_sq type = GF Timelevels = 1 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" interpolator="matter"' "half of magnectic pressure: contraction of b_a b^a " ################################################## ### variables in the local tensor basis File [modified]: schedule.ccl Delta lines: +6 -0 =================================================================== --- trunk/schedule.ccl 2011-11-01 23:25:12 UTC (rev 294) +++ trunk/schedule.ccl 2011-11-02 13:44:20 UTC (rev 295) @@ -82,6 +82,12 @@ { STORAGE:divB } + if (calculate_bcom) + { + STORAGE:GRHydro::bcom + STORAGE:GRHydro::bcom0 + STORAGE:GRHydro::bcom_sq + } } STORAGE:GRHydro_eos_scalars STORAGE:GRHydro_minima

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 294)
by knarf＠cct.lsu.edu 01 Nov '11

01 Nov '11

User: knarf Date: 2011/11/01 06:25 PM Removed: /trunk/src/include/ Modified: /trunk/ interface.ccl Log: remove (now) unused old Tmunu interface) File Changes: Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +0 -3 =================================================================== --- trunk/interface.ccl 2011-11-01 23:22:50 UTC (rev 293) +++ trunk/interface.ccl 2011-11-01 23:25:12 UTC (rev 294) @@ -9,9 +9,6 @@ implements: GRHydro inherits: ADMBase, Boundary, SpaceMask, ADMMacros, Tmunubase, HydroBase -#INCLUDES SOURCE: GRHydro_EMTensor.inc in CalcTmunu.inc -#INCLUDES: GRHydro_EMTensor_temps.inc in CalcTmunu_temps.inc - USES INCLUDE: FishEye.h USES INCLUDE: SpaceMask.h USES INCLUDE HEADER: carpet.hh

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 293)
by knarf＠cct.lsu.edu 01 Nov '11

01 Nov '11

User: knarf Date: 2011/11/01 06:22 PM Modified: /trunk/src/ GRHydro_PPM.F90, GRHydro_PPMM.F90 Log: remove unused macros SpaceMask_CheckStateBitsF90_1D File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.F90 Delta lines: +0 -4 =================================================================== --- trunk/src/GRHydro_PPM.F90 2011-11-01 23:21:17 UTC (rev 292) +++ trunk/src/GRHydro_PPM.F90 2011-11-01 23:22:50 UTC (rev 293) @@ -47,10 +47,6 @@ @@*/ -#define SpaceMask_CheckStateBitsF90_1D(mask,i,type_bits,state_bits) \ - (iand(mask((i)),(type_bits)).eq.(state_bits)) - - subroutine SimplePPM_1d(handle,poly,& nx,dx,rho,velx,vely,velz,eps,press,rhominus,& velxminus,velyminus,velzminus,epsminus,rhoplus,velxplus,velyplus,& File [modified]: GRHydro_PPMM.F90 Delta lines: +0 -4 =================================================================== --- trunk/src/GRHydro_PPMM.F90 2011-11-01 23:21:17 UTC (rev 292) +++ trunk/src/GRHydro_PPMM.F90 2011-11-01 23:22:50 UTC (rev 293) @@ -30,10 +30,6 @@ @@*/ -#define SpaceMask_CheckStateBitsF90_1D(mask,i,type_bits,state_bits) \ - (iand(mask((i)),(type_bits)).eq.(state_bits)) - - subroutine SimplePPM_1dM(handle,poly,nx,dx,& rho,velx,vely,velz,Bvcx,Bvcy,Bvcz,psidc,eps,press,& rhominus,velxminus,velyminus,velzminus,Bvcxminus,Bvcyminus,Bvczminus,psidcminus,epsminus,&

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 292)
by knarf＠cct.lsu.edu 01 Nov '11

01 Nov '11

User: knarf Date: 2011/11/01 06:21 PM Modified: /trunk/ schedule.ccl /trunk/src/ GRHydro_CalcUpdate.F90, GRHydro_Con2Prim.F90, GRHydro_Con2PrimM.F90, GRHydro_RegisterMask.c, GRHydro_UpdateMask.F90, GRHydro_UpdateMaskM.F90 Log: dump the atmosphere bitmask in favor of the atmosphere integer mask which is used/allocated already anyway File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_CalcUpdate.F90 Delta lines: +8 -25 =================================================================== --- trunk/src/GRHydro_CalcUpdate.F90 2011-11-01 23:17:55 UTC (rev 291) +++ trunk/src/GRHydro_CalcUpdate.F90 2011-11-01 23:21:17 UTC (rev 292) @@ -41,14 +41,6 @@ CCTK_INT :: i,j,k,itracer CCTK_REAL :: idx, alp_l, alp_r, Bvec_l, Bvec_r - CCTK_INT :: type_bits, atmosphere, not_atmosphere - - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",& - "in_atmosphere") - call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",& - "not_in_atmosphere") - idx = 1.d0 / CCTK_DELTA_SPACE(flux_direction) if (CCTK_EQUALS(method_type, "RSA FV")) then @@ -127,28 +119,19 @@ if (wk_atmosphere .eq. 1) then - if ( (atmosphere_mask(i,j,k) .eq. 1) .or. & - (SpaceMask_CheckStateBitsF90(space_mask,i,j,k,type_bits,atmosphere)) ) then + if (atmosphere_mask(i,j,k) .ne. 0) then !!$ We are in the atmosphere so the momentum flux must vanish srhs(i,j,k,:) = 0.d0 - if ( ( (atmosphere_mask(i-1,j ,k ) .eq. 1) .and. & - (atmosphere_mask(i+1,j ,k ) .eq. 1) .and. & - (atmosphere_mask(i ,j-1,k ) .eq. 1) .and. & - (atmosphere_mask(i ,j+1,k ) .eq. 1) .and. & - (atmosphere_mask(i ,j ,k-1) .eq. 1) .and. & - (atmosphere_mask(i ,j ,k+1) .eq. 1) & - ) .or. & - ( (SpaceMask_CheckStateBitsF90(space_mask,i-1,j ,k ,type_bits,atmosphere)) .and. & - (SpaceMask_CheckStateBitsF90(space_mask,i+1,j ,k ,type_bits,atmosphere)) .and. & - (SpaceMask_CheckStateBitsF90(space_mask,i ,j-1,k ,type_bits,atmosphere)) .and. & - (SpaceMask_CheckStateBitsF90(space_mask,i ,j+1,k ,type_bits,atmosphere)) .and. & - (SpaceMask_CheckStateBitsF90(space_mask,i ,j ,k-1,type_bits,atmosphere)) .and. & - (SpaceMask_CheckStateBitsF90(space_mask,i ,j ,k+1,type_bits,atmosphere)) & - ) & - ) then + if ( (atmosphere_mask(i-1,j ,k ) .ne. 0) .and. & + (atmosphere_mask(i+1,j ,k ) .ne. 0) .and. & + (atmosphere_mask(i ,j-1,k ) .ne. 0) .and. & + (atmosphere_mask(i ,j+1,k ) .ne. 0) .and. & + (atmosphere_mask(i ,j ,k-1) .ne. 0) .and. & + (atmosphere_mask(i ,j ,k+1) .ne. 0) & + ) then !!$ All neighbours are also atmosphere so all rhs vanish File [modified]: GRHydro_Con2Prim.F90 Delta lines: +4 -33 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2011-11-01 23:17:55 UTC (rev 291) +++ trunk/src/GRHydro_Con2Prim.F90 2011-11-01 23:21:17 UTC (rev 292) @@ -61,9 +61,6 @@ logical :: epsnegative character*256 :: warnline - CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits - CCTK_REAL :: local_min_tracer CCTK_REAL :: local_perc_ptol @@ -98,10 +95,6 @@ vup => vel end if - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) @@ -126,7 +119,7 @@ do i = 1, nx !do not compute if in atmosphere or in excised region - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & (hydro_excision_mask(i,j,k) .ne. 0)) cycle epsnegative = .false. @@ -1041,9 +1034,6 @@ logical :: epsnegative character(len=100) warnline - CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits - CCTK_REAL :: local_min_tracer ! save memory when MP is not used @@ -1081,11 +1071,6 @@ call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& GRHydro_rho_min,epsmin,xtemp,xye,pmin,epsmin,keyerr,anyerr) - - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) @@ -1101,7 +1086,7 @@ do i = GRHydro_stencil, nx - GRHydro_stencil + 1 !do not compute if in atmosphere or in an excised region - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle gxxl = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) @@ -1278,9 +1263,6 @@ integer :: i, j, k, itracer, nx, ny, nz CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det - CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits - CCTK_REAL :: local_min_tracer ! character(len=400) :: warnline @@ -1309,10 +1291,6 @@ end if - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) @@ -1334,7 +1312,7 @@ do i = 1, nx !do not compute if in atmosphere or in an excised region - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),\ @@ -1672,9 +1650,6 @@ CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr - CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits - ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 @@ -1695,10 +1670,6 @@ g33 => gzz end if - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) @@ -1708,7 +1679,7 @@ do i = GRHydro_stencil, nx - GRHydro_stencil + 1 !do not compute if in atmosphere or in an excised region - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle gxxl = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +5 -32 =================================================================== --- trunk/src/GRHydro_Con2PrimM.F90 2011-11-01 23:17:55 UTC (rev 291) +++ trunk/src/GRHydro_Con2PrimM.F90 2011-11-01 23:21:17 UTC (rev 292) @@ -55,9 +55,6 @@ CCTK_INT :: epsnegative character(len=100) warnline - CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits - CCTK_REAL :: local_min_tracer, local_gam, local_pgam,local_K,sc ! begin EOS Omni vars @@ -67,10 +64,6 @@ xpress=0.0d0;xtemp=0.0d0;xye=0.0d0;xeps=0.0d0 ! end EOS Omni vars - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) @@ -100,7 +93,7 @@ do i = 1, nx !do not compute if in atmosphere or in excised region - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & (hydro_excision_mask(i,j,k) .ne. 0)) cycle epsnegative = 0 @@ -297,9 +290,6 @@ CCTK_INT :: epsnegative character(len=100) warnline - CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits - CCTK_REAL :: local_min_tracer ! begin EOS Omni vars @@ -320,10 +310,6 @@ 1.d0,1.0d0,xtemp,xye,local_gam,keyerr,anyerr) local_gam=local_gam+1.0 - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) @@ -339,7 +325,7 @@ do i = GRHydro_stencil, nx - GRHydro_stencil + 1 !do not compute if in atmosphere or in an excised region - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) @@ -531,9 +517,6 @@ integer :: i, j, k, itracer, nx, ny, nz CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det,b2 - CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits - CCTK_INT :: epsnegative CCTK_REAL :: local_min_tracer, local_gam, local_pgam,local_K, sc @@ -546,10 +529,6 @@ xpress=0.0d0;xtemp=0.0d0;xye=0.0d0;xeps=0.0d0 ! end EOS Omni vars - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) @@ -572,7 +551,7 @@ do i = 1, nx !do not compute if in atmosphere or in an excised region - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle 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)) @@ -661,10 +640,8 @@ CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr CCTK_REAL :: b2minus, b2plus + CCTK_INT :: epsnegative - CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits, epsnegative - CCTK_REAL :: local_gam, local_pgam,local_K,scplus,scminus ! begin EOS Omni vars @@ -674,10 +651,6 @@ xpress=0.0d0;xtemp=0.0d0;xye=0.0d0;xeps=0.0d0 ! end EOS Omni vars - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) @@ -687,7 +660,7 @@ do i = GRHydro_stencil, nx - GRHydro_stencil + 1 !do not compute if in atmosphere or in an excised region - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) File [modified]: GRHydro_RegisterMask.c Delta lines: +0 -21 =================================================================== --- trunk/src/GRHydro_RegisterMask.c 2011-11-01 23:17:55 UTC (rev 291) +++ trunk/src/GRHydro_RegisterMask.c 2011-11-01 23:21:17 UTC (rev 292) @@ -25,7 +25,6 @@ At the moment, the recognized states and values are - Hydro_Atmosphere (in_atmosphere, not_in_atmosphere) Hydro_RiemannProblem (trivial, not_trivial) @enddesc @@ -44,15 +43,8 @@ int ierr; - const char *atmosphere_list[2] = {"in_atmosphere","not_in_atmosphere"}; const char *rp_list[2] = {"trivial","not_trivial"}; - ierr = SpaceMask_RegisterType("Hydro_Atmosphere", 2, atmosphere_list); - if (ierr) - { - CCTK_WARN(0, "Failed to register the atmosphere with the mask!"); - } - ierr = SpaceMask_RegisterType("Hydro_RiemannProblemX", 2, rp_list); if (ierr) { @@ -74,16 +66,3 @@ return 0; } -void GRHydro_SetupDescriptors(CCTK_ARGUMENTS) -{ - - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - *atmosphere_field_descriptor=SpaceMask_GetTypeBits("Hydro_Atmosphere"); - *atmosphere_atmosp_descriptor=SpaceMask_GetStateBits("Hydro_Atmosphere", - "in_atmosphere"); - *atmosphere_normal_descriptor = SpaceMask_GetStateBits("Hydro_Atmosphere", - "not_in_atmosphere"); -} - File [modified]: GRHydro_UpdateMask.F90 Delta lines: +5 -59 =================================================================== --- trunk/src/GRHydro_UpdateMask.F90 2011-11-01 23:17:55 UTC (rev 291) +++ trunk/src/GRHydro_UpdateMask.F90 2011-11-01 23:21:17 UTC (rev 292) @@ -34,11 +34,6 @@ CCTK_INT :: i,j,k CCTK_REAL :: frac - CCTK_INT :: type_bits, atmosphere - - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - frac = CCTK_DELTA_TIME if(evolve_temper.ne.1.and.evolve_Y_e.ne.1) then @@ -47,16 +42,13 @@ do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) if ( GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0) .or. & - (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \ - type_bits, atmosphere)) .or. & + (atmosphere_mask(i,j,k) .ne. 0) .or. & (tau(i,j,k) + frac * taurhs(i,j,k) .le. 0.d0) .or. & (dens(i,j,k) + frac * densrhs(i,j,k) .le. 0.d0) ) then densrhs(i,j,k) = 0.0d0 srhs(i,j,k,:) = 0.0d0 taurhs(i,j,k) = 0.0d0 atmosphere_mask(i,j,k) = 1 - - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) end if end do end do @@ -69,16 +61,13 @@ do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) if ( GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0) .or. & - (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \ - type_bits, atmosphere)) .or. & + (atmosphere_mask(i,j,k) .ne. 0) .or. & (dens(i,j,k) + frac * densrhs(i,j,k) .le. 0.d0) ) then y_e_con_rhs(i,j,k) = 0.0d0 densrhs(i,j,k) = 0.0d0 srhs(i,j,k,:) = 0.0d0 taurhs(i,j,k) = 0.0d0 atmosphere_mask(i,j,k) = 1 - - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) end if end do end do @@ -110,8 +99,6 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - CCTK_INT :: type_bits, not_atmosphere, i, j, k - ! Initialize all rhs variables and the mask. ! The former vars need to be initialized since there is ! no rhs computation in CCTK_INITIAL or POSTINITIAL. @@ -130,20 +117,6 @@ atmosphere_mask = 0 - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(not_atmosphere, & - &"Hydro_Atmosphere", "not_in_atmosphere") - - do k = 1, cctk_lsh(3) - do j = 1, cctk_lsh(2) - do i = 1, cctk_lsh(1) - - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, not_atmosphere) - - end do - end do - end do - call CCTK_INFO("Setting up the atmosphere mask: all points are not_atmosphere") end subroutine GRHydro_SetupMask @@ -169,23 +142,14 @@ DECLARE_CCTK_ARGUMENTS - CCTK_INT :: type_bits, atmosphere, i, j, k + CCTK_INT :: i,j,k - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, & - &"Hydro_Atmosphere", "in_atmosphere") - do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - if (rho(i,j,k) .le. GRHydro_rho_min) then - - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) atmosphere_mask(i,j,k) = 1 - end if - end do end do end do @@ -225,7 +189,6 @@ CCTK_INT :: i, j, k CCTK_REAL :: det, psi4pt - CCTK_INT :: type_bits, atmosphere, not_atmosphere ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -258,20 +221,12 @@ vup => vel end if - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",& - "in_atmosphere") - call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",& - "not_in_atmosphere") !$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) - if ( (atmosphere_mask(i, j, k) .eq. 1) & - &.or. (SpaceMask_CheckStateBitsF90(space_mask,i, j, k, type_bits,\ - atmosphere)) & - &) then + if (atmosphere_mask(i, j, k) .ne. 0) then rho(i,j,k) = GRHydro_rho_min velx(i,j,k) = 0.0d0 @@ -305,9 +260,7 @@ tau(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)) if (wk_atmosphere .eq. 0) then - atmosphere_mask(i, j, k) = 0 - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\ - not_atmosphere) + atmosphere_mask(i, j, k) = 0 end if endif @@ -344,7 +297,6 @@ CCTK_REAL :: det, psi4pt CCTK_REAL :: rho_min - CCTK_INT :: type_bits, atmosphere, not_atmosphere CCTK_INT :: eos_handle @@ -428,12 +380,6 @@ eos_handle = GRHydro_polytrope_handle - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",& - "in_atmosphere") - call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",& - "not_in_atmosphere") - rho_min = GRHydro_rho_min if (initial_atmosphere_factor .gt. 0) then rho_min = rho_min * initial_atmosphere_factor File [modified]: GRHydro_UpdateMaskM.F90 Delta lines: +1 -25 =================================================================== --- trunk/src/GRHydro_UpdateMaskM.F90 2011-11-01 23:17:55 UTC (rev 291) +++ trunk/src/GRHydro_UpdateMaskM.F90 2011-11-01 23:21:17 UTC (rev 292) @@ -44,27 +44,12 @@ CCTK_INT :: i, j, k CCTK_REAL :: det, psi4pt - CCTK_INT :: type_bits, atmosphere, not_atmosphere - - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",& - "in_atmosphere") - call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",& - "not_in_atmosphere") - do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - if ( (atmosphere_mask(i, j, k) .eq. 1) & - &.or. (SpaceMask_CheckStateBitsF90(space_mask,i, j, k, type_bits,\ - atmosphere)) & - &) then + if ( atmosphere_mask(i, j, k) .eq. 1) then -!!$ write(*,*) 'Resetting at ',i,j,k, atmosphere_mask(i, j, k), & -!!$ & (SpaceMask_CheckStateBitsF90(space_mask,i, j, k, type_bits,\ -!!$ atmosphere)) - rho(i,j,k) = GRHydro_rho_min vel(i,j,k,1) = 0.0d0 vel(i,j,k,2) = 0.0d0 @@ -81,8 +66,6 @@ Bvec(i,j,k,1),Bvec(i,j,k,2),Bvec(i,j,k,3),w_lorentz(i,j,k)) if (wk_atmosphere .eq. 0) then atmosphere_mask(i, j, k) = 0 - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\ - not_atmosphere) end if end if @@ -105,7 +88,6 @@ CCTK_INT :: i, j, k CCTK_REAL :: det, psi4pt - CCTK_INT :: type_bits, atmosphere, not_atmosphere CCTK_INT :: eos_handle ! begin EOS Omni vars @@ -117,12 +99,6 @@ eos_handle = GRHydro_polytrope_handle - call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") - call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",& - "in_atmosphere") - call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",& - "not_in_atmosphere") - do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) Directory: /trunk/ ================== File [modified]: schedule.ccl Delta lines: +1 -5 =================================================================== --- trunk/schedule.ccl 2011-11-01 23:17:55 UTC (rev 291) +++ trunk/schedule.ccl 2011-11-01 23:21:17 UTC (rev 292) @@ -496,11 +496,6 @@ } "Set the atmosphere mask" } -schedule GRHydro_SetupDescriptors AT CCTK_Initial BEFORE HydroBase_Initial -{ - LANG: C -} "Get and store the mask descriptors" - # If using mesh refinement @@ -1140,3 +1135,4 @@ schedule GROUP SetTmunu AT POSTPOSTINITIAL AFTER Con2Prim BEFORE ADMConstraintsGroup { } "Calculate the stress-energy tensor" +

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 291)
by jfaber＠einsteintoolkit.org 01 Nov '11

01 Nov '11

User: jfaber Date: 2011/11/01 06:17 PM Modified: /trunk/ param.ccl Log: Adding calc_bcom parameter File Changes: Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +4 -0 =================================================================== --- trunk/param.ccl 2011-11-01 23:02:41 UTC (rev 290) +++ trunk/param.ccl 2011-11-01 23:17:55 UTC (rev 291) @@ -537,6 +537,10 @@ { } "no" +boolean calculate_bcom "Calculate the comoving contravariant magnetic 4-vector b^a?" +{ +} "no" + ############### Parameters for Tmunu damping (useful in atmosphere) ##################### CCTK_REAL Tmunu_damping_radius_min "damping radius at which we start to damp with a tanh function"

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 290)
by jfaber＠einsteintoolkit.org 01 Nov '11

01 Nov '11

User: jfaber Date: 2011/11/01 06:02 PM Modified: /trunk/ interface.ccl Log: bcom and related quantities added to variable list EOSOMni support for Con2prim added, but commented out File Changes: Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +21 -4 =================================================================== --- trunk/interface.ccl 2011-11-01 22:22:31 UTC (rev 289) +++ trunk/interface.ccl 2011-11-01 23:02:41 UTC (rev 290) @@ -88,6 +88,23 @@ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) +#void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec,CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ +# CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ +# CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ +# CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ +# CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ +# CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ +# CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ +# CCTK_REAL INOUT Bvecsq, \ +# CCTK_REAL INOUT w_lorentz, \ +# CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ +# CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ +# CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ +# CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ +# CCTK_REAL INOUT det, \ +# CCTK_INT OUT epsnegative, \ +# CCTK_REAL OUT retval) + void FUNCTION 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, \ @@ -335,10 +352,6 @@ real Y_e_con type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "Conserved electron fraction" -#real bcom[3] type = GF Timelevels = 3 tags='Prolongation="none" tensortypealias="D" tensorweight=+1.0 interpolator="matter"' "comoving magnetic field components" - -#real bcom0 type = GF Timelevels = 3 tags='Prolongation="none" tensortypealias="Scalar"' "0 component of the comoving magnetic field" - real GRHydro_tracers[number_of_tracers] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar"' { tracer @@ -357,6 +370,10 @@ real divB type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Magnetic field constraint" +real bcom[3] type = GF Timelevels = 1 tags='Prolongation="none" tensortypealias="U"' "b^i: comoving contravariant magnetic field 4-vector spatial components" +real bcom0 type = GF Timelevels = 1 tags='Prolongation="none" tensortypealias="Scalar"' "b^0 component of the comoving contravariant magnetic field 4-vector" +real bcom_sq type = GF Timelevels = 1 tags='Prolongation="none" tensortypealias="Scalar"' "half of magnectic pressure: contraction of b_a b^a " + ################################################## ### variables in the local tensor basis ##################################################

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 131)
by jfaber＠einsteintoolkit.org 01 Nov '11

01 Nov '11

User: jfaber Date: 2011/11/01 06:01 PM Modified: /trunk/ interface.ccl Log: EOSomni support added but commented out File Changes: Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +17 -1 =================================================================== --- trunk/interface.ccl 2011-10-25 18:45:14 UTC (rev 130) +++ trunk/interface.ccl 2011-11-01 23:01:55 UTC (rev 131) @@ -14,7 +14,6 @@ CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ CCTK_REAL OUT det) - SUBROUTINE Prim2ConPoly(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, \ @@ -116,6 +115,23 @@ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) +#void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ +# CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ +# CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ +# CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ +# CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ +# CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ +# CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ +# CCTK_REAL INOUT Bvecsq, \ +# CCTK_REAL INOUT w_lorentz, \ +# CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ +# CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ +# CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ +# CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ +# CCTK_REAL INOUT det, \ +# CCTK_INT OUT epsnegative, \ +# CCTK_REAL OUT retval) + void FUNCTION Con2PrimPolyM(CCTK_INT INOUT handle, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ CCTK_REAL INOUT sc, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 289)
by bcmsma＠astro.rit.edu 01 Nov '11

01 Nov '11

User: bmundim Date: 2011/11/01 05:22 PM Modified: /trunk/src/ GRHydro_Con2PrimM_pt_EOSOmni.c Log: RIT dev: Minor change in sign convention of an equation. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM_pt_EOSOmni.c Delta lines: +7 -7 =================================================================== --- trunk/src/GRHydro_Con2PrimM_pt_EOSOmni.c 2011-11-01 22:21:25 UTC (rev 288) +++ trunk/src/GRHydro_Con2PrimM_pt_EOSOmni.c 2011-11-01 22:22:31 UTC (rev 289) @@ -903,12 +903,12 @@ //Eq. 5: Qdotn + B^2(1+vsq)/2 - QdotBsq/2/W^2 + vsq W + rho_0(vsq) + u = 0 //rho0 = D * sqrt(1-vsq) - resid[1] = lgp->Qdotn + lgp->half_Bsq*(1.0+vsq) - 0.5*QB2Winv2 + vsq*W + rho0 + u; + resid[1] = -lgp->Qdotn - lgp->half_Bsq*(1.0+vsq) + 0.5*QB2Winv2 - vsq*W - rho0 - u; drho0_dv = -0.5*lgp->D / LorInv; - j10 = vsq+QB2Winv3; - j11 = lgp->half_Bsq + W + drho0_dv; + j10 = -vsq-QB2Winv3; + j11 = -lgp->half_Bsq - W - drho0_dv; //Eq. 6: u+ p - W(1-vsq) + rho0 = 0 => p-W = -W vsq - rho0 - u @@ -927,18 +927,18 @@ jac[0][2] = 0.0; jac[1][0] = j10; jac[1][1] = j11; - jac[1][2] = 1.0; + jac[1][2] = -1.0; jac[2][0] = vsq-1.0; jac[2][1] = dpress_dv + W + drho0_dv; jac[2][2] = dpress_du + 1.0; - c00 = j11*jac[2][2]-jac[2][1]; + c00 = j11*jac[2][2]+jac[2][1]; c01 = -1*j01*jac[2][2]; c02 = j01*jac[1][2]; - c10 = jac[2][0]-j10*jac[2][2]; + c10 = -jac[2][0]-j10*jac[2][2]; c11 = j00*jac[2][2]; - c12 = -j00; + c12 = j00; c20 = j10*jac[2][1]-j11*jac[2][0]; c21 = j01*jac[2][0]-j00*jac[2][1];

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 288)
by bcmsma＠astro.rit.edu 01 Nov '11

01 Nov '11

User: bmundim Date: 2011/11/01 05:21 PM Added: /trunk/src/ GRHydro_Con2PrimM_pt_EOSOmni.c Log: RIT dev: New more general EOS Omni interface for Con2Prim. File Changes: Directory: /trunk/src/ ====================== File [added]: GRHydro_Con2PrimM_pt_EOSOmni.c Delta lines: <Unavailable>

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits