13 real(8),
parameter:: pi2=8.0*atan(1.0)
16 real(8):: a(3), phase(3)
17 real(8):: t_turn, t_save=0.0
19 character(len=64):: solver
20 procedure(single_solenoidal),
pointer:: selected=>null()
25 logical,
save:: first_time=.true.
31 SUBROUTINE init (self, solver, id, mesh)
33 character(len=64):: solver
35 class(
mesh_t),
pointer:: mesh(:)
36 character(len=32),
save:: type=
'void' 37 integer,
save :: seed=22
38 namelist /force_params/
type, seed
40 call trace_begin(
'force_t%init')
45 read (io%input, force_params, iostat=iostat)
46 if (io%master)
write (*, force_params)
50 self%solver = trim(solver)
51 call self%random%init (seed)
53 case (
'single_solenoidal')
54 self%selected => single_solenoidal
56 nullify(self%selected)
58 print*,
'UNKNOWN FORCING TYPE' 66 SUBROUTINE dealloc (self)
68 END SUBROUTINE dealloc
73 FUNCTION single_solenoidal (self, time, d, p, Ux, Uy, Uz, m)
RESULT (ff)
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
81 integer :: ix, iy, iz, iostat
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
87 call trace_begin(
'force_t%single_solenoidal')
93 read (io%input, force_solenoidal_params, iostat=iostat)
94 if (io%master)
write (*, force_solenoidal_params)
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)
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))
119 END FUNCTION single_solenoidal
Template module for mesh.
Simple forcing module for tests.