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