← Back to branch summary

~airpollution/fluidity/fluidity_airpollution

« back to all changes in this revision

Viewing changes to femtools/Streamfunction.F90

  • Committer: ziyouzhj
  • Date: 2013-12-09 16:51:29 UTC
  • Revision ID: ziyouzhj@gmail.com-20131209165129-ucoetc3u0atyy05c
airpolution

Show diffs side-by-side

added added

removed removed

Lines of Context:
femtools/Streamfunction.F90
 
1
!    Copyright (C) 2010 Imperial College London and others.
 
2
!    
 
3
!    Please see the AUTHORS file in the main source directory for a full list
 
4
!    of copyright holders.
 
5
!
 
6
!    Prof. C Pain
 
7
!    Applied Modelling and Computation Group
 
8
!    Department of Earth Science and Engineering
 
9
!    Imperial College London
 
10
!
 
11
!    amcgsoftware@imperial.ac.uk
 
12
!    
 
13
!    This library is free software; you can redistribute it and/or
 
14
!    modify it under the terms of the GNU Lesser General Public
 
15
!    License as published by the Free Software Foundation,
 
16
!    version 2.1 of the License.
 
17
!
 
18
!    This library is distributed in the hope that it will be useful,
 
19
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
 
20
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
21
!    Lesser General Public License for more details.
 
22
!
 
23
!    You should have received a copy of the GNU Lesser General Public
 
24
!    License along with this library; if not, write to the Free Software
 
25
!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
 
26
!    USA
 
27
 
 
28
#include "fdebug.h"
 
29
 
 
30
module streamfunction
 
31
  use state_module
 
32
  use fields
 
33
  use sparse_tools
 
34
  use spud
 
35
  use global_parameters, only: OPTION_PATH_LEN
 
36
  use sparsity_patterns
 
37
  use solvers
 
38
  use boundary_conditions
 
39
  use vector_tools
 
40
  use transform_elements
 
41
  use eventcounter
 
42
  use sparsity_patterns_meshes
 
43
  USE parallel_fields
 
44
  USE Parallel_Tools
 
45
  
 
46
  implicit none
 
47
 
 
48
  private
 
49
  public :: calculate_stream_function_multipath_2d
 
50
 
 
51
  type(integer_vector), dimension(:), allocatable, save :: flux_face_list
 
52
  real, dimension(:,:), allocatable, save :: flux_normal
 
53
 
 
54
  integer, save :: last_adapt=-1
 
55
 
 
56
contains
 
57
 
 
58
  subroutine find_stream_paths(X, streamfunc)
 
59
    !!< Find the paths through the mesh over which the velocity will be
 
60
    !!< integrated to calculate the flux.
 
61
    type(vector_field), intent(in) :: X
 
62
    type(scalar_field), intent(inout) :: streamfunc
 
63
 
 
64
    type(ilist) :: tmp_face_list
 
65
    
 
66
    integer :: bc_count
 
67
    character(len=OPTION_PATH_LEN) :: option_path
 
68
    real, dimension(2) :: start, end, dx,  face_c
 
69
 
 
70
    integer, dimension(:), pointer :: neigh
 
71
    integer :: face, p, ele1, ele2, e2, i
 
72
    real :: dx2, c1, c2
 
73
 
 
74
    bc_count=get_boundary_condition_count(streamfunc)
 
75
       
 
76
    if(.not.(allocated(flux_face_list))) then
 
77
       allocate(flux_face_list(bc_count))
 
78
       CALL nullify(flux_face_list) 
 
79
       allocate(flux_normal(2, bc_count))
 
80
    end if
 
81
 
 
82
 
 
83
    bc_loop: do i=1, bc_count
 
84
       call get_boundary_condition(streamfunc, i, option_path=option_path)
 
