12 real(8),
parameter:: pi2=8.0*atan(1.0)
17 procedure(single_solenoidal),
pointer:: selected=>null()
18 character(len=64):: solver
23 procedure:: single_solenoidal
25 integer,
save:: seed=22
31 SUBROUTINE init (self, solver, id, mesh)
33 character(len=64):: solver
35 class(
mesh_t),
dimension(:),
pointer:: mesh
38 character(len=32),
save:: type=
'none' 39 logical,
save:: first_time=.true.
40 namelist /force_params/
type, seed
42 call trace_begin(
'force_t%init')
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)
53 case (
'single_solenoidal')
54 self%selected => single_solenoidal
56 nullify(self%selected)
58 self%solver = trim(solver)
60 call self%random%init (seed)
67 SUBROUTINE dealloc (self)
69 END SUBROUTINE dealloc
75 FUNCTION single_solenoidal (self, time, d, p, Ux, Uy, Uz, m)
RESULT (ff)
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
83 integer :: ix, iy, iz, iostat
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
89 call trace_begin(
'force_t%single_solenoidal')
95 read (io%input, force_solenoidal_params, iostat=iostat)
96 if (io%master)
write (*, force_solenoidal_params)
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()
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))
120 END FUNCTION single_solenoidal
Define code units, in terms of (here) CGS units.
Template module for mesh.
Simple forcing module for tests.