← Back to branch summary

~fluidity-core/fluidity/bug_1183080

« back to all changes in this revision

Viewing changes to assemble/Momentum_CG.F90

  • Committer: Tim Greaves
  • Date: 2013-04-19 11:53:43 UTC
  • mfrom: (3463.2.94 turbulence-clean)
  • Revision ID: tim.greaves@imperial.ac.uk-20130419115343-qq0m21v905aachuy
This merge is a replay of a previous merge which was reviewed in detail by
Christian Jacobs. There have been no changes since this was approved, hence I'm
fast-tracking this through as the concerns that led to the previous uncommit
have been resolved.

This merge brings in general updates from Jonathan Bull to his turbulence codes
and updates the BFS examples and various related tests to reflect this.

Show diffs side-by-side

added added

removed removed

Lines of Context:
assemble/Momentum_CG.F90
147
147
    ! LES coefficients and options
148
148
    real :: smagorinsky_coefficient
149
149
    character(len=OPTION_PATH_LEN) :: length_scale_type
150
 
    logical :: have_lilly, have_eddy_visc, backscatter
151
 
    logical :: have_strain, have_filtered_strain, have_filter_width
 
150
    logical :: have_eddy_visc, have_filter_width, have_coeff
152
151
 
153
152
    ! Temperature dependent viscosity coefficients:
154
153
    real :: reference_viscosity
249
248
      ! For 4th order:
250
249
      type(tensor_field):: grad_u
251
250
      ! For Germano Dynamic LES:
252
 
      type(vector_field), pointer :: tnu
253
 
      type(tensor_field), pointer :: leonard
254
 
      real                        :: alpha
 
251
      type(vector_field), pointer :: fnu, tnu
 
252
      type(tensor_field), pointer :: leonard, strainprod
 
253
      real                        :: alpha, gamma
255
254
 
256
255
      ! for temperature dependent viscosity :
257
256
      type(scalar_field), pointer :: temperature
382
381
         &"/continuous_galerkin/les_model")
383
382
      if (have_les) then
384
383
         ! Set everything to false initially, then set to true if present
385
 
         have_lilly=.false.; have_eddy_visc=.false.; backscatter=.false.
386
 
         have_strain=.false.; have_filtered_strain=.false.; have_filter_width=.false.
 
384
         have_eddy_visc=.false.; have_filter_width=.false.; have_coeff=.false.
387
385
 
388
386
         les_option_path=(trim(u%option_path)//"/prognostic/spatial_discretisation"//&
389
387
                 &"/continuous_galerkin/les_model")
401
399
            if(have_eddy_visc) then
402
400
               ! Initialise the eddy viscosity field. Calling this subroutine works because
403
401
               ! you can't have 2 different types of LES model for the same material phase.
404
 
               call les_init_diagnostic_fields(state, have_eddy_visc, .false., .false., .false.)
 
402
               call les_init_diagnostic_fields(state, have_eddy_visc, .false., .false.)
405
403
            end if
406
404
         end if
407
405
         if (les_fourth_order) then
415
413
                 smagorinsky_coefficient)
416
414
         end if
417
415
         if(dynamic_les) then
418
 
           ! Are we using the Lilly (1991) modification?
