← Back to branch summary

~fluidity-core/fluidity/bug_1183080

« back to all changes in this revision

Viewing changes to assemble/LES.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/LES.F90
48
48
 
49
49
contains
50
50
 
51
 
  subroutine les_init_diagnostic_fields(state, have_eddy_visc, have_strain, have_filtered_strain, have_filter_width)
 
51
  subroutine les_init_diagnostic_fields(state, have_eddy_visc, have_filter_width, have_coeff)
52
52
 
53
53
    ! Arguments
54
 
    type(state_type), intent(inout) :: state
55
 
    logical, intent(in) :: have_eddy_visc, have_strain, have_filtered_strain, have_filter_width
 
54
    type(state_type), intent(inout)             :: state
 
55
    logical, intent(in)                         :: have_eddy_visc, have_filter_width, have_coeff
56
56
    
57
57
    ! Local variables
58
 
    logical, dimension(4) :: have_diagnostic_field
59
 
    character(len=FIELD_NAME_LEN), dimension(4) :: diagnostic_field_names
60
 
    type(tensor_field), pointer :: field
61
 
    integer :: i
 
58
    logical, dimension(2)                       :: have_diagnostic_tfield
 
59
    logical, dimension(1)                       :: have_diagnostic_sfield
 
60
    character(len=FIELD_NAME_LEN), dimension(2) :: diagnostic_tfield_names
 
61
    character(len=FIELD_NAME_LEN), dimension(1) :: diagnostic_sfield_names
 
62
    type(tensor_field), pointer                 :: tfield
 
63
    type(scalar_field), pointer                 :: sfield
 
64
    integer                                     :: i
62
65
 
63
66
    ewrite(2,*) "Initialising optional LES diagnostic fields"
64
67
    
65
 
    have_diagnostic_field = (/have_eddy_visc, have_strain, have_filtered_strain, have_filter_width/)
66
 
    diagnostic_field_names(1) = "EddyViscosity"
67
 
    diagnostic_field_names(2) = "StrainRate"
68
 
    diagnostic_field_names(3) = "FilteredStrainRate"
69
 
    diagnostic_field_names(4) = "FilterWidth"
 
68
    have_diagnostic_tfield = (/have_eddy_visc, have_filter_width/)
 
69
    diagnostic_tfield_names(1) = "EddyViscosity"
 
70
    diagnostic_tfield_names(2) = "FilterWidth"
70
71
    
71
 
    diagnostic_field_loop: do i = 1, size(diagnostic_field_names)
72
 
      if(have_diagnostic_field(i)) then
73
 
         field => extract_tensor_field(state, diagnostic_field_names(i))
74
 
         call zero(field)
75
 
      end if
76
 
    end do diagnostic_field_loop
 
72
    diagnostic_tfield_loop: do i = 1, size(diagnostic_tfield_names)
 
73
      if(have_diagnostic_tfield(i)) then
 
74
         tfield => extract_tensor_field(state, diagnostic_tfield_names(i))
 
75
         call zero(tfield)
 
76
      end if
 
77
    end do diagnostic_tfield_loop
 
78
 
 
79
    have_diagnostic_sfield = (/have_coeff/)
 
80
    diagnostic_sfield_names(1) = "SmagorinskyCoefficient"
 
81
 
 
82
    diagnostic_sfield_loop: do i = 1, size(diagnostic_sfield_names)
 
83
      if(have_diagnostic_sfield(i)) then
 
84
         sfield => extract_scalar_field(state, diagnostic_sfield_names(i))
 
85
         call zero(sfield)
 
86
      end if
 
87
    end do diagnostic_sfield_loop
77
88
 
78
89
  end subroutine les_init_diagnostic_fields
79
90
 
80
91
  subroutine les_assemble_diagnostic_fields(state, nu, ele, detwei, &
81
 
                 mesh_size_gi, strain_gi, t_strain_gi, les_tensor_gi, &
82
 
                 have_eddy_visc, have_strain, have_filtered_strain, have_filter_width)
 
92
                 mesh_size_gi,les_tensor_gi, les_coef_gi, &
 
93
                 have_eddy_visc, have_filter_width, have_coeff)
83
94
 
84
95
    ! Arguments
85
96
    type(state_type), intent(inout)                             :: state
86
97
    type(vector_field), intent(in)                              :: nu
87
98
    integer, intent(in)                                         :: ele
