← Back to branch summary

~airpollution/fluidity/fluidity_airpollution

« back to all changes in this revision

Viewing changes to assemble/Assemble_CMC.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/Assemble_CMC.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
 
 
28
#include "fdebug.h"
 
29
 
 
30
module assemble_CMC
 
31
 
 
32
  use fldebug
 
33
  use state_module
 
34
  use sparse_tools
 
35
  use sparse_tools_petsc
 
36
  use spud
 
37
  use fields
 
38
  use fields_base
 
39
  use fefields
 
40
  use sparse_matrices_fields
 
41
  use field_options
 
42
  use global_parameters, only: OPTION_PATH_LEN
 
43
  use linked_lists
 
44
  
 
45
  implicit none 
 
46
 
 
47
  private
 
48
  
 
49
  public :: assemble_cmc_dg, assemble_masslumped_cmc, &
 
50
            assemble_masslumped_ctm, repair_stiff_nodes, zero_stiff_nodes, &
 
51
            assemble_diagonal_schur, assemble_scaled_pressure_mass_matrix
 
52
 
 
53
contains
 
54
 
 
55
    subroutine assemble_cmc_dg(CMC, CTP, CT, inverse_mass)
 
56
      !!< Assemble the pressure matrix C^T M^{-1} C for a DG mesh.
 
57
      !!< This currently does not support rotations.
 
58
      type(csr_matrix), intent(inout) :: CMC
 
59
      type(block_csr_matrix), intent(in) :: CTP, CT
 
60
      type(block_csr_matrix), intent(in):: inverse_mass
 
61
 
 
62
      type(csr_matrix) :: CM ! Temporary half way matrix
 
63
      type(csr_matrix) :: CMC_tmp ! Temporary accumulator.
 
64
      type(csr_matrix) :: CT_block ! The current CT dimension
 
65
      type(csr_matrix) :: CTP_block ! The current CTP dimension
 
66
 
 
67
      integer :: dim, dim2
 
68
 
 
69
      call zero(CMC)
 
70
      
 
71
      if(inverse_mass%diagonal) then
 
72
 
 
73
         do dim=1, blocks(CT,2)
 
74
            
 
75
            ! This is a borrowed reference so no need to deallocate.
 
76
            CT_block=block(CT, 1, dim)
 
77
            CTP_block=block(CTP, 1, dim)
 
78
            
 
79
            CM=matmul_T(CTP_block, block(inverse_mass,dim,dim), CT%sparsity)
 
80
            
 
81
            CMC_tmp=matmul_T(CM, CT_block, model=CMC%sparsity)
 
82
            
 
83
            call addto(CMC, CMC_tmp)
 
84
            
 
85
            call deallocate(CM)
 
86
            call deallocate(CMC_tmp)
 
87
 
 
88
         end do
 
89
 
 
90
      else
 
91
 
 
92
         do dim=1, blocks(CT,2)
 
93
            do dim2 = 1, blocks(CT,2)
 
94
               
 
95
               ! This is a borrowed reference so no need to deallocate.
 
96
               CT_block=block(CT, 1, dim)
 
97
               CTP_block=block(CTP, 1, dim2)
 
98
               
 
99
               CM=matmul_T(CTP_block, block(inverse_mass,dim2,dim), CT%sparsity)
 
100
               
 
101
               CMC_tmp=matmul_T(CM, CT_block, model=CMC%sparsity)
 
102
               
 
103
               call addto(CMC, CMC_tmp)
 
104
               
 
105
               call deallocate(CM)
 
106
               call deallocate(CMC_tmp)
 
107
            end do
 
108
         end do
 
109
 
 
110
      end if
 
111
 
 
112
      ewrite_minmax(cmc)
 
113
 
 
114
    end subroutine assemble_cmc_dg
 
115
 
 
116
    subroutine assemble_masslumped_cmc(cmc_m, ctp_m, inverse_masslump, ct_m)
 
117
      !!< Assemble the pressure matrix C_P^T M_l^{-1} C.
 
118
      !!< This currently does not support rotations.
 
119
 
 
120
      ! this subroutine is designed to start to replace cmc_wrapper and getcmc
 
121
 
 
122
      type(csr_matrix), intent(inout) :: cmc_m
 
123
      type(block_csr_matrix), intent(in) :: ctp_m, ct_m
 
124
      type(vector_field), intent(in) :: inverse_masslump
 
125
 
 
126
      ewrite(1,*) 'Entering assemble_masslumped_cmc'
 
127
 
 
128
      call mult_div_vector_div_T(cmc_m, ctp_m, inverse_masslump, ct_m)
 
129
      
 
130
      ewrite_minmax(cmc_m)
 
131
 
 
132
    end subroutine assemble_masslumped_cmc
 
