← Back to branch summary

~fluidity-core/fluidity/sea-ice-branch

« back to all changes in this revision

Viewing changes to femtools/Diagnostic_Fields.F90

  • Committer: Simon Mouradian
  • Date: 2012-10-19 10:35:59 UTC
  • mfrom: (3520.32.371 fluidity)
  • Revision ID: simon.mouradian06@imperial.ac.uk-20121019103559-y36qa47phc69q8sc
merge from trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
femtools/Diagnostic_Fields.F90
41
41
  use fefields, only: compute_lumped_mass, compute_cv_mass
42
42
  use MeshDiagnostics
43
43
  use spud
44
 
  use CV_Shape_Functions
 
44
  use CV_Shape_Functions, only: make_cv_element_shape, make_cvbdy_element_shape
45
45
  use CV_Faces
46
46
  use CVTools
47
47
  use CV_Upwind_Values
48
 
  use CV_Face_Values
 
48
  use CV_Face_Values, only: evaluate_face_val, theta_val
49
49
  use cv_options
50
50
  use parallel_tools
51
51
  use sparsity_patterns
477
477
    type(scalar_field), intent(inout) :: grn
478
478
 
479
479
    type(vector_field), pointer :: U, X
480
 
    integer :: ele, gi, stat
 
480
    integer :: ele, gi, stat, a, b
481
481
    ! Transformed quadrature weights.
482
482
    real, dimension(ele_ngi(GRN, 1)) :: detwei
483
483
    ! Inverse of the local coordinate change matrix.
525
525
       ! The matmul is as given by dham
526
526
       GRN_q=ele_val_at_quad(U, ele)
527
527
       vis_q=ele_val_at_quad(viscosity, ele)
 
528
 
 
529
       ! for full and partial stress form we need to set the off diagonal terms of the viscosity tensor to zero
 
530
       ! to be able to invert it 
 
531
       if (have_option(trim(U%option_path)//&
 
532
            &"/prognostic/spatial_discretisation/continuous_galerkin"//&
 
533
            &"/stress_terms/stress_form") .or. &
 
534
            have_option(trim(U%option_path)//&
 
535
            &"/prognostic/spatial_discretisation/continuous_galerkin"//&
 
536
            &"/stress_terms/partial_stress_form")) then
 
537
 
 
538
          do a=1,size(vis_q,1)
 
539
             do b=1,size(vis_q,2)
 
540
                if(a.eq.b) cycle
 
541
                vis_q(a,b,:) = 0.0
 
542
             end do
 
543
          end do
 
544
       end if
 
545
 
528
546
       do gi=1, size(detwei)
529
547
          GRN_q(:,gi)=matmul(GRN_q(:,gi), J(:,:,gi))
530
548
          GRN_q(:,gi)=matmul(inverse(vis_q(:,:,gi)), GRN_q(:,gi))
3102
3120
!
3103
3121
      type(state_type), intent(in) :: state
3104
3122
      type(vector_field), intent(inout) :: bed_shear_stress
3105
 
      type(vector_field), pointer :: U
3106
 
      integer, dimension(:), allocatable:: faceglobalnodes
3107
 
      integer :: j,snloc,ele,sele,globnod
 
3123
      type(scalar_field) :: masslump
 
3124
      type(vector_field), pointer :: U, X, bed_shear_stress_surface
 
3125
      type(tensor_field), pointer :: visc
 
3126
      integer, dimension(:), allocatable :: faceglobalnodes
 
3127
      integer :: i,j,snloc,ele,sele,globnod,face
3108
3128
      real :: speed,density,drag_coefficient
 
3129
      real, dimension(mesh_dim(bed_shear_stress)) :: invmass
 
3130
      type(mesh_type) :: surface_mesh
 
3131
      integer, dimension(:), pointer :: surface_node_list
3109
3132
!
 
3133
      
 
3134
      ! assumes constant density
