diff --git a/src/problems.F90 b/src/problems.F90 index 93a2aaf..bec4d7f 100644 --- a/src/problems.F90 +++ b/src/problems.F90 @@ -147,9 +147,6 @@ module problems case("current_sheet") setup_problem => setup_problem_current_sheet - case("jet") - setup_problem => setup_problem_jet - case default setup_problem => setup_problem_user @@ -1550,194 +1547,6 @@ module problems !------------------------------------------------------------------------------- ! end subroutine setup_problem_current_sheet -! -!=============================================================================== -! -! subroutine SETUP_PROBLEM_JET: -! ---------------------------- -! -! Subroutine sets the initial conditions for the relativistic jet -! problem. -! -! Arguments: -! -! pdata - pointer to the datablock structure of the currently initialized -! block; -! -!=============================================================================== -! - subroutine setup_problem_jet(pdata) - -! include external procedures and variables -! - use blocks , only : block_data - use coordinates, only : im, jm, km - use coordinates, only : ax, ay, az - use equations , only : prim2cons - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp - use parameters , only : get_parameter_real - -! local variables are not implicit by default -! - implicit none - -! input arguments -! - type(block_data), pointer, intent(inout) :: pdata - -! default parameter values -! - real(kind=8), save :: djet = 1.00d-01 - real(kind=8), save :: damb = 1.00d+01 - real(kind=8), save :: bamb = 1.00d-08 - real(kind=8), save :: pres = 1.00d-02 - real(kind=8), save :: vjet = 0.99d+00 - real(kind=8), save :: bjet = 1.00d-03 - real(kind=8), save :: ljet = 1.00d-00 - real(kind=8), save :: rjet = 1.00d+00 - real(kind=8), save :: rjet2 = 1.00d+00 - -! local saved parameters -! - logical , save :: first = .true. - -! local variables -! - integer :: i, j, k - real(kind=8) :: dx, dy, dz, rm, rr - -! local arrays -! - real(kind=8), dimension(nv,im) :: q, u - real(kind=8), dimension(nv) :: qj - real(kind=8), dimension(im) :: x - real(kind=8), dimension(jm) :: y - real(kind=8), dimension(km) :: z -! -!------------------------------------------------------------------------------- -! -#ifdef PROFILE -! start accounting time for the problem setup -! - call start_timer(imu) -#endif /* PROFILE */ - -! prepare problem constants during the first subroutine call -! - if (first) then - -! get problem parameters -! - call get_parameter_real("djet" , djet) - call get_parameter_real("damb" , damb) - call get_parameter_real("pres" , pres) - call get_parameter_real("bamb" , bamb) - call get_parameter_real("bjet" , bjet) - call get_parameter_real("vjet" , vjet) - call get_parameter_real("ljet" , ljet) - call get_parameter_real("rjet" , rjet) - -! calculate Rjet² -! - rjet2 = rjet * rjet - -! reset the first execution flag -! - first = .false. - - end if ! first call - -! set the conditions inside the jet radius -! - qj(idn) = djet - if (ipr > 0) qj(ipr) = pres - qj(ivx) = vjet - qj(ivy) = 0.0d+00 - qj(ivz) = 0.0d+00 - if (ibx > 0) then - qj(ibx) = 0.0d+00 - qj(iby) = 0.0d+00 - qj(ibz) = bjet - qj(ibp) = 0.0d+00 - end if ! ibx > 0 - -! prepare block coordinates -! - x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im) - dx = x(2) - x(1) - y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm) - dy = y(2) - y(1) -#if NDIMS == 3 - z(1:km) = pdata%meta%zmin + az(pdata%meta%level,1:km) - dz = z(2) - z(1) -#else /* NDIMS == 3 */ - z(1:km) = 0.0d+00 - dz = 0.0d+00 -#endif /* NDIMS == 3 */ - rm = dy * dy + dz * dz - -! iterate over all positions in the YZ plane -! - do k = 1, km - do j = 1, jm - -! calculate radius -! - rr = y(j) * y(j) + z(k) * z(k) - -! set the ambient density, pressure, and velocity -! - q(idn,1:im) = damb - if (ipr > 0) q(ipr,1:im) = pres - q(ivx,1:im) = 0.0d+00 - q(ivy,1:im) = 0.0d+00 - q(ivz,1:im) = 0.0d+00 - -! if magnetic field is present, set it to be uniform with the desired strength -! and orientation -! - if (ibx > 0) then - q(ibx,1:im) = 0.0d+00 - q(iby,1:im) = 0.0d+00 - q(ibz,1:im) = bamb - q(ibp,1:im) = 0.0d+00 - end if ! ibx > 0 - -! set the jet injection -! - if (rr <= max(rm, rjet2)) then - do i = 1, im - if (x(i) <= max(dx, ljet)) then - q(1:nv,i) = qj(1:nv) - end if - end do ! i = 1, im - end if ! R < Rjet - -! convert the primitive variables to conservative ones -! - call prim2cons(im, q(1:nv,1:im), u(1:nv,1:im)) - -! copy the conserved variables to the current block -! - pdata%u(1:nv,1:im,j,k) = u(1:nv,1:im) - -! copy the primitive variables to the current block -! - pdata%q(1:nv,1:im,j,k) = q(1:nv,1:im) - - end do ! j = 1, jm - end do ! k = 1, km - -#ifdef PROFILE -! stop accounting time for the problems setup -! - call stop_timer(imu) -#endif /* PROFILE */ - -!------------------------------------------------------------------------------- -! - end subroutine setup_problem_jet !=============================================================================== ! diff --git a/src/shapes.F90 b/src/shapes.F90 index 30f9fb0..1f462d1 100644 --- a/src/shapes.F90 +++ b/src/shapes.F90 @@ -137,9 +137,6 @@ module shapes case("blast") update_shapes => update_shapes_blast - case("jet") - update_shapes => update_shapes_jet - ! no shape update ! case default @@ -550,176 +547,6 @@ module shapes !------------------------------------------------------------------------------- ! end subroutine update_shapes_blast -! -!=============================================================================== -! -! subroutine UPDATE_SHAPES_JET: -! ---------------------------- -! -! Subroutine resets the primitive and conserved variables within a defined -! shape for the jet problem. -! -! Arguments: -! -! pdata - pointer to the data block structure of the currently initialized -! block; -! time - time at the moment of update; -! dt - time step since the last update; -! -!=============================================================================== -! - subroutine update_shapes_jet(pdata, time, dt) - -! include external procedures and variables -! - use blocks , only : block_data - use coordinates , only : im, jm, km - use coordinates , only : ax, ay, az - use equations , only : prim2cons - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp - use parameters , only : get_parameter_real - -! local variables are not implicit by default -! - implicit none - -! subroutine arguments -! - type(block_data), pointer, intent(inout) :: pdata - real(kind=8) , intent(in) :: time, dt - -! default parameter values -! - real(kind=8), save :: djet = 1.00d-01 - real(kind=8), save :: pres = 1.00d-02 - real(kind=8), save :: vjet = 0.99d+00 - real(kind=8), save :: bjet = 1.00d-03 - real(kind=8), save :: ljet = 1.00d+00 - real(kind=8), save :: rjet = 1.00d+00 - real(kind=8), save :: rjet2 = 1.00d+00 - -! local saved parameters -! - logical , save :: first = .true. - -! local variables -! - integer :: i, j, k - real(kind=8) :: dx, dy, dz, rm, rr - -! local arrays -! - real(kind=8), dimension(nv,im) :: q, u - real(kind=8), dimension(nv) :: qj, uj - real(kind=8), dimension(im) :: x - real(kind=8), dimension(jm) :: y - real(kind=8), dimension(km) :: z -! -!------------------------------------------------------------------------------- -! -#ifdef PROFILE -! start accounting time for the shape update -! - call start_timer(imu) -#endif /* PROFILE */ - -! prepare problem constants during the first subroutine call -! - if (first) then - -! get problem parameters -! - call get_parameter_real("djet" , djet) - call get_parameter_real("pres" , pres) - call get_parameter_real("ljet" , ljet) - call get_parameter_real("bjet" , bjet) - call get_parameter_real("vjet" , vjet) - call get_parameter_real("rjet" , rjet) - -! calculate Rjet² -! - rjet2 = rjet * rjet - -! reset the first execution flag -! - first = .false. - - end if ! first call - -! set the conditions inside the jet radius -! - qj(idn) = djet - if (ipr > 0) qj(ipr) = pres - qj(ivx) = vjet - qj(ivy) = 0.0d+00 - qj(ivz) = 0.0d+00 - if (ibx > 0) then - qj(ibx) = 0.0d+00 - qj(iby) = 0.0d+00 - qj(ibz) = bjet - qj(ibp) = 0.0d+00 - end if ! ibx > 0 - call prim2cons(1, qj(1:nv), uj(1:nv)) - -! prepare block coordinates -! - x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im) - dx = x(2) - x(1) - y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm) - dy = y(2) - y(1) -#if NDIMS == 3 - z(1:km) = pdata%meta%zmin + az(pdata%meta%level,1:km) - dz = z(2) - z(1) -#else /* NDIMS == 3 */ - z(1:km) = 0.0d+00 - dz = 0.0d+00 -#endif /* NDIMS == 3 */ - rm = dy * dy + dz * dz - -! iterate over all positions in the YZ plane -! - do k = 1, km - do j = 1, jm - -! copy the primitive variable vector -! - q(1:nv,1:im) = pdata%q(1:nv,1:im,j,k) - u(1:nv,1:im) = pdata%u(1:nv,1:im,j,k) - -! calculate radius -! - rr = y(j) * y(j) + z(k) * z(k) - - if (rr <= max(rm, rjet2)) then - do i = 1, im - if (x(i) <= max(dx, ljet)) then - q(1:nv,i) = qj(1:nv) - u(1:nv,i) = uj(1:nv) - end if - end do ! i = 1, im - end if ! R < Rjet - -! copy the primitive variables to the current block -! - pdata%q(1:nv,1:im,j,k) = q(1:nv,1:im) - -! copy the conserved variables to the current block -! - pdata%u(1:nv,1:im,j,k) = u(1:nv,1:im) - - end do ! j = 1, jm - end do ! k = 1, km - -#ifdef PROFILE -! stop accounting time for the shape update -! - call stop_timer(imu) -#endif /* PROFILE */ - -!------------------------------------------------------------------------------- -! - end subroutine update_shapes_jet !=============================================================================== !