DISPATCH
force_mod.f90
1 !===============================================================================
2 !> This is a start on a module that could become a general turbulence forcing
3 !> module, where the standard forcing we have been using in RAMSES is implemented
4 !===============================================================================
5 MODULE force_mod
6  USE io_mod
7  USE trace_mod
8  USE vector_mod
9  USE mesh_mod
10  USE random_mod
11  implicit none
12  private
13  real(8), parameter:: pi2=8.0*atan(1.0)
14  type, public:: force_t
15  integer:: id
16  real(8):: a(3), phase(3)
17  real(8):: t_turn, t_save=0.0
18  type(random_t):: random
19  character(len=64):: solver
20  procedure(single_solenoidal), pointer:: selected=>null()
21  contains
22  procedure:: init
23  procedure:: dealloc
24  end type
25  logical, save:: first_time=.true.
26 CONTAINS
27 
28 !===============================================================================
29 !> Setup a uniform initial state, with density=d0 and B=B0
30 !===============================================================================
31 SUBROUTINE init (self, solver, id, mesh)
32  class(force_t):: self
33  character(len=64):: solver
34  integer:: id, iostat
35  class(mesh_t), pointer:: mesh(:)
36  character(len=32), save:: type='void'
37  integer, save :: seed=22
38  namelist /force_params/ type, seed
39  !----------------------------------------------------------------------------
40  call trace_begin('force_t%init')
41  !$omp critical (input_cr)
42  if (first_time) then
43  first_time = .false.
44  rewind(io%input)
45  read (io%input, force_params, iostat=iostat)
46  if (io%master) write (*, force_params)
47  end if
48  !$omp end critical (input_cr)
49  self%id = id
50  self%solver = trim(solver)
51  call self%random%init (seed)
52  select case (type)
53  case ('single_solenoidal')
54  self%selected => single_solenoidal
55  case ('void')
56  nullify(self%selected)
57  case default
58  print*,'UNKNOWN FORCING TYPE'
59  end select
60  call trace_end
61 END SUBROUTINE init
62 
63 !===============================================================================
64 !> Empty procedure; nothing to deallocate
65 !===============================================================================
66 SUBROUTINE dealloc (self)
67  class(force_t):: self
68 END SUBROUTINE dealloc
69 
70 !===============================================================================
71 !> Forcing on a single 3-D wavenumber, changing amplitude and phase periodically
72 !===============================================================================
73 FUNCTION single_solenoidal (self, time, d, p, Ux, Uy, Uz, m) RESULT (ff)
74  class(force_t) :: self
75  real(8) :: time
76  real, dimension(:,:,:), pointer :: d, ux, uy, uz
77  real, dimension(:,:,:,:), pointer :: p
78  class(mesh_t), dimension(:), pointer :: m
79  real, dimension(m(1)%gn,m(2)%gn,m(3)%gn,3) :: ff
80  !.............................................................................
81  integer :: ix, iy, iz, iostat
82  real(8) :: kx, ky, kz
83  logical, save :: first_time=.true.
84  real, save :: t_turn=0.3, a0(3)=0.1, k(3)=1.0
85  namelist /force_solenoidal_params/ k, a0, t_turn
86  !-----------------------------------------------------------------------------
87  call trace_begin('force_t%single_solenoidal')
88  if (first_time) then
89  !$omp critical (input_cr)
90  if (first_time) then
91  first_time = .false.
92  rewind(io%input)
93  read (io%input, force_solenoidal_params, iostat=iostat)
94  if (io%master) write (*, force_solenoidal_params)
95  end if
96  !$omp end critical (input_cr)
97  end if
98  if (time >= self%t_save) then
99  self%t_save = self%t_save + t_turn
100  self%a = a0*2./3.*(1.0+0.1*self%random%ran3())/t_turn
101  self%phase = pi2*0.5*self%random%ran3()
102  if (self%id==1) print *,'new force at', real(time), self%a, real(self%phase/pi2)
103  end if
104  associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r, l=>m%li, u=>m%ui)
105  do iz=m(3)%lb,m(3)%ub
106  kz = pi2*k(3)*(r3(iz)+m(3)%p)/m(3)%b + self%phase(3)
107  do iy=m(2)%lb,m(2)%ub
108  ky = pi2*k(2)*(r2(iy)+m(2)%p)/m(2)%b + self%phase(2)
109  do ix=m(1)%lb,m(1)%ub
110  kx = pi2*k(1)*(r1(ix)+m(1)%p)/m(1)%b + self%phase(1)
111  ff(ix,iy,iz,1) = self%a(1)*(sin(ky)+sin(kz))
112  ff(ix,iy,iz,2) = self%a(2)*(sin(kz)+sin(kx))
113  ff(ix,iy,iz,3) = self%a(3)*(sin(kx)+sin(ky))
114  end do
115  end do
116  end do
117  end associate
118  call trace_end
119 END FUNCTION single_solenoidal
120 
121 END MODULE
Template module for mesh.
Definition: mesh_mod.f90:4
Definition: io_mod.f90:4
Simple forcing module for tests.
Definition: force_mod.f90:4