88
 
    real, dimension(ele_ngi(nu,ele)), intent(in)                :: detwei
89
 
    real, dimension(nu%dim,nu%dim,ele_ngi(nu,ele)),intent(in)   :: strain_gi, t_strain_gi, mesh_size_gi, les_tensor_gi
90
 
    logical, intent(in) :: have_eddy_visc, have_strain, have_filtered_strain, have_filter_width
 
99
    real, dimension(ele_ngi(nu,ele)), intent(in)                :: les_coef_gi, detwei
 
100
    real, dimension(nu%dim,nu%dim,ele_ngi(nu,ele)),intent(in)   :: mesh_size_gi, les_tensor_gi
 
101
    logical, intent(in) :: have_eddy_visc, have_filter_width, have_coeff
91
102
    
92
103
    ! Local variables
93
 
    type(tensor_field), pointer                                 :: field
 
104
    type(tensor_field), pointer                                 :: tfield
 
105
    type(scalar_field), pointer                                 :: sfield
94
106
    real, dimension(nu%dim,nu%dim,ele_loc(nu,ele))              :: tensor_loc
 
107
    real, dimension(ele_loc(nu,ele))                            :: scalar_loc
95
108
 
96
 
    ! Eddy viscosity field m_ij
 
109
    ! Eddy viscosity
97
110
    if(have_eddy_visc) then
98
 
      field => extract_tensor_field(state, "EddyViscosity")
 
111
      tfield => extract_tensor_field(state, "EddyViscosity")
99
112
      tensor_loc=shape_tensor_rhs(ele_shape(nu, ele), les_tensor_gi, detwei)
100
 
      call addto(field, ele_nodes(nu, ele), tensor_loc)
101
 
    end if
102
 
 
103
 
    ! Strain rate field S1
104
 
    if(have_strain) then
105
 
      field => extract_tensor_field(state, "StrainRate")
106
 
      tensor_loc=shape_tensor_rhs(ele_shape(nu, ele), strain_gi, detwei)
107
 
      call addto(field, ele_nodes(nu, ele), tensor_loc)
108
 
    end if
109
 
 
110
 
    ! Filtered strain rate field S2
111
 
    if(have_filtered_strain) then
112
 
      field => extract_tensor_field(state, "FilteredStrainRate")
113
 
      tensor_loc=shape_tensor_rhs(ele_shape(nu, ele), t_strain_gi, detwei)
114
 
      call addto(field, ele_nodes(nu, ele), tensor_loc)
 
113
      call addto(tfield, ele_nodes(nu, ele), tensor_loc)
115
114
    end if
116
115
 
117
116
    ! Filter width
118
117
    if(have_filter_width) then
119
 
      field => extract_tensor_field(state, "FilterWidth")
 
118
      tfield => extract_tensor_field(state, "FilterWidth")
120
119
      tensor_loc=shape_tensor_rhs(ele_shape(nu, ele), mesh_size_gi, detwei)
121
 
      call addto(field, ele_nodes(nu, ele), tensor_loc)
 
120
      call addto(tfield, ele_nodes(nu, ele), tensor_loc)
 
121
    end if
 
122
 
 
123
    ! Smagorinsky Coefficient
 
124
    if(have_coeff) then
 
125
      sfield => extract_scalar_field(state, "SmagorinskyCoefficient")
 
126
      scalar_loc=shape_rhs(ele_shape(nu, ele), les_coef_gi*detwei)
 
127
      call addto(sfield, ele_nodes(nu, ele), scalar_loc)
122
128
    end if
123
129
 
124
130
  end subroutine les_assemble_diagnostic_fields
125
131
 
126
 
  subroutine les_solve_diagnostic_fields(state, have_eddy_visc, have_strain, have_filtered_strain, have_filter_width)
 
132
  subroutine les_solve_diagnostic_fields(state, have_eddy_visc, have_filter_width, have_coeff)
127
133
 
128
134
    ! Arguments
129
135
    type(state_type), intent(inout) :: state
130
 
    logical, intent(in) :: have_eddy_visc, have_strain, have_filtered_strain, have_filter_width
 
136
    logical, intent(in) :: have_eddy_visc, have_filter_width, have_coeff
131
137
    
132
138
    ! Local variables
133
 
    logical, dimension(4) :: have_diagnostic_field
134
 
    character(len=FIELD_NAME_LEN), dimension(4) :: diagnostic_field_names
