MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_rotating_frame.t
Go to the documentation of this file.
1 !> Module for including rotating frame in (magneto)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,dtfactor,ixI^L,ixO^L,wCT,w,x)
43  use mod_usr_methods
44  use mod_geometry
46  use mod_comm_lib, only: mpistop
47 
48  integer, intent(in) :: ixI^L, ixO^L
49  double precision, intent(in) :: qdt, dtfactor,x(ixI^S,1:ndim)
50  double precision, intent(in) :: wCT(ixI^S,1:nw)
51  double precision, intent(inout) :: w(ixI^S,1:nw)
52  integer :: idir
53 
54  ! .. local ..
55  double precision :: rotating_terms(ixI^S), frame_omega(ixI^S)
56  double precision :: work(ixI^S)
57 
58  select case (coordinate)
59  case (cylindrical)
60 
61  ! S[mrad] = 2*mphi*Omegaframe + rho*r*Omegaframe**2
62  rotating_terms(ixo^s) = omega_frame**2 * x(ixo^s,r_) * wct(ixo^s,iw_rho)
63 
64  if (phi_ > 0) then
65  rotating_terms(ixo^s) = rotating_terms(ixo^s) + 2.d0 * omega_frame *wct(ixo^s,iw_mom(phi_))
66  end if
67 
68  if(local_timestep) then
69  w(ixo^s, iw_mom(r_)) = w(ixo^s, iw_mom(r_)) + block%dt(ixo^s)*dtfactor * rotating_terms(ixo^s)
70  else
71  w(ixo^s, iw_mom(r_)) = w(ixo^s, iw_mom(r_)) + qdt * rotating_terms(ixo^s)
72  endif
73  ! S[mphi] = -2*mrad*Omegaframe
74  if (phi_ > 0) then
75  rotating_terms(ixo^s) = - 2.0d0*omega_frame * wct(ixo^s,iw_mom(r_))
76  if(local_timestep) then
77  w(ixo^s, iw_mom(phi_)) = w(ixo^s, iw_mom(phi_)) + block%dt(ixo^s)*dtfactor * rotating_terms(ixo^s)
78  else
79  w(ixo^s, iw_mom(phi_)) = w(ixo^s, iw_mom(phi_)) + qdt * rotating_terms(ixo^s)
80  endif
81  end if
82 
83  ! S[etot] = mrad*r*Omegaframe**2
84  if (phys_energy .and. (.not.phys_internal_e)) then
85  if(local_timestep) then
86  w(ixo^s, iw_e) = w(ixo^s, iw_e) + block%dt(ixo^s) *dtfactor * omega_frame**2 * x(ixo^s,r_) * wct(ixo^s,iw_mom(r_))
87  else
88  w(ixo^s, iw_e) = w(ixo^s, iw_e) + qdt * omega_frame**2 * x(ixo^s,r_) * wct(ixo^s,iw_mom(r_))
89  endif
90  endif
91 
92  case (spherical)
93  frame_omega(ixo^s) = omega_frame{^nooned * dsin(x(ixo^s,2))}
94 
95  ! S[mrad] = 2*mphi*Omegaframe + rho*r*Omegaframe**2
96  rotating_terms(ixo^s) = frame_omega(ixo^s)**2 * x(ixo^s,r_) * wct(ixo^s,iw_rho)
97 
98  if (phi_ > 0) then
99  rotating_terms(ixo^s) = rotating_terms(ixo^s) + &
100  2.d0 * frame_omega(ixo^s) * wct(ixo^s,iw_mom(phi_))
101  end if
102  if(local_timestep) then
103  w(ixo^s, iw_mom(r_)) = w(ixo^s, iw_mom(r_)) + block%dt(ixo^s)*dtfactor * rotating_terms(ixo^s)
104  else
105  w(ixo^s, iw_mom(r_)) = w(ixo^s, iw_mom(r_)) + qdt * rotating_terms(ixo^s)
106  endif
107  {^nooned
108  ! S[mtheta] = cot(theta) * S[mrad], reuse above rotating_terms
109  if(local_timestep) then
110  w(ixo^s, iw_mom(2)) = w(ixo^s, iw_mom(2)) + block%dt(ixo^s)*dtfactor * rotating_terms(ixo^s)/ tan(x(ixo^s, 2))
111  else
112  w(ixo^s, iw_mom(2)) = w(ixo^s, iw_mom(2)) + qdt * rotating_terms(ixo^s)/ tan(x(ixo^s, 2))
113  endif
114  ! S[mphi] = -2*Omegaframe * (mrad + cot(theta)*mtheta)
115  if (phi_ > 0) then
116  rotating_terms(ixo^s) = -2.d0*frame_omega(ixo^s)* wct(ixo^s, iw_mom(r_))&
117  - 2.d0*wct(ixo^s, iw_mom(2)) * frame_omega(ixo^s)/ tan(x(ixo^s, 2))
118  if(local_timestep) then
119  w(ixo^s, iw_mom(3)) = w(ixo^s, iw_mom(3)) + block%dt(ixo^s)*dtfactor * rotating_terms(ixo^s)
120  else
121  w(ixo^s, iw_mom(3)) = w(ixo^s, iw_mom(3)) + qdt * rotating_terms(ixo^s)
122  endif
123  end if
124  }
125 
126  ! S[etot] = r*Omegaframe**2 * (mrad + cot(theta)*mtheta)
127  if (phys_energy .and. (.not.phys_internal_e)) then
128  work(ixo^s) = frame_omega(ixo^s)**2 * x(ixo^s,r_) * wct(ixo^s,iw_mom(r_))
129  {^nooned
130  work(ixo^s) = work(ixo^s) + frame_omega(ixo^s)**2 * x(ixo^s,r_) * wct(ixo^s, iw_mom(2))/ tan(x(ixo^s, 2))
131  }
132  if(local_timestep) then
133  w(ixo^s, iw_e) = w(ixo^s, iw_e) + block%dt(ixo^s)*dtfactor* work(ixo^s)
134  else
135  w(ixo^s, iw_e) = w(ixo^s, iw_e) + qdt * work(ixo^s)
136  endif
137  endif
138 
139  case default
140  call mpistop("Rotating frame not implemented in this geometry")
141  end select
142 
143  end subroutine rotating_frame_add_source
144 
145 end module mod_rotating_frame
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer coordinate
Definition: mod_geometry.t:7
integer, parameter spherical
Definition: mod_geometry.t:11
integer, parameter cylindrical
Definition: mod_geometry.t:10
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
integer, parameter unitpar
file handle for IO
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
logical local_timestep
each cell has its own timestep or not
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
logical phys_internal_e
Solve internal energy instead of total energy.
Definition: mod_physics.t:33
logical phys_energy
Solve energy equation or not.
Definition: mod_physics.t:27
Module for including rotating frame in (magneto)hydrodynamics simulations The rotation vector is assu...
subroutine rotating_frame_params_read(files)
Read this module's parameters from a file.
subroutine rotating_frame_add_source(qdt, dtfactor, 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_init()
Initialize the module.
Module with all the methods that users can customize in AMRVAC.