3110
3135
      call get_option(trim(bed_shear_stress%option_path)//"/diagnostic/density", density)
3111
 
      call get_option(trim(bed_shear_stress%option_path)//"/diagnostic/drag_coefficient", drag_coefficient)
3112
 
      
3113
 
      call zero(bed_shear_stress)      
3114
 
      U => extract_vector_field(state, "Velocity")
3115
 
      snloc = face_loc(U, 1)
3116
 
      allocate( faceglobalnodes(1:snloc) )
3117
 
      do sele=1,surface_element_count(U)
3118
 
        ele = face_ele(U, sele)
3119
 
        faceglobalnodes = face_global_nodes(U, sele)
3120
 
        do j = 1,snloc
3121
 
           globnod = faceglobalnodes(j)
3122
 
           speed = norm2(node_val(U, globnod))
3123
 
           call set(bed_shear_stress, globnod, density*drag_coefficient*speed * node_val(U, globnod))
3124
 
        end do
3125
 
      end do
3126
 
      deallocate( faceglobalnodes )
 
3136
 
 
3137
      ! calculate using drag coefficient
 
3138
      if (have_option(trim(bed_shear_stress%option_path)//&
 
3139
           &"/diagnostic/calculation_method/drag_coefficient")) then
 
3140
 
 
3141
         call zero(bed_shear_stress) 
 
3142
 
 
3143
         call get_option(trim(bed_shear_stress%option_path)//&
 
3144
              & "/diagnostic/calculation_method/drag_coefficient",&
 
3145
              & drag_coefficient)
 
3146
 
 
3147
         U => extract_vector_field(state, "Velocity")
 
3148
         snloc = face_loc(U, 1)
 
3149
         allocate( faceglobalnodes(1:snloc) )
 
3150
         do sele=1,surface_element_count(U)
 
3151
            ele = face_ele(U, sele)
 
3152
            faceglobalnodes = face_global_nodes(U, sele)
 
3153
            do j = 1,snloc
 
3154
               globnod = faceglobalnodes(j)
 
3155
               speed = norm2(node_val(U, globnod))
 
3156
               call set(bed_shear_stress, globnod, density*drag_coefficient*speed * node_val(U, globnod))
 
3157
            end do
 
3158
         end do
 
3159
         deallocate( faceglobalnodes )
 
3160
 
 
3161
      ! calculate using velocity gradient
 
3162
      else if (have_option(trim(bed_shear_stress%option_path)//&
 
3163
           &"/diagnostic/calculation_method/wall_treatment")) then
 
3164
 
 
3165
         ! not handled here - handled in Weak_BCs
 
3166
 
 
3167
      ! calculate using velocity gradient
 
3168
      else if (have_option(trim(bed_shear_stress%option_path)//&
 
3169
           &"/diagnostic/calculation_method/velocity_gradient")) then
 
3170
 
 
3171
         call zero(bed_shear_stress) 
 
3172
 
 
3173
         visc => extract_tensor_field(state, "Viscosity")
 
3174
         U    => extract_vector_field(state, "Velocity")
 
3175
         X    => extract_vector_field(state, "Coordinate")
 
3176
         
 
3177
         ! In the CG case we need to calculate a global lumped mass
 
3178
         if(continuity(bed_shear_stress)>=0) then
 
3179
            call allocate(masslump, bed_shear_stress%mesh, 'Masslump')
 
3180
            call zero(masslump)
 
3181
         end if
 
3182
 
 
3183
         do face = 1, surface_element_count(bed_shear_stress)
 
3184
            call calculate_bed_shear_stress_ele(bed_shear_stress, masslump, face, X, U,&
 
3185
                 & visc, density)
 
3186
         end do
 
3187
            
 
3188
         ! In the CG case we globally apply inverse mass
 
3189
         if(continuity(bed_shear_stress)>=0) then
 
3190
            where (masslump%val/=0.0)
 
3191
               masslump%val=1./masslump%val
 
3192
            end where
 
3193
            call scale(bed_shear_stress, masslump)
 
3194
            call deallocate(masslump)
 
3195
         end if 
 
3196
         
 
3197
      else
 
3198
         FLAbort('Unknown bed shear stress calculation method')
 
3199
      end if
 
3200
 
3127
3201
   end subroutine calculate_bed_shear_stress
3128
3202
 
 
3203
   subroutine calculate_bed_shear_stress_ele(bed_shear_stress, masslump, face, X, U, visc&
 
3204
        &, density)
 
3205
 
 
3206
     type(vector_field), intent(inout) :: bed_shear_stress
 
3207
     type(scalar_field), intent(inout) :: masslump
 
3208
     type(vector_field), intent(in), pointer :: X, U
 
3209
     type(tensor_field), intent(in), pointer :: visc
 
3210
     integer, intent(in) :: face
 
3211
     real, intent(in) :: density
 
3212
 
 
3213
     integer :: i, j, i_gi, ele, dim
 
3214
     type(element_type) :: augmented_shape
 
3215
     type(element_type), pointer :: f_shape, shape, X_f_shape, X_shape
 
3216
     real, dimension(X%dim, X%dim, ele_ngi(X, face_ele(X, face))) :: invJ
 
3217
     real, dimension(X%dim, X%dim, face_ngi(X, face)) :: f_invJ  
 
3218
     real, dimension(face_ngi(X, face)) :: detwei
 
3219
     real, dimension(X%dim, face_ngi(X, face)) :: normal, normal_shear_at_quad, X_ele
 
3220
     real, dimension(X%dim) :: abs_normal
 
3221
     real, dimension(ele_loc(X, face_ele(X, face)), face_ngi(X, face), X%dim) :: ele_dshape_at_face_quad
 
3222
     real, dimension(X%dim, X%dim, face_ngi(X, face)) :: grad_U_at_quad, visc_at_quad, shear_at_quad  
 
3223
     real, dimension(X%dim, face_loc(U, face)) :: normal_shear_at_loc
 
3224
     real, dimension(face_loc(X, face), face_loc(U, face)) :: mass
 
3225
 
 
3226
     ele    = face_ele(X, face) ! ele number for volume mesh
 
3227
     dim    = mesh_dim(bed_shear_stress) ! field dimension 
 
3228
 
 
3229
     ! get shape functions
 
3230
     f_shape => face_shape(U, face)     
 
3231
     shape   => ele_shape(U, ele)     
 
3232
     
 
3233
     ! generate shape functions that include quadrature points on the face required
 
3234
     ! check that the shape does not already have these first
 
3235
     assert(shape%degree == 1)
 
3236
     if(associated(shape%dn_s)) then
 
3237
        augmented_shape = shape
 
3238
        call incref(augmented_shape)
 
3239
     else
 
3240
        augmented_shape = make_element_shape(shape%loc, shape%dim, shape%degree, &
 
3241
             & shape%quadrature, quad_s = f_shape%quadrature)
 
3242
     end if
 
3243
    
 
3244
     ! assumes that the jacobian is the same for all quadrature points
 
3245
     ! this is not valid for spheres!
 
3246
     call compute_inverse_jacobian(ele_val(X, ele), ele_shape(X, ele), invj = invJ)
 
3247
     assert(ele_numbering_family(shape) == FAMILY_SIMPLEX)
 
3248
     f_invJ = spread(invJ(:, :, 1), 3, size(f_invJ, 3))
 
3249
      
 
3250
     call transform_facet_to_physical(X, face, detwei_f = detwei, normal = normal)
 
3251
    
 
3252
     ! Evaluate the volume element shape function derivatives at the surface
 
3253
     ! element quadrature points
 
3254
     ele_dshape_at_face_quad = eval_volume_dshape_at_face_quad(augmented_shape, &
 
3255
          & local_face_number(X, face), f_invJ)
 
3256
 
 
3257
     ! Calculate grad of U at the surface element quadrature points
 
3258
     do i=1, dim
 
3259
        do j=1, dim
 
3260
           grad_U_at_quad(i, j, :) = &
 
3261
                & matmul(ele_val(U, j, ele), ele_dshape_at_face_quad(:,:,i))
 
3262
        end do
 
3263
     end do
 
3264
 
 
3265
     visc_at_quad = face_val_at_quad(visc, face)
 
3266
     X_ele = face_val_at_quad(X, face)
 
3267
     do i_gi = 1, face_ngi(X, face)
 
3268
        ! Multiply by visosity
 
3269
        shear_at_quad(:,:,i_gi) = matmul(grad_U_at_quad(:,:,i_gi), visc_at_quad(:,:,i_gi))
 
3270
 
 
3271
        ! Get absolute of normal vector
 
3272
        do i = 1,dim
 
3273
           abs_normal(i) = abs(normal(i,i_gi))
 
3274
        end do
 
3275
 
 
3276
        ! Multiply by surface normal (dim,sgi) to obtain shear in direction normal
 
3277
        ! to surface - transpose (because fluidity stores data in row-major order??)
 
3278
        normal_shear_at_quad(:,i_gi) = matmul(transpose(shear_at_quad(:,:,i_gi)),&
 
3279
             & abs_normal) 
 
3280
     end do  
 
3281
 
 
3282
     normal_shear_at_loc = shape_vector_rhs(f_shape, normal_shear_at_quad, density *&
 
3283
          & detwei)
 
3284
 
 
3285
     mass = shape_shape(f_shape, f_shape, detwei)
 
3286
     ! In the CG case we need to calculate a global lumped mass
 
3287
     if(continuity(bed_shear_stress)>=0) then
 
3288
        call addto(masslump, face_global_nodes(bed_shear_stress,face), sum(mass,1))
 
3289
     else ! In the DG case we will apply the inverse mass locally.
 
3290
        do i = 1,dim
 
3291
           normal_shear_at_loc(i,:) = matmul(inverse(mass), normal_shear_at_loc(i,:))
 
3292
        end do
 
3293
     end if
 
3294
 
 
3295
     ! add to bed_shear_stress field
 
3296
     call addto(bed_shear_stress, face_global_nodes(bed_shear_stress,face), normal_shear_at_loc)
 
3297
 
 
3298
     call deallocate(augmented_shape)
 
3299
          
 
3300
   end subroutine calculate_bed_shear_stress_ele
 
3301
 
3129
3302
   subroutine calculate_max_bed_shear_stress(state, max_bed_shear_stress)
3130
3303
!
3131
3304
      type(state_type), intent(in) :: state