← Back to branch summary

~airpollution/fluidity/fluidity_airpollution

« back to all changes in this revision

Viewing changes to assemble/Drag.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:
assemble/Drag.F90
 
1
!    Copyright (C) 2006 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
#include "fdebug.h"
 
28
module drag_module
 
29
use fields_data_types
 
30
use fields
 
31
use state_module
 
32
use boundary_conditions
 
33
use fetools
 
34
use global_parameters, only : OPTION_PATH_LEN
 
35
use spud
 
36
use sparse_tools
 
37
use sparse_tools_petsc
 
38
implicit none
 
39
 
 
40
private
 
41
public drag_surface
 
42
 
 
43
contains
 
44
 
 
45
subroutine drag_surface(bigm, rhs, state, density)
 
46
!!< Applies quadratic or linear drag at boundaries with bc of "drag" type
 
47
!!< This version applies the proper surface integral.
 
48
   type(petsc_csr_matrix), intent(inout):: bigm
 
49
   type(vector_field), intent(inout):: rhs
 
50
   type(state_type), intent(in):: state
 
51
   type(scalar_field), intent(in) :: density
 
52
   
 
53
   type(vector_field), pointer:: velocity, nl_velocity, position, old_velocity
 
54
   type(scalar_field), pointer:: drag_coefficient, distance_top, distance_bottom
 
55
   character(len=OPTION_PATH_LEN) bctype
 
56
   real, dimension(:), allocatable:: face_detwei, coefficient, density_face_gi
 
57
   real, dimension(:,:), allocatable:: drag_mat
 
58
   real dt, theta, gravity_magnitude
 
59
   integer, dimension(:), allocatable:: faceglobalnodes
 
60
   integer, dimension(:), pointer:: surface_element_list
 
61
   integer i, j, k, nobcs, stat
 
62
   integer snloc, sele, sngi
 
63
   logical:: parallel_dg, have_distance_bottom, have_distance_top, have_gravity, manning_strickler
 
64
     
 
65
   ewrite(1,*) 'Inside drag_surface'
 
66
   
 
67
   velocity => extract_vector_field(state, "Velocity")
 
68
   ! velocity at the beginning of the time step
 
69
   old_velocity => extract_vector_field(state, "OldVelocity")
 
70
   ! velocity weighted between old and new with theta
 
71
   nl_velocity => extract_vector_field(state, "NonlinearVelocity")
 
72
   position => extract_vector_field(state, "Coordinate")
 
73
   distance_bottom => extract_scalar_field(state, "DistanceToBottom", stat)
 
74
   have_distance_bottom = stat == 0
 
75
   distance_top => extract_scalar_field(state, "DistanceToTop", stat)
 
76
   have_distance_top = stat == 0
 
77
   
 
78
   call get_option("/timestepping/timestep", dt)
 
79
   call get_option(trim(velocity%option_path)//"/prognostic/temporal_discretisation/theta", &
 
80
                      theta)
 
81
   call get_option("/physical_parameters/gravity/magnitude", gravity_magnitude, &
 
82
          stat=stat)
 
83
   have_gravity = stat == 0
 
84
   parallel_dg=continuity(velocity)<0 .and. IsParallel()
 
85
                      
 
86
   sngi=face_ngi(velocity, 1)
 
87
   snloc=face_loc(velocity,1)
 
88
   
 
89
   allocate(faceglobalnodes(1:snloc), &
 
90
     face_detwei(1:sngi), coefficient(1:sngi), &
 
91
     drag_mat(1:snloc,1:snloc), density_face_gi(1:sngi))
 
92
   
 
93
   nobcs=option_count(trim(velocity%option_path)//'/prognostic/boundary_conditions')
 
94
   do i=1, nobcs
 
95
      call get_boundary_condition(velocity, i, type=bctype, &
 
96
         surface_element_list=surface_element_list)
 
97
      if (bctype=='drag') then
 
98
         manning_strickler=have_option(trim(velocity%option_path)//&
 
99
               '/prognostic/boundary_conditions['//int2str(i-1)//']/type[0]/quadratic_drag/manning-strickler')
 
100
         if (manning_strickler) then
 
101
            if (.not. have_distance_bottom .or. .not. have_distance_top .or. .not. have_gravity) then
 
102
               ewrite(-1,*) "Manning-strickler drag needs DistanceToTop and DistanceToBottom fields and gravity."
 
103
               FLExit("Turn on ocean_boundaries underneath geometry.")
 
104
            end if
 
105
         end if
 
106
         drag_coefficient => extract_scalar_surface_field(velocity, i, "DragCoefficient")
 
107
         do j=1, size(surface_element_list)
 
108
           
 
109
            sele=surface_element_list(j)
 
110
            if (parallel_dg) then
 
111
              if (.not. element_owned(velocity, face_ele(velocity, sele))) cycle
 
112
            end if
 
113
 
 
114
            call transform_facet_to_physical(position, sele, face_detwei)
 
115
            
 
116
            faceglobalnodes=face_global_nodes(nl_velocity, sele)
 
117
            
 
118
            if(have_option(trim(velocity%option_path)//&
 
119
               '/prognostic/boundary_conditions['//int2str(i-1)//']/type[0]/linear_drag')) then
 
120
              ! drag coefficient: C_D
 
121
              coefficient=ele_val_at_quad(drag_coefficient, j)
 
122
            else ! default to quadratic_drag
 
123
              ! drag coefficient: C_D * |u|
 
124
              coefficient=ele_val_at_quad(drag_coefficient, j)* &
 
125
                sqrt(sum(face_val_at_quad(nl_velocity, sele)**2, dim=1))
 
126
              if (manning_strickler) then
 
127
                 ! The manning-strickler formulation takes the form n**2g|u|u/(H**0.3333), where H is the water level, g is gravity and n is the Manning coefficient
 
128
                 ! Note that distance_bottom+distance_top is the current water level H
 
129
                 coefficient=ele_val_at_quad(drag_coefficient, j)*gravity_magnitude*coefficient/((face_val_at_quad(distance_bottom, sele)+face_val_at_quad(distance_top, sele))**(1./3.))
 
130
               end if
 
131
            end if
 
132
               
 
133
            ! density to turn this into a momentum absorption term
 
134
            ! (of course this will just be 1 with boussinesq)
 
135
            density_face_gi = face_val_at_quad(density, sele)
 
136
               
 
137
            drag_mat=shape_shape(face_shape(velocity, sele), &
 
138
               face_shape(velocity, sele), coefficient*face_detwei*density_face_gi)
 
139
               
 
140
            do k=1, velocity%dim
 
141
               call addto(bigm, k, k, faceglobalnodes, faceglobalnodes, &
 
142
                  dt*theta*drag_mat)
 
143
               call addto(rhs, k, faceglobalnodes, &
 
144
                  -matmul(drag_mat, face_val(old_velocity, k, sele)) )
 
145
            end do
 
146
            
 
147
         end do
 
148
            
 
149
      end if
 
150
   end do
 
151
     
 
152
   deallocate(faceglobalnodes, face_detwei, coefficient, drag_mat)
 
153
   
 
154
end subroutine drag_surface
 
155
 
 
156
end module drag_module