-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathwall_density.f90
159 lines (125 loc) · 4.47 KB
/
wall_density.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
! project : wormlike chain SCFT
! program : SCFT
! source : wall_density.f90
! type : module
! author : Jiuzhou Tang ([email protected])
!purpose :: This module is for building the wall density profile according to
! Glenn's so called MASK skill.
!!!! Warning, this module currently only works for 2D model
module wall
implicit none
!integer,public :: periodic ! periodic sine like wall density for the bottom wall
double precision,public :: periodic ! zhoupan change
double precision,public :: wall_thickness
double precision,public :: zoverdelta
double precision,public :: wave_intensity
integer,public :: wall_at
real*8,public,allocatable :: R_wall_top(:)
real*8,public,allocatable :: R_wall_bottom(:)
real*8,public,allocatable :: R_wall_sum(:)
contains
subroutine Wall_builder()
USE nrtype,only :Pi
USE global_para
USE mpi
USE control
USE constants
USE utility
USE mmpi
USE mpi_fftw3_operation
implicit none
integer :: n,i,j,k,errorw
real*8 :: arg,r,delta_pattern,temp_thickness
!build the upper wall
if (allocated(R_wall_top)) then
else
allocate(R_wall_top(0:local_size-1))
endif
if (allocated(R_wall_bottom)) then
else
allocate(R_wall_bottom(0:local_size-1))
endif
if (allocated(R_wall_sum)) then
else
allocate(R_wall_sum(0:local_size-1))
endif
write(*,*) "wall_thickness=",wall_thickness
write(*,*) "zoverdelta",zoverdelta
wall_at=1
do i=0,LOCAL_SIZE-1
r=(kD(i)%x-1)*dy
arg=-1.0*wall_thickness
arg=arg+r*(1.0-2.0*wall_at)
!arg=arg+wall_at*SIDEx*dx !yiqiande wall
arg=arg+wall_at*SIDEy*dy!zhoupan changes the wall
arg=arg*zoverdelta
R_wall_top(i)=0.5*(1.0-tanh(arg))
enddo
wall_at=0
!build the bottom sine like wall
do i=0,LOCAL_SIZE-1
r=(kD(i)%x-1)*dy
temp_thickness=wall_thickness
temp_thickness=temp_thickness*(wave_intensity*cos(2*Pi*periodic* &
((kD(i)%y-1+local_x_start)/(SIDEx*1.0d0)))+1.0)
!arg=-1.0*wall_thickness
arg=-1.0*temp_thickness
arg=arg+r*(1.0-2.0*wall_at)
arg=arg+wall_at*SIDEy*dy
arg=arg*zoverdelta
R_wall_bottom(i)=0.5*(1.0-tanh(arg))
!R_wall_bottom(i)=R_wall_bottom(i)*(wave_intensity*sin(2*Pi*periodic* &
! ((kD(i)%y-1+local_x_start)/(SIDEx*1.0d0)))+1.0)
!R_wall_top(i)=0
!R_wall_bottom(i)=0
!R_wall_bottom(i)= R_wall_top(i) !zhoupan change to be plate
R_wall_sum(i)=R_wall_top(i)+R_wall_bottom(i)
enddo
call wall_dump()
end subroutine Wall_builder
subroutine wall_dump()
USE nrtype,only :DP
USE global_para
USE mpi
USE control
USE constants
USE utility
USE mmpi
USE mpi_fftw3_operation
implicit none
integer :: i,j,k,aaa,bbb,K_i,k_j,K_k
character(len=30)::aa
character(len=30)::bb
aaa=myid
write(aa,*) aaa
open(unit=23,file= 'Wall_top'//trim(adjustl(aa)) // '.dat',status='replace')
k=0
do K_i=0,local_nx-1
do k_j=0,SIDEy-1
write(23,'(2f10.6,f12.8)') (K_i+local_x_start)*dx,K_j*dy,R_wall_top(k)
k=k+1
enddo
enddo
close(23)
open(unit=23,file= 'Wall_bottom'//trim(adjustl(aa)) // '.dat',status='replace')
k=0
do K_i=0,local_nx-1
do k_j=0,SIDEy-1
write(23,'(2f10.6,f12.8)') (K_i+local_x_start)*dx,K_j*dy,R_wall_bottom(k)
k=k+1
enddo
enddo
close(23)
open(unit=23,file= 'Wall_sum'//trim(adjustl(aa)) // '.dat',status='replace')
k=0
do K_i=0,local_nx-1
do k_j=0,SIDEy-1
write(23,'(2f10.6,f12.8)') (K_i+local_x_start)*dx,K_j*dy,R_wall_sum(k)
k=k+1
enddo
write(23,'(2f10.6,f12.8)')
enddo
close(23)
call mp_barrier()
end subroutine wall_dump
end module wall