85
       
 
86
       if (have_option(trim(option_path)//"/primary_boundary")) then
 
87
          ! No path on the primary boundary.
 
88
 
 
89
          if (associated(flux_face_list(i)%ptr)) then
 
90
             deallocate(flux_face_list(i)%ptr)
 
91
          end if
 
92
          allocate(flux_face_list(i)%ptr(0))
 
93
 
 
94
          cycle bc_loop
 
95
       end if
 
96
       
 
97
       ! We must be on a secondary boundary.
 
98
       call get_option(trim(option_path)//"/secondary_boundary/primary_point"&
 
99
            &, start) 
 
100
       call get_option(trim(option_path)//"/secondary_boundary/secondary_point"&
 
101
            &, end) 
 
102
       
 
103
       dx=-abs(start-end)
 
104
       
 
105
       dx2=dot_product(dx,dx)
 
106
 
 
107
       do ele1=1, element_count(streamfunc)
 
108
          neigh=>ele_neigh(X, ele1)
 
109
          
 
110
          do e2=1,size(neigh)
 
111
             ele2=neigh(e2)
 
112
             ! Don't do boundaries
 
113
             if (ele2<=0) cycle
 
114
             ! Do each edge only once
 
115
             if (ele1>ele2) cycle
 
116
 
 
117
             !for parallel check that we own the node
 
118
             if (.not.element_owned(streamfunc, ele1)) cycle
 
119
 
 
120
             face=ele_face(streamfunc, ele1, ele2)
 
121
 
 
122
             face_c=sum(face_val(X,face),2)/face_loc(X,face)
 
123
 
 
124
             p=dot_product(face_c,dx)/dx2
 
125
 
 
126
             ! If the face is not within the limits of the line, don't do it.
 
127
             if (p<0.or.p>1) cycle
 
128
 
 
129
             c1=cross_product2(sum(ele_val(X,ele1),2)/ele_loc(X,ele1)-start, dx)
 
130
             c2=cross_product2(sum(ele_val(X,ele2),2)/ele_loc(X,ele2)-start, dx)
 
131
                
 
132
             if(C1<0 .and. C2>=0) then
 
133
                continue
 
134
             else if (C1>=0 .and. C2<0) then
 
135
                continue
 
136
             else
 
137
                cycle
 
138
             end if
 
139
 
 
140
             call insert(tmp_face_list, face)
 
141
 
 
142
          end do
 
143
 
 
144
       end do
 
145
 
 
146
       if (associated(flux_face_list(i)%ptr)) then
 
147
          deallocate(flux_face_list(i)%ptr)
 
148
       end if
 
149
       allocate(flux_face_list(i)%ptr(tmp_face_list%length))
 
150
       
 
151
       flux_face_list(i)%ptr=list2vector(tmp_face_list)
 
152
       call flush_list(tmp_face_list)
 
153
 
 
154
       ! Work out the orthonormal to the line.
 
155
       dx=start-end
 
156
       dx=dx/sqrt(dx2)
 
157
       flux_normal(:,i)=(/-dx(2), dx(1)/)
 
158
       
 
159
    end do bc_loop
 
160
    
 
161
  end subroutine find_stream_paths
 
162
 
 
163
  function boundary_value(X, U, bc_num)
 
164
    !!< Calculate the value of the streamfunction on the boundary provided
 
165
    !!< by integrating the velocity flux across a line between this
 
166
    !!< boundary and the primary boundary. 
 
167
    real :: boundary_value
 
168
    type(vector_field), intent(in) :: X, U
 
169
    integer, intent(in) :: bc_num
 
170
    
 
171
    integer :: face, i
 
172
 
 
173
    boundary_value=0.0
 
174
 
 
175
    do i=1, size(flux_face_list(bc_num)%ptr)
 
176
       face=flux_face_list(bc_num)%ptr(i)
 
177
       
 
178
       boundary_value=boundary_value + face_flux(face, X, U, bc_num)
 
179
              
 
180
    end do
 
181
 
 
182
    ! for parallel so each partitition calculates the bit of the flux that it owns and they sum along the boundary so they all have the correct bd. condition
 
183
    call allsum(boundary_value)
 
184
    
 
185
  contains
 
186
 
 
187
    function face_flux(face, X, U, bc_num)
 
188
      real :: face_flux
 
189
      integer, intent(in) :: face, bc_num
 
190
      type(vector_field), intent(in) :: X, U
 
191
      
 
192
      real, dimension(face_ngi(U, face)) :: detwei
 
193
      real, dimension(U%dim,face_ngi(U, face)) :: normal, U_quad
 
194
      integer :: gi
 
195
 
 
196
      call transform_facet_to_physical(X, face, detwei_f=detwei, normal=normal)
 
197
 
 
198
      U_quad=face_val_at_quad(U,face)
 
199
      
 
200
      face_flux=0.0
 
201
 
 
202
      do gi=1, size(detwei)
 
203
         face_flux=face_flux+abs(dot_product(flux_normal(:,bc_num),normal(:,gi)))&
 
204
              * dot_product(U_quad(:,gi),flux_normal(:,bc_num))&
 
205
              * detwei(gi)
 
206
      end do
 
207
 
 
208
    end function face_flux
 
209
 
 
210
  end function boundary_value
 
211
 
 
212
  subroutine calculate_stream_function_multipath_2d(state, streamfunc)
 
213
    !!< Calculate the stream function for a 
 
214
    type(state_type), intent(inout) :: state
 
215
    type(scalar_field), intent(inout) :: streamfunc
 
216
    
 
217
    integer :: i, ele, stat
 
218
    type(vector_field), pointer :: X, U
 
219
    type(csr_sparsity), pointer :: psi_sparsity
 
220
    type(csr_matrix) :: psi_mat
 
221
    type(scalar_field) :: rhs
 
222
    real :: flux_val
 
223
    type(scalar_field), pointer :: surface_field
 
224
 
 
225
    integer :: mesh_movement
 
226
    integer, save :: last_mesh_movement = -1
 
227
 
 
228
    X => extract_vector_field(state, "Coordinate", stat)
 
229
    if(present_and_nonzero(stat)) return
 
230
    U => extract_vector_field(state, "Velocity", stat)
 
231
    if(present_and_nonzero(stat)) return    
 
232
 
 
233
    if (X%dim/=2) then
 
234
       FLExit("Streamfunction is only valid in 2d")
 
235
    end if
 
236
    ! No discontinuous stream functions.
 
237
    if (continuity(streamfunc)<0) then
 
238
       FLExit("Streamfunction must be a continuous field")    
 
239
    end if
 
240
 
 
241
    if (last_adapt<eventcount(EVENT_ADAPTIVITY)) then
 
242
       last_adapt=eventcount(EVENT_ADAPTIVITY)
 
243
       call find_stream_paths(X, streamfunc)
 
244
    end if
 
245
    
 
246
    
 
247
    
 
248
    psi_mat = extract_csr_matrix(state, "StreamFunctionMatrix", stat = stat)
 
249
    if(stat == 0) then
 
250
      mesh_movement = eventcount(EVENT_MESH_MOVEMENT)
 
251
      if(mesh_movement /= last_mesh_movement) then
 
252
        stat = 1
 
253
      end if
 
254
    end if
 
255
    if(stat /= 0) then
 
256
      psi_sparsity => get_csr_sparsity_firstorder(state, streamfunc%mesh, streamfunc%mesh)
 
257
    
 
258
      call allocate(psi_mat, psi_sparsity, name="StreamFunctionMatrix")
 
259
      call zero(psi_mat)
 
260
 
 
261
      call allocate(rhs, streamfunc%mesh, "StreamFunctionRHS")
 
262
      call zero(rhs)
 
263
 
 
264
      do ele=1, element_count(streamfunc)
 
265
       
 
266
        call calculate_streamfunc_ele(rhs, ele, X, U, psi_mat = psi_mat)
 
267
 
 
268
      end do
 
269
 
 
270
      call insert(state, psi_mat, psi_mat%name)
 
271
   else
 
272
      call incref(psi_mat)
 
273
   
 
274
      call allocate(rhs, streamfunc%mesh, "StreamFunctionRHS")
 
275
      call zero(rhs)
 
276
 
 
277
      do ele=1, element_count(streamfunc)
 
278
       
 
279
        call calculate_streamfunc_ele(rhs, ele, X, U)
 
280
 
 
281
      end do
 
282
   end if
 
283
   last_mesh_movement = eventcount(EVENT_MESH_MOVEMENT)
 
284
 
 
285
    do i = 1, get_boundary_condition_count(streamfunc)
 
286
 
 
287
       surface_field=>extract_surface_field(streamfunc, i, "value")
 
288
 
 
289
       flux_val=boundary_value(X,U,i)
 
290
     
 
291
       call set(surface_field, flux_val)
 
292
 
 
293
    end do
 
294
 
 
295
    call zero(streamfunc)
 
296
 
 
297
    call apply_dirichlet_conditions(psi_mat, rhs, streamfunc)
 
298
 
 
299
    call petsc_solve(streamfunc, psi_mat, rhs)
 
300
 
 
301
    call deallocate(rhs)
 
302
    call deallocate(psi_mat)
 
303
 
 
304
  contains
 
305
    
 
306
    subroutine calculate_streamfunc_ele(rhs, ele, X, U, psi_mat)
 
307
      type(scalar_field), intent(inout) :: rhs
 
308
      type(vector_field), intent(in) :: X,U
 
309
      integer, intent(in) :: ele
 
310
      type(csr_matrix), optional, intent(inout) :: psi_mat
 
311
      
 
312
      ! Transformed gradient function for velocity.
 
313
      real, dimension(ele_loc(U, ele), ele_ngi(U, ele), mesh_dim(U)) :: du_t
 
314
      ! Ditto for the stream function, psi
 
315
      real, dimension(ele_loc(rhs, ele), ele_ngi(rhs, ele), mesh_dim(rhs))&
 
316
           & :: dpsi_t 
 
317
 
 
318
      ! Local vorticity_matrix
 
319
      real, dimension(2, ele_loc(rhs, ele), ele_loc(U, ele)) ::&
 
320
           & lvorticity_mat
 
321
      ! Local vorticity
 
322
      real, dimension(ele_loc(rhs, ele)) :: lvorticity
 
323
 
 
324
      ! Variable transform times quadrature weights.
 
325
      real, dimension(ele_ngi(U,ele)) :: detwei
 
326
      
 
327
      type(element_type), pointer :: U_shape, psi_shape
 
328
      integer, dimension(:), pointer :: psi_ele
 
329
      integer :: i
 
330
 
 
331
      U_shape=> ele_shape(U, ele)
 
332
      psi_shape=> ele_shape(rhs, ele)
 
333
      psi_ele=>ele_nodes(rhs, ele)
 
334
      
 
335
      ! Transform U derivatives and weights into physical space.
 
336
      call transform_to_physical(X, ele, U_shape, dshape=du_t, detwei=detwei)
 
337
      ! Ditto psi.
 
338
      call transform_to_physical(X, ele, psi_shape, dshape=dpsi_t)
 
339
 
 
340
      if(present(psi_mat)) then
 
341
        call addto(psi_mat, psi_ele, psi_ele, &
 
342
             dshape_dot_dshape(dpsi_t, dpsi_t, detwei))
 
343
      end if
 
344
 
 
345
      lvorticity_mat=shape_curl_shape_2d(psi_shape, du_t, detwei)
 
346
      
 
347
      lvorticity=0.0
 
348
      do i=1,2
 
349
         lvorticity=lvorticity &
 
350
              +matmul(lvorticity_mat(i,:,:), ele_val(U, i, ele))
 
351
      end do
 
352
      
 
353
      call addto(rhs, psi_ele, -lvorticity)
 
354
      
 
355
    end subroutine calculate_streamfunc_ele
 
356
 
 
357
  end subroutine calculate_stream_function_multipath_2d
 
358
  
 
359
end module streamfunction