419
 
           have_lilly = have_option(trim(les_option_path)//"/dynamic_les/enable_lilly")
420
 
 
421
 
           ! Whether or not to allow backscatter (negative eddy viscosity)
422
 
           backscatter = have_option(trim(les_option_path)//"/dynamic_les/enable_backscatter")
423
 
 
 
416
           ! Scalar or tensor filter width
 
417
           call get_option(trim(les_option_path)//"/dynamic_les/length_scale_type", length_scale_type)
424
418
           ! Initialise optional diagnostic fields
425
419
           have_eddy_visc = have_option(trim(les_option_path)//"/dynamic_les/tensor_field::EddyViscosity")
426
 
           have_strain = have_option(trim(les_option_path)//"/dynamic_les/tensor_field::StrainRate")
427
 
           have_filtered_strain = have_option(trim(les_option_path)//"/dynamic_les/tensor_field::FilteredStrainRate")
428
420
           have_filter_width = have_option(trim(les_option_path)//"/dynamic_les/tensor_field::FilterWidth")
429
 
           call les_init_diagnostic_fields(state, have_eddy_visc, have_strain, have_filtered_strain, have_filter_width)
 
421
           have_coeff = have_option(trim(les_option_path)//"/dynamic_les/scalar_field::SmagorinskyCoefficient")
 
422
           call les_init_diagnostic_fields(state, have_eddy_visc, have_filter_width, have_coeff)
430
423
 
431
424
           ! Initialise necessary local fields.
432
425
           ewrite(2,*) "Initialising compulsory dynamic LES fields"
433
 
           if(have_option(trim(les_option_path)//"/dynamic_les/vector_field::FilteredVelocity")) then
434
 
             tnu => extract_vector_field(state, "FilteredVelocity")
 
426
           if(have_option(trim(les_option_path)//"/dynamic_les/vector_field::FirstFilteredVelocity")) then
 
427
             fnu => extract_vector_field(state, "FirstFilteredVelocity")
 
428
           else
 
429
             allocate(fnu)
 
430
             call allocate(fnu, u%dim, u%mesh, "FirstFilteredVelocity")
 
431
           end if
 
432
           call zero(fnu)
 
433
           if(have_option(trim(les_option_path)//"/dynamic_les/vector_field::TestFilteredVelocity")) then
 
434
             tnu => extract_vector_field(state, "TestFilteredVelocity")
435
435
           else
436
436
             allocate(tnu)
437
 
             call allocate(tnu, u%dim, u%mesh, "FilteredVelocity")
 
437
             call allocate(tnu, u%dim, u%mesh, "TestFilteredVelocity")
438
438
           end if
439
439
           call zero(tnu)
440
440
           allocate(leonard)
441
441
           call allocate(leonard, u%mesh, "LeonardTensor")
442
442
           call zero(leonard)
 
443
           allocate(strainprod)
 
444
           call allocate(strainprod, u%mesh, "StrainProduct")
 
445
           call zero(strainprod)
443
446
 
444
 
           ! Get (test filter)/(mesh filter) size ratio alpha. Default value is 2.
 
447
           ! Get (first filter)/(mesh size) ratio alpha. Default value is 2.
445
448
           call get_option(trim(les_option_path)//"/dynamic_les/alpha", alpha, default=2.0)
 
449
           ! Get (test filter)/(first filter) size ratio alpha. Default value is 2.
 
450
           call get_option(trim(les_option_path)//"/dynamic_les/gama", gamma, default=2.0)
446
451
 
447
452
           ! Calculate test-filtered velocity field and Leonard tensor field.
448
453
           ewrite(2,*) "Calculating test-filtered velocity and Leonard tensor"
449
 
           call leonard_tensor(nu, x, tnu, leonard, alpha, les_option_path)
 
454
           call leonard_tensor(nu, x, fnu, tnu, leonard, strainprod, alpha, gamma, les_option_path)
450
455
 
451
456
           ewrite_minmax(leonard)
 
457
           ewrite_minmax(strainprod)
452
458
         else
 
459
            fnu => dummyvector
453
460
            tnu => dummyvector
454
461
            leonard => dummytensor
 
462
            strainprod => dummytensor
455
463
         end if
456
464
      else
457
465
         les_second_order=.false.; les_fourth_order=.false.; wale=.false.; dynamic_les=.false.
458
 
         tnu => dummyvector; leonard => dummytensor
 
466
         fnu => dummyvector; tnu => dummyvector; leonard => dummytensor; strainprod => dummytensor
459
467
      end if
460
468
      
461
469
 
678
686
              density, ct_rhs, &
679
687
              source, absorption, buoyancy, gravity, &
680
688
              viscosity, grad_u, &
681
 
              tnu, leonard, alpha, &
 
689
              fnu, tnu, leonard, strainprod, alpha, gamma, &
682
690
              gp, surfacetension, &
683
691
              assemble_ct_matrix_here, depth, &
684
692
              alpha_u_field, abs_wd, temperature, nvfrac, &
815
823
      ewrite_minmax(rhs)
816
824
      
817
825
      if(les_second_order .or. dynamic_les) then
818
 
         call les_solve_diagnostic_fields(state, have_eddy_visc, have_strain, have_filtered_strain, have_filter_width)
 
826
         call les_solve_diagnostic_fields(state, have_eddy_visc, have_filter_width, have_coeff)
819
827
      end if
820
828
 
821
829
      if (les_fourth_order) then
823
831
      end if
824
832
 
825
833
      if (dynamic_les) then
826
 
        if(.not. have_option(trim(les_option_path)//"/dynamic_les/vector_field::FilteredVelocity")) then
 
834
        if(.not. have_option(trim(les_option_path)//"/dynamic_les/vector_field::FirstFilteredVelocity")) then
827
835
          call deallocate(tnu); deallocate(tnu)
828
836
        end if
 
837
        if(.not. have_option(trim(les_option_path)//"/dynamic_les/vector_field::TestFilteredVelocity")) then
 
838
          call deallocate(fnu); deallocate(fnu)
 
839
        end if
829
840
        call deallocate(leonard); deallocate(leonard)
 
841
        call deallocate(strainprod); deallocate(strainprod)
830
842
      end if
831
843
 
832
844
      call deallocate(dummytensor)
1115
1127
                                            density, ct_rhs, &
1116
1128
                                            source, absorption, buoyancy, gravity, &
1117
1129
                                            viscosity, grad_u, &
1118
 
                                            tnu, leonard, alpha, &
 
1130
                                            fnu, tnu, leonard, strainprod, alpha, gamma, &
1119
1131
                                            gp, surfacetension, &
1120
1132
                                            assemble_ct_matrix_here, depth, &
1121
1133
                                            alpha_u_field, abs_wd, temperature, nvfrac, supg_shape)
1143
1155
      type(tensor_field), intent(in) :: viscosity, grad_u
1144
1156
 
1145
1157
      ! Fields for Germano Dynamic LES Model
1146
 
      type(vector_field), intent(in)    :: tnu
1147
 
      type(tensor_field), intent(in)    :: leonard
1148
 
      real, intent(in)                  :: alpha
 
1158
      type(vector_field), intent(in)    :: fnu, tnu
 
1159
      type(tensor_field), intent(in)    :: leonard, strainprod
 
1160
      real, intent(in)                  :: alpha, gamma
1149
1161
 
1150
1162
      type(scalar_field), intent(in) :: gp
1151
1163
      type(tensor_field), intent(in) :: surfacetension
1354
1366
      ! Viscous terms
1355
1367
      if(have_viscosity .or. have_les) then
1356
1368
        call add_viscosity_element_cg(state, ele, test_function, u, oldu_val, nu, x, viscosity, grad_u, &
1357
 
           tnu, leonard, alpha, du_t, detwei, big_m_tensor_addto, rhs_addto, temperature, density, nvfrac)
 
1369
           fnu, tnu, leonard, strainprod, alpha, gamma, du_t, detwei, big_m_tensor_addto, rhs_addto, temperature, density, nvfrac)
1358
1370
      end if
1359
1371
      
1360
1372
      ! Coriolis terms
1970
1982
    end subroutine add_absorption_element_cg
1971
1983
      
1972
1984
    subroutine add_viscosity_element_cg(state, ele, test_function, u, oldu_val, nu, x, viscosity, grad_u, &
1973
 
        tnu, leonard, alpha, &
 
1985
        fnu, tnu, leonard, strainprod, alpha, gamma, &
1974
1986
         du_t, detwei, big_m_tensor_addto, rhs_addto, temperature, density, nvfrac)
1975
1987
      type(state_type), intent(inout) :: state
1976
1988
      integer, intent(in) :: ele
1982
1994
      type(tensor_field), intent(in) :: grad_u
1983
1995
 
1984
1996
      ! Fields for Germano Dynamic LES Model
1985
 
      type(vector_field), intent(in)    :: tnu
1986
 
      type(tensor_field), intent(in)    :: leonard
1987
 
      real, intent(in)                  :: alpha
 
1997
      type(vector_field), intent(in)    :: fnu, tnu
 
1998
      type(tensor_field), intent(in)    :: leonard, strainprod
 
1999
      real, intent(in)                  :: alpha, gamma
1988
2000
 
1989
2001
      ! Local quantities specific to Germano Dynamic LES Model
1990
 
      real                                           :: numerator, denominator
1991
 
      real, dimension(x%dim, x%dim, ele_ngi(u,ele))  :: strain_gi, t_strain_gi
 
2002
      real, dimension(x%dim, x%dim, ele_ngi(u,ele))  :: strain_gi, t_strain_gi, strainprod_gi
1992
2003
      real, dimension(x%dim, x%dim, ele_ngi(u,ele))  :: mesh_size_gi, leonard_gi
 
2004
      real, dimension(x%dim, x%dim, ele_ngi(u,ele))  :: f_tensor, t_tensor
 
2005
      real, dimension(x%dim, x%dim)                  :: mij
1993
2006
      real, dimension(ele_ngi(u, ele))               :: strain_mod, t_strain_mod
 
2007
      real, dimension(ele_ngi(u, ele))               :: f_scalar, t_scalar
1994
2008
      type(element_type)                             :: shape_nu
1995
2009
      integer, dimension(:), pointer                 :: nodes_nu
1996
2010
 
2090
2104
            ! you can't have 2 different types of LES model for the same material_phase.
2091
2105
            if(have_eddy_visc) then
2092
2106
              call les_assemble_diagnostic_fields(state, u, ele, detwei, &
2093
 
                   les_tensor_gi, les_tensor_gi, les_tensor_gi, les_tensor_gi, &
2094
 
                 have_eddy_visc, .false., .false., .false.)
 
2107
                   les_tensor_gi, les_tensor_gi, les_coef_gi, &
 
2108
                 have_eddy_visc, .false., .false.)
2095
2109
            end if
2096
2110
 
2097
2111
         ! Fourth order Smagorinsky model
2112
2126
            nodes_nu => ele_nodes(nu, ele)
2113
2127
            les_tensor_gi=0.0
2114
2128
 
2115
 
            ! Get strain S1 for unfiltered velocity (dim,dim,ngi)
2116
 
            strain_gi = les_strain_rate(du_t, ele_val(nu, ele))
2117
 
            ! Get strain S2 for test-filtered velocity (dim,dim,ngi)
 
2129
            ! Get strain S1 for first-filtered velocity
 
2130
            strain_gi = les_strain_rate(du_t, ele_val(fnu, ele))
 
2131
            ! Get strain S2 for test-filtered velocity
2118
2132
            t_strain_gi = les_strain_rate(du_t, ele_val(tnu, ele))
2119
 
            ! Filter width G1 associated with mesh size (units length^2)
 
2133
            ! Mesh size (units length^2)
2120
2134
            mesh_size_gi = length_scale_tensor(du_t, ele_shape(u, ele))
2121
 
            ! Leonard tensor L at gi
 
2135
            ! Leonard tensor and strain product at Gauss points
2122
2136
            leonard_gi = ele_val_at_quad(leonard, ele)
 
2137
            strainprod_gi = ele_val_at_quad(strainprod, ele)
2123
2138
 
2124
2139
            do gi=1, ele_ngi(nu, ele)
2125
 
              ! Get strain modulus |S1| for unfiltered velocity (ngi)
2126
 
              strain_mod(gi) = sqrt( 2*sum(strain_gi(:,:,gi)*strain_gi(:,:,gi) ) )
2127
 
              ! Get strain modulus |S2| for test-filtered velocity (ngi)
2128
 
              t_strain_mod(gi) = sqrt( 2*sum(t_strain_gi(:,:,gi)*t_strain_gi(:,:,gi) ) )
 
2140
               ! Get strain modulus |S1| for first-filtered velocity
 
2141
               strain_mod(gi) = sqrt( 2*sum(strain_gi(:,:,gi)*strain_gi(:,:,gi) ) )
 
2142
               ! Get strain modulus |S2| for test-filtered velocity
 
2143
               t_strain_mod(gi) = sqrt( 2*sum(t_strain_gi(:,:,gi)*t_strain_gi(:,:,gi) ) )
2129
2144
            end do
2130
2145
 
2131
 
            ! If sum of strain components = 0, don't use dynamic LES model
2132
 
            if(abs(sum(strain_gi(:,:,:))) < epsilon(0.0)) then
2133
 
              les_tensor_gi = 0.0
2134
 
              t_strain_mod = 0.0
2135
 
            else
2136
 
              ! Choose original Germano model or Lilly's (1991) modification from options
2137
 
              if(.not. have_lilly) then
2138
 
                do gi=1, ele_ngi(nu, ele)
2139
 
                  ! |S1|*L.S1
2140
 
                  numerator = sum( leonard_gi(:,:,gi)*strain_gi(:,:,gi) )*strain_mod(gi)
2141
 
 
2142
 
                  ! alpha^2*|S2|*S2.S1
2143
 
                  ! This term is WRONG until I find a way of filtering the strain rate product. The difference may be quite small though.
2144
 
                  denominator = -alpha**2*t_strain_mod(gi)*sum(t_strain_gi(:,:,gi)*strain_gi(:,:,gi))
2145
 
 
2146
 
                  ! Dynamic eddy viscosity m_ij = C*S1
2147
 
                  les_tensor_gi(:,:,gi) = numerator/denominator
2148
 
 
2149
 
                  ! Whether or not to allow negative eddy viscosity (backscattering)
2150
 
                  ! but do not allow (viscosity+eddy_viscosity) < 0.
2151
 
                  if(any(les_tensor_gi(:,:,gi) < 0.0)) then
2152
 
                    if(backscatter) then
2153
 
                      les_tensor_gi(:,:,gi) = max(les_tensor_gi(:,:,gi), epsilon(0.0) - viscosity_gi(:,:,gi))
2154
 
                    else
2155
 
                      les_tensor_gi(:,:,gi) = max(les_tensor_gi(:,:,gi),0.0)
2156
 
                    end if
2157
 
                  end if
2158
 
                end do
2159
 
              else if(have_lilly) then
2160
 
                do gi=1, ele_ngi(nu, ele)
2161
 
                  ! |S1|*L.S1
2162
 
                  numerator = sum(leonard_gi(:,:,gi)*t_strain_gi(:,:,gi))*strain_mod(gi)
2163
 
                  ! alpha^2*|S2|^2*S2.S2
2164
 
                  ! This term is WRONG until I find a way of filtering the strain rate product. The difference may be quite small though.
2165
 
                  denominator = -alpha**2*(t_strain_mod(gi))**2*sum(t_strain_gi(:,:,gi)*t_strain_gi(:,:,gi))
2166
 
                  ! Dynamic eddy viscosity m_ij
2167
 
                  les_tensor_gi(:,:,gi) = numerator/denominator
2168
 
 
2169
 
                  ! Whether or not to allow negative eddy viscosity (backscattering)
2170
 
                  ! but do not allow (viscosity+eddy_viscosity) < 0.
2171
 
                  if(any(les_tensor_gi(:,:,gi) < 0.0)) then
2172
 
                    if(backscatter) then
2173
 
                      les_tensor_gi(:,:,gi) = max(les_tensor_gi(:,:,gi), epsilon(0.0) - viscosity_gi(:,:,gi))
2174
 
                    else
2175
 
                      les_tensor_gi(:,:,gi) = max(les_tensor_gi(:,:,gi),0.0)
2176
 
                    end if
2177
 
                  end if
2178
 
                end do
2179
 
              end if
2180
 
            end if
 
2146
            select case(length_scale_type)
 
2147
               case("scalar")
 
2148
                  ! Scalar first filter width G1 = alpha^2*meshsize (units length^2)
 
2149
                  f_scalar = alpha**2*length_scale_scalar(x, ele)
 
2150
                  ! Combined width G2 = (1+gamma^2)*G1
 
2151
                  t_scalar = (1.0+gamma**2)*f_scalar
 
2152
                  do gi=1, ele_ngi(nu, ele)
 
2153
                    ! Tensor M_ij = (|S2|*S2)G2 - ((|S1|S1)^f2)G1
 
2154
                    mij = t_strain_mod(gi)*t_strain_gi(:,:,gi)*t_scalar(gi) - strainprod_gi(:,:,gi)*f_scalar(gi)
 
2155
                    ! Model coeff C_S = -(L_ij M_ij) / 2(M_ij M_ij)
 
2156
                    les_coef_gi(gi) = -0.5*sum(leonard_gi(:,:,gi)*mij) / sum(mij*mij)
 
2157
                    ! Constrain C_S to be between 0 and 0.04.
 
2158
                    les_coef_gi(gi) = min(max(les_coef_gi(gi),0.0), 0.04)
 
2159
                    ! Isotropic tensor dynamic eddy viscosity = -2C_S|S1|.alpha^2.G1
 
2160
                    les_tensor_gi(:,:,gi) = 2*alpha**2*les_coef_gi(gi)*strain_mod(gi)*f_scalar(gi)
 
2161
                  end do
 
2162
               case("tensor")
 
2163
                  ! First filter width G1 = alpha^2*mesh size (units length^2)
 
2164
                  f_tensor = alpha**2*mesh_size_gi
 
2165
                  ! Combined width G2 = (1+gamma^2)*G1
 
2166
                  t_tensor = (1.0+gamma**2)*f_tensor
 
2167
                  do gi=1, ele_ngi(nu, ele)
 
2168
                    ! Tensor M_ij = (|S2|*S2).G2 - ((|S1|S1)^f2).G1
 
2169
                    mij = t_strain_mod(gi)*t_strain_gi(:,:,gi)*t_tensor(:,:,gi) - strainprod_gi(:,:,gi)*f_tensor(:,:,gi)
 
2170
                    ! Model coeff C_S = -(L_ij M_ij) / 2(M_ij M_ij)
 
2171
                    les_coef_gi(gi) = -0.5*sum(leonard_gi(:,:,gi)*mij) / sum(mij*mij)
 
2172
                    ! Constrain C_S to be between 0 and 0.04.
 
2173
                    les_coef_gi(gi) = min(max(les_coef_gi(gi),0.0), 0.04)
 
2174
                    ! Anisotropic tensor dynamic eddy viscosity m_ij = -2C_S|S1|.alpha^2.G1
 
2175
                    les_tensor_gi(:,:,gi) = 2*alpha**2*les_coef_gi(gi)*strain_mod(gi)*f_tensor(:,:,gi)
 
2176
                  end do
 
2177
            end select
2181
2178
 
2182
2179
            ! Assemble diagnostic fields
2183
2180
            call les_assemble_diagnostic_fields(state, nu, ele, detwei, &
2184
 
                 mesh_size_gi, strain_gi, t_strain_gi, les_tensor_gi, &
2185
 
                 have_eddy_visc, have_strain, have_filtered_strain, have_filter_width)
 
2181
                 mesh_size_gi, les_tensor_gi, les_coef_gi, &
 
2182
                 have_eddy_visc, have_filter_width, have_coeff)
2186
2183
 
2187
2184
         else
2188
2185
            FLAbort("Unknown LES model")