133
 
 
134
    subroutine assemble_diagonal_schur(schur_diagonal_matrix,u,inner_m,ctp_m,ct_m)
 
135
      !!< Assemble the matrix C_P^T * [(Big_m)_diagonal]^-1 * C. 
 
136
      !!< This is used as a preconditioner for the full projection solve
 
137
      !!< when using the full momentum matrix.
 
138
 
 
139
      ! Fluidity velocity vector:
 
140
      type(vector_field), intent(in):: u
 
141
      ! Divergence matrices:
 
142
      type(block_csr_matrix), intent(in) :: ctp_m, ct_m
 
143
      ! Inner matrix of Schur complement:
 
144
      type(petsc_csr_matrix), intent(inout) :: inner_m   
 
145
      ! Product matrix:
 
146
      type(csr_matrix), intent(inout):: schur_diagonal_matrix
 
147
      ! Diagonal of inner_m - required to set up preconditioner matrix for Stokes problems
 
148
      ! (i.e. DiagonalSchurComplement):
 
149
      type(vector_field) :: inner_m_diagonal   
 
150
 
 
151
      integer :: i
 
152
 
 
153
      ewrite(1,*) 'Entering assemble_diagonal_schur'
 
154
 
 
155
      call zero(schur_diagonal_matrix)
 
156
      call allocate(inner_m_diagonal,u%dim,u%mesh,"Diagonal_inner_m")
 
157
      call zero(inner_m_diagonal)
 
158
      call extract_diagonal(inner_m,inner_m_diagonal)
 
159
 
 
160
      ewrite_minmax(inner_m_diagonal)
 
161
      if(any(inner_m_diagonal%val < 0)) then
 
162
        ewrite(-1,*) 'Inner_m_diagonal has negative values'
 
163
        FLExit("Negative values in the diagonal schur complement preconditioner")
 
164
 
 
165
      end if
 
166
 
 
167
      call mult_div_invvector_div_T(schur_diagonal_matrix, ctp_m, inner_m_diagonal, ct_m)
 
168
      ewrite_minmax(schur_diagonal_matrix)
 
169
      call deallocate(inner_m_diagonal)
 
170
 
 
171
    end subroutine assemble_diagonal_schur
 
172
 
 
173
    subroutine assemble_scaled_pressure_mass_matrix(state, scaled_pressure_mass_matrix, p_mesh, dt)
 
174
 
 
175
      ! This routine assembles the scaled_pressure_mass_matrix at the
 
176
      ! quadrature points. It is scaled by the inverse of viscosity.
 
177
 
 
178
      type(state_type), intent(in) :: state   
 
179
      ! Scaled pressure mass matrix - already allocated in Momentum_Eq:
 
180
      type(csr_matrix), intent(inout) :: scaled_pressure_mass_matrix
 
181
 
 
182
      ! Pressure mesh (for free surface this should be the extended pressure mesh)
 
183
      type(mesh_type), intent(in) :: p_mesh
 
184
      real, intent(in) :: dt
 
185
 
 
186
      ! Viscosity tensor:
 
187
      type(tensor_field), pointer :: viscosity     
 
188
      ! Viscosity component:
 
189
      type(scalar_field) :: viscosity_component
 
190
      ! Positions:
 
191
      type(vector_field), pointer :: positions
 
192
 
 
193
      ! Relevant declerations for mass matrix calculation:
 
194
      integer :: ele
 
195
      real, dimension(:), allocatable :: detwei
 
196
      type(element_type), pointer :: p_shape
 
197
      real, dimension(:,:), allocatable :: mass_matrix
 
198
      real, dimension(:), allocatable :: mu_gi
 
199
 
 
200
      ewrite(1,*) 'Entering assemble_scaled_pressure_mass_matrix'    
 
201
 
 
202
      ! Positions:
 
203
      positions => extract_vector_field(state, "Coordinate")
 
204
 
 
205
      ! Extract viscosity tensor from state:
 
206
      viscosity => extract_tensor_field(state,'Viscosity')
 
207
 
 
208
      ! Extract first component of viscosity tensor from full tensor:
 
209
      viscosity_component = extract_scalar_field(viscosity,1,1)
 
210
 
 
211
      ! Initialise and assemble scaled pressure mass matrix:
 
212
      allocate(detwei(ele_ngi(p_mesh, 1)), &
 
213
               mass_matrix(ele_loc(p_mesh, 1), ele_loc(p_mesh, 1)), &
 
214
               mu_gi(ele_ngi(viscosity_component, 1)))
 
215
 
 
216
      call zero(scaled_pressure_mass_matrix)
 
217
      do ele = 1, ele_count(p_mesh)
 
218
        p_shape => ele_shape(p_mesh, ele)
 
219
        mu_gi = ele_val_at_quad(viscosity_component, ele)
 
220
        call transform_to_physical(positions, ele, detwei=detwei)
 