135
 
    type(tensor_field), pointer :: field
136
 
    integer :: i
137
 
    type(vector_field), pointer :: u
138
 
    type(csr_matrix), pointer :: mass_matrix
139
 
    type(scalar_field), pointer :: lumped_mass
140
 
    type(scalar_field) :: inv_lumped_mass
141
 
    logical :: lump_mass = .false.
142
 
    logical :: use_submesh = .false.
 
139
    logical, dimension(2)                       :: have_diagnostic_tfield
 
140
    logical, dimension(1)                       :: have_diagnostic_sfield
 
141
    character(len=FIELD_NAME_LEN), dimension(2) :: diagnostic_tfield_names
 
142
    character(len=FIELD_NAME_LEN), dimension(1) :: diagnostic_sfield_names
 
143
    type(tensor_field), pointer                 :: tfield
 
144
    type(scalar_field), pointer                 :: sfield
 
145
    integer                                     :: i
 
146
    type(vector_field), pointer                 :: u
 
147
    type(csr_matrix), pointer                   :: mass_matrix
 
148
    type(scalar_field), pointer                 :: lumped_mass
 
149
    type(scalar_field)                          :: inv_lumped_mass
 
150
    logical                                     :: lump_mass = .false.
 
151
    logical                                     :: use_submesh = .false.
143
152
    
144
153
    ewrite(2,*) "Solving for optional LES diagnostic fields"
145
154
        
146
155
    u => extract_vector_field(state, "Velocity")
147
156
    
148
 
    have_diagnostic_field = (/have_eddy_visc, have_strain, have_filtered_strain, have_filter_width/)
149
 
    diagnostic_field_names(1) = "EddyViscosity"
150
 
    diagnostic_field_names(2) = "StrainRate"
151
 
    diagnostic_field_names(3) = "FilteredStrainRate"
152
 
    diagnostic_field_names(4) = "FilterWidth"
153
 
    
154
 
    diagnostic_field_loop: do i = 1, size(diagnostic_field_names)
155
 
      if(have_diagnostic_field(i)) then
156
 
         field => extract_tensor_field(state, diagnostic_field_names(i))
 
157
    have_diagnostic_tfield = (/have_eddy_visc, have_filter_width/)
 
158
    diagnostic_tfield_names(1) = "EddyViscosity"
 
159
    diagnostic_tfield_names(2) = "FilterWidth"
 
160
    
 
161
    diagnostic_tfield_loop: do i = 1, size(diagnostic_tfield_names)
 
162
      if(have_diagnostic_tfield(i)) then
 
163
         tfield => extract_tensor_field(state, diagnostic_tfield_names(i))
 
164
         lump_mass = have_option(trim(tfield%option_path)//"/diagnostic/mass_matrix"//&
 
165
            &"/use_lumped_mass_matrix")
 
166
         use_submesh = have_option(trim(tfield%option_path)//"/diagnostic/mass_matrix"//&
 
167
            &"/use_lumped_mass_matrix/use_submesh") ! For P2 meshes.
 
168
            
 
169
         if(lump_mass) then
 
170
            if(use_submesh) then
 
171
               lumped_mass => get_lumped_mass_on_submesh(state, tfield%mesh)
 
172
            else
 
173
               lumped_mass => get_lumped_mass(state, tfield%mesh)
 
174
            end if
 
175
            call allocate(inv_lumped_mass, tfield%mesh)
 
176
            call invert(lumped_mass, inv_lumped_mass)
 
177
            call scale(tfield, inv_lumped_mass)
 
178
            call deallocate(inv_lumped_mass)
 
179
         else
 
180
            mass_matrix => get_mass_matrix(state, tfield%mesh)
 
181
            call petsc_solve(tfield, mass_matrix, tfield, option_path=u%option_path)
 
182
         end if
 
183
      end if
 
184
    end do diagnostic_tfield_loop
 
185
    
 
186
    have_diagnostic_sfield = (/have_coeff/)
 
187
    diagnostic_sfield_names(1) = "SmagorinskyCoefficient"
 
188
    
 
189
    diagnostic_sfield_loop: do i = 1, size(diagnostic_sfield_names)
 
190
      if(have_diagnostic_sfield(i)) then
 
191
         sfield => extract_scalar_field(state, diagnostic_sfield_names(i))
 
192
         lump_mass = have_option(trim(sfield%option_path)//"/diagnostic/mass_matrix"//&
 
193
            &"/use_lumped_mass_matrix")
 
