MPI-AMRVAC  3.0
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_rotating_frame.t
Go to the documentation of this file.
1 !> Module for including rotating frame in hydrodynamics simulations
2 !> The rotation vector is assumed to be along z direction
3 !>(both in cylindrical and spherical)
4 
6  implicit none
7 
8  !> Rotation frequency of the frame
9  double precision :: omega_frame
10 
11  !> Index of the density (in the w array)
12  integer, private, parameter :: rho_ = 1
13 
14 contains
15  !> Read this module's parameters from a file
16  subroutine rotating_frame_params_read(files)
18  character(len=*), intent(in) :: files(:)
19  integer :: n
20 
21  namelist /rotating_frame_list/ omega_frame
22 
23  do n = 1, size(files)
24  open(unitpar, file=trim(files(n)), status="old")
25  read(unitpar, rotating_frame_list, end=111)
26 111 close(unitpar)
27  end do
28 
29  end subroutine rotating_frame_params_read
30 
31  !> Initialize the module
32  subroutine rotating_frame_init()
34  integer :: nwx,idir
35 
37 
38  end subroutine rotating_frame_init
39 
40  !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
41  subroutine rotating_frame_add_source(qdt,ixI^L,ixO^L,wCT,w,x)
43  use mod_usr_methods
44  use mod_geometry
45 
46  integer, intent(in) :: ixI^L, ixO^L
47  double precision, intent(in) :: qdt, x(ixI^S,1:ndim)
48  double precision, intent(in) :: wCT(ixI^S,1:nw)
49  double precision, intent(inout) :: w(ixI^S,1:nw)
50  integer :: idir
51 
52  ! .. local ..
53 
54  double precision :: rotating_terms(ixI^S),frame_omega(ixI^S)
55 
56  select case (coordinate)
57  case (cylindrical)
58  rotating_terms(ixo^s) = omega_frame**2 * x(ixo^s,r_) * wct(ixo^s,iw_rho)
59  if (phi_ > 0) then
60  rotating_terms(ixo^s) = rotating_terms(ixo^s) + 2.d0 * omega_frame *wct(ixo^s,iw_mom(phi_))
61  end if
62  w(ixo^s, iw_mom(r_)) = w(ixo^s, iw_mom(r_)) + qdt * rotating_terms(ixo^s)
63 
64  if(phi_>0 .and. .not. angmomfix) then
65  rotating_terms(ixo^s) = - two * omega_frame * wct(ixo^s,iw_mom(r_))
66  w(ixo^s, iw_mom(phi_)) = w(ixo^s, iw_mom(phi_)) + qdt * rotating_terms(ixo^s)
67  end if
68  case (spherical)
69  !call mpistop("Rotating frame not implemented yet with spherical geometries")
70  ! source[mrad] = 2mphi*vframe/r + rho*vframe**2/r
71 
72  call rotating_frame_omega(x,ixi^l,ixo^l,frame_omega)
73  rotating_terms(ixo^s) = frame_omega(ixo^s)**2 * x(ixo^s,r_) * wct(ixo^s,iw_rho)
74  if (phi_ > 0) then
75  rotating_terms(ixo^s) = rotating_terms(ixo^s) + &
76  2.d0 * frame_omega(ixo^s) * wct(ixo^s,iw_mom(phi_))
77  end if
78  w(ixo^s, iw_mom(r_)) = w(ixo^s, iw_mom(r_)) + qdt * rotating_terms(ixo^s)
79 
80  {^nooned
81  ! source[mtheta] = cot(theta)*(2mphi*vframe/r + rho*vframe**2/r )
82  rotating_terms(ixo^s) = frame_omega(ixo^s)**2 * x(ixo^s,r_) * wct(ixo^s,iw_rho)
83  if (phi_>0 ) then
84  rotating_terms(ixo^s) = rotating_terms(ixo^s) + &
85  2.d0*wct(ixo^s, iw_mom(phi_))* frame_omega(ixo^s)
86  end if
87  w(ixo^s, iw_mom(2)) = w(ixo^s, iw_mom(2)) + qdt * rotating_terms(ixo^s)/ tan(x(ixo^s, 2))
88 
89  if (phi_>0 .and. .not. angmomfix) then
90  ! source[mphi]=-2*mr*vframe/r-2*cot(theta)*mtheta*vframe/r
91  rotating_terms(ixo^s) = -2.d0*frame_omega(ixo^s)* wct(ixo^s, iw_mom(r_))&
92  - 2.d0*wct(ixo^s, iw_mom(2)) * frame_omega(ixo^s)/ tan(x(ixo^s, 2))
93  w(ixo^s, iw_mom(3)) = w(ixo^s, iw_mom(3)) + qdt * rotating_terms(ixo^s)
94  end if
95  }
96  case default
97  call mpistop("Rotating frame not implemented in this geometries")
98  end select
99 
100  end subroutine rotating_frame_add_source
101 
102 
103  subroutine rotating_frame_velocity(x,ixI^L,ixO^L,frame_vel)
105  use mod_usr_methods
106  use mod_geometry
107 
108  integer, intent(in) :: ixI^L, ixO^L
109  double precision, intent(in) :: x(ixI^S,1:ndim)
110  double precision, intent(out) :: frame_vel(ixI^S)
111 
112  frame_vel(ixo^s) = x(ixo^s,r_) * omega_frame
113  {^nooned
114  if (coordinate == spherical .and. ndim > 1) then
115  frame_vel(ixo^s) = frame_vel(ixo^s) * dsin(x(ixo^s,2))
116  end if
117  }
118 
119  end subroutine rotating_frame_velocity
120 
121  subroutine rotating_frame_omega(x,ixI^L,ixO^L,frame_omega)
123  use mod_usr_methods
124  use mod_geometry
125 
126  integer, intent(in) :: ixI^L, ixO^L
127  double precision, intent(in) :: x(ixI^S,1:ndim)
128  double precision, intent(out) :: frame_omega(ixI^S)
129 
130  frame_omega(ixo^s) = omega_frame
131  {^nooned
132  if (coordinate == spherical .and. ndim > 1) then
133  frame_omega(ixo^s) = frame_omega(ixo^s)* dsin(x(ixo^s,2))
134  end if
135  }
136 
137  end subroutine rotating_frame_omega
138 end module mod_rotating_frame
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:194
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer coordinate
Definition: mod_geometry.t:6
integer, parameter spherical
Definition: mod_geometry.t:10
integer, parameter cylindrical
Definition: mod_geometry.t:9
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, parameter unitpar
file handle for IO
logical angmomfix
Enable to strictly conserve the angular momentum (works both in cylindrical and spherical coordinates...
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
Module for including rotating frame in hydrodynamics simulations The rotation vector is assumed to be...
subroutine rotating_frame_params_read(files)
Read this module's parameters from a file.
subroutine rotating_frame_add_source(qdt, ixIL, ixOL, wCT, w, x)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
double precision omega_frame
Rotation frequency of the frame.
subroutine rotating_frame_omega(x, ixIL, ixOL, frame_omega)
subroutine rotating_frame_init()
Initialize the module.
subroutine rotating_frame_velocity(x, ixIL, ixOL, frame_vel)
Module with all the methods that users can customize in AMRVAC.