221
        mass_matrix = shape_shape(p_shape, p_shape, detwei/(mu_gi*dt))
 
222
        call addto(scaled_pressure_mass_matrix, ele_nodes(p_mesh, ele),&
 
223
             ele_nodes(p_mesh, ele), mass_matrix)
 
224
      end do
 
225
 
 
226
      ewrite_minmax(scaled_pressure_mass_matrix)
 
227
 
 
228
      deallocate(detwei, mass_matrix, mu_gi)
 
229
 
 
230
    end subroutine assemble_scaled_pressure_mass_matrix
 
231
      
 
232
    subroutine repair_stiff_nodes(cmc_m, stiff_nodes_list)
 
233
    
 
234
      type(csr_matrix), intent(inout) :: cmc_m
 
235
      type(ilist), intent(inout) :: stiff_nodes_list
 
236
      
 
237
      integer :: row
 
238
      real, pointer :: row_diag
 
239
      integer, dimension(:), pointer :: row_m
 
240
      real, dimension(:), pointer :: row_val
 
241
      real :: tolerance
 
242
      
 
243
      ewrite(1,*) 'in repair_stiff_nodes()'
 
244
      
 
245
      tolerance = maxval(cmc_m%val)*epsilon(0.0)
 
246
      ewrite(2,*) 'tolerance = ', tolerance
 
247
      
 
248
      call flush_list(stiff_nodes_list)
 
249
      
 
250
      do row = 1, size(cmc_m, 1)
 
251
        row_diag=>diag_val_ptr(cmc_m, row)
 
252
        row_m=>row_m_ptr(cmc_m, row)
 
253
        row_val=>row_val_ptr(cmc_m, row)
 
254
        if(row_diag<tolerance) then
 
255
          ewrite(2,*) 'before, row, row_diag, sum(row_val) = ', row, row_diag, sum(abs(row_val))
 
256
          where(cmc_m%sparsity%colm==row) cmc_m%val = 0.0
 
257
          call zero_row(cmc_m, row)
 
258
          call addto_diag(cmc_m, row, 1.0)
 
259
          call insert(stiff_nodes_list, row)
 
260
        end if
 
261
      end do
 
262
      
 
263
      call print_list(stiff_nodes_list, 2)
 
264
 
 
265
    end subroutine repair_stiff_nodes
 
266
      
 
267
    subroutine zero_stiff_nodes(rhs, stiff_nodes_list)
 
268
    
 
269
      type(scalar_field), intent(inout) :: rhs
 
270
      type(ilist), intent(in) :: stiff_nodes_list
 
271
      
 
272
      real, dimension(:), pointer :: row_val
 
273
      
 
274
      ewrite(1,*) 'in zero_stiff_nodes()'
 
275
      
 
276
      if(stiff_nodes_list%length>0) then
 
277
        ewrite(2,*) 'before node_val = ', node_val(rhs, list2vector(stiff_nodes_list))
 
278
        call set(rhs, list2vector(stiff_nodes_list), spread(0.0, 1, stiff_nodes_list%length))
 
279
      end if
 
280
 
 
281
    end subroutine zero_stiff_nodes
 
282
    
 
283
    subroutine assemble_masslumped_ctm(ctm_m, ctp_m, masslump)
 
284
      !!< Assemble the matrix C_P^T M_l^{-1}
 
285
      !!< This currently does not support rotations.
 
286
 
 
287
      type(block_csr_matrix), intent(in) :: ctp_m
 
288
      type(scalar_field), intent(in) :: masslump
 
289
      type(block_csr_matrix), intent(inout) :: ctm_m
 
290
 
 
291
      type(csr_matrix) :: lctm_m_block
 
292
 
 
293
      integer :: dim, row
 
294
      real, dimension(:), pointer :: row_val
 
295
      integer, dimension(:), pointer :: row_indices
 
296
 
 
297
      ewrite(1,*) 'Entering assemble_masslumped_ctm'
 
298
 
 
299
      call zero(ctm_m)
 
300
 
 
301
      do dim = 1, ctm_m%blocks(2)
 
302
 
 
303
        lctm_m_block = block(ctm_m, 1, dim)
 
304
 
 
305
        do row = 1, size(ctp_m, 1)
 
306
          row_indices=>row_m_ptr(ctp_m, row)
 
307
          row_val=>row_val_ptr(ctp_m, 1, dim, row)
 
308
          call set(lctm_m_block, (/row/), row_indices, &
 
309
                  spread((row_val/node_val(masslump, row_indices)), 1, 1))
 
310
        end do
 
311
 
 
312
      end do
 
313
 
 
314
      ewrite_minmax(ctm_m)
 
315
 
 
316
    end subroutine assemble_masslumped_ctm
 
317
    
 
318
end module assemble_cmc