194
         use_submesh = have_option(trim(sfield%option_path)//"/diagnostic/mass_matrix"//&
 
195
            &"/use_lumped_mass_matrix/use_submesh") ! For P2 meshes.
 
196
            
 
197
         if(lump_mass) then
 
198
            if(use_submesh) then
 
199
               lumped_mass => get_lumped_mass_on_submesh(state, sfield%mesh)
 
200
            else
 
201
               lumped_mass => get_lumped_mass(state, sfield%mesh)
 
202
            end if
 
203
            call allocate(inv_lumped_mass, sfield%mesh)
 
204
            call invert(lumped_mass, inv_lumped_mass)
 
205
            call scale(sfield, inv_lumped_mass)
 
206
            call deallocate(inv_lumped_mass)
 
207
         else
 
208
            mass_matrix => get_mass_matrix(state, sfield%mesh)
 
209
            call petsc_solve(sfield, mass_matrix, sfield, option_path=u%option_path)
 
210
         end if
 
211
      end if
 
212
    end do diagnostic_sfield_loop
157
213
 
158
 
         lump_mass = have_option(trim(field%option_path)//"/diagnostic/mass_matrix"//&
159
 
            &"/use_lumped_mass_matrix")
160
 
         use_submesh = have_option(trim(field%option_path)//"/diagnostic/mass_matrix"//&
161
 
            &"/use_lumped_mass_matrix/use_submesh") ! For P2 meshes.
162
 
            
163
 
         if(lump_mass) then
164
 
            if(use_submesh) then
165
 
               lumped_mass => get_lumped_mass_on_submesh(state, field%mesh)
166
 
            else
167
 
               lumped_mass => get_lumped_mass(state, field%mesh)
168
 
            end if
169
 
            call allocate(inv_lumped_mass, field%mesh)
170
 
            call invert(lumped_mass, inv_lumped_mass)
171
 
            call scale(field, inv_lumped_mass)
172
 
            call deallocate(inv_lumped_mass)
173
 
         else
174
 
            mass_matrix => get_mass_matrix(state, field%mesh)
175
 
            call petsc_solve(field, mass_matrix, field, option_path=u%option_path)
176
 
         end if
177
 
         
178
 
      end if
179
 
    end do diagnostic_field_loop
180
 
    
181
214
  end subroutine les_solve_diagnostic_fields
182
215
  
183
 
  subroutine leonard_tensor(nu, positions, tnu, leonard, alpha, path)
 
216
  subroutine leonard_tensor(nu, positions, fnu, tnu, leonard, strainprod, alpha, gamma, path)
184
217
 
185
218
    ! Unfiltered velocity
186
 
    type(vector_field), pointer               :: nu
187
 
    type(vector_field), intent(in)            :: positions
188
 
    ! Filtered velocity
189
 
    type(vector_field), pointer               :: tnu
190
 
    ! Leonard tensor field
191
 
    type(tensor_field), pointer               :: leonard
192
 
    ! Scale factor: test filter/mesh size
193
 
    real, intent(in)                          :: alpha
194
 
    character(len=OPTION_PATH_LEN), intent(in):: path
 
219
    type(vector_field), pointer                           :: nu
 
220
    type(vector_field), intent(in)                        :: positions
 
221
    ! Filtered velocities
 
222
    type(vector_field), pointer                           :: fnu, tnu
 
223
    ! Leonard tensor and strain product
 
224
    type(tensor_field), pointer                           :: leonard, strainprod
 
225
    ! Scale factors
 
226
    real, intent(in)                                      :: alpha, gamma
 
227
    character(len=OPTION_PATH_LEN), intent(in)            :: path
195
228
    ! Local quantities
196
 
    type(tensor_field), pointer               :: ui_uj, tui_tuj
197
 
    character(len=OPTION_PATH_LEN)            :: lpath
198
 
    integer                                   :: i
199
 
    real, dimension(:), allocatable           :: u_loc
200
 
    real, dimension(:,:), allocatable         :: t_loc
 
229
    type(tensor_field), pointer                           :: ui_uj, tui_tuj
 
230
    character(len=OPTION_PATH_LEN)                        :: lpath
 
231
    integer                                               :: i, ele, gi
 
232
    real, dimension(:), allocatable                       :: u_loc
 
233
    real, dimension(:,:), allocatable                     :: t_loc
 
234
    real, dimension(ele_loc(nu,1), ele_ngi(nu,1), nu%dim) :: du_t
 
235
    real, dimension(ele_ngi(nu,1))                        :: detwei
 
236
    real, dimension(nu%dim, nu%dim, ele_ngi(nu,1))        :: strain_gi, strain_prod_gi
 
237
    type(element_type)                                    :: shape_nu
201
238
 
202
239
    ! Path is to level above solver options
203
240
    lpath = (trim(path)//"/dynamic_les")
204
241
    ewrite(2,*) "filter factor alpha: ", alpha
205
 
 
206
 
    ewrite_minmax(nu)
207
 
    ewrite_minmax(tnu)
208
 
 
209
 
    call anisotropic_smooth_vector(nu, positions, tnu, alpha, lpath)
210
 
 
211
 
    ewrite_minmax(nu)
 
242
    ewrite(2,*) "filter factor gamma: ", gamma
 
243
 
 
244
    ! First filter operator returns u^f:
 
245
    call anisotropic_smooth_vector(nu, positions, fnu, alpha, lpath)
 
246
    ! Test filter operator needs the ratio of test filter to mesh size and returns u^ft:
 
247
    call anisotropic_smooth_vector(fnu, positions, tnu, alpha*gamma, lpath)
 
248
    ewrite_minmax(nu)
 
249
    ewrite_minmax(fnu)
212
250
    ewrite_minmax(tnu)
213
251
 
214
252
    ! Velocity products (ui*uj)
223
261
 
224
262
    ! Get cross products of velocities
225
263
    do i=1, node_count(nu)
226
 
      ! Mesh filter ^r
227
 
      u_loc = node_val(nu,i)
 
264
      u_loc = node_val(fnu,i)
228
265
      t_loc = outer_product(u_loc, u_loc)
229
266
      call set( ui_uj, i, t_loc )
230
 
      ! Test filter ^t
231
267
      u_loc = node_val(tnu,i)
232
 
      ! Calculate (test-filtered velocity) products: (ui^rt*uj^rt)
 
268
      ! Calculate (test-filtered velocity) products: (ui^ft*uj^ft)
233
269
      t_loc = outer_product(u_loc, u_loc)
234
270
      call set( tui_tuj, i, t_loc )
235
271
    end do
236
272
 
237
 
    ! Calculate test-filtered (velocity products): (ui^r*uj^r)^t
238
 
    call anisotropic_smooth_tensor(ui_uj, positions, leonard, alpha, lpath)
 
273
    ! Calculate test-filtered (velocity products): (ui^f*uj^f)^t
 
274
    call anisotropic_smooth_tensor(ui_uj, positions, leonard, alpha*gamma, lpath)
239
275
 
240
276
    ! Leonard tensor field
241
277
    call addto( leonard, tui_tuj, -1.0 )
242
278
 
 
279
    ! Zero tensor field for reuse in strain product assembly
 
280
    call zero(ui_uj)
 
281
 
 
282
    do i=1, element_count(nu)
 
283
      shape_nu = ele_shape(nu, i)
 
284
      ! Assuming no FE stabilisation is used with LES so we can use velocity shape.
 
285
      call transform_to_physical(positions, i, shape_nu, dshape=du_t, detwei=detwei)
 
286
      ! Strain rate of first filtered velocity S1^f
 
287
      strain_gi = les_strain_rate(du_t, ele_val(fnu, i))
 
288
      do gi=1, ele_ngi(nu, ele)
 
289
        ! Strain product = strain modulus*strain rate: |S1^f|S1^f
 
290
        strain_prod_gi(:,:,gi) = sqrt(2*sum(strain_gi(:,:,gi)*strain_gi(:,:,gi))) * strain_gi(:,:,gi)
 
291
      end do
 
292
      ! Assemble local tensor field
 
293
      call addto(ui_uj, ele_nodes(nu,i), shape_tensor_rhs(ele_shape(nu,i), strain_prod_gi, detwei))
 
294
    end do
 
295
 
 
296
    ! Filter strain product with test filter: (|S1^f|S1^f)^t
 
297
    call anisotropic_smooth_tensor(ui_uj, positions, strainprod, alpha*gamma, lpath)
 
298
 
243
299
    ! Deallocates
244
300
    deallocate(u_loc, t_loc)
245
301
    call deallocate(ui_uj)