|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
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 |
!
|
|
|
2392
by tmb1
Changing isntances of Chris' person email address to the group |
11 |
! amcgsoftware@imperial.ac.uk
|
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
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 differential_operator_diagnostics |
|
31 |
||
|
1562
by maddison
Applying the usual Intel module bug workaround |
32 |
! these 5 need to be on top and in this order, so as not to confuse silly old intel compiler
|
33 |
use quadrature |
|
34 |
use elements |
|
35 |
use sparse_tools |
|
36 |
use fields |
|
37 |
use state_module |
|
38 |
!
|
|
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
39 |
use diagnostic_source_fields |
|
773
by maddison
Adding a finite_element_divergence diagnostic algorithm. Bugfix in div algorithm and removing some invalid options from the shallow water schema |
40 |
use divergence_matrix_cg |
|
1093
by maddison
Applying usual tricks to make the divergence matrix caching work with mesh movement |
41 |
use eventcounter |
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
42 |
use field_derivatives |
43 |
use field_options |
|
44 |
use fldebug |
|
|
825
by maddison
Adding finite_element_transpose diagnostic algorithm |
45 |
use geostrophic_pressure |
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
46 |
use global_parameters, only : OPTION_PATH_LEN |
|
773
by maddison
Adding a finite_element_divergence diagnostic algorithm. Bugfix in div algorithm and removing some invalid options from the shallow water schema |
47 |
use solvers |
48 |
use sparse_matrices_fields |
|
49 |
use sparsity_patterns_meshes |
|
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
50 |
use spud |
51 |
use state_fields_module |
|
52 |
||
53 |
implicit none
|
|
54 |
|
|
55 |
private
|
|
56 |
|
|
|
1488
by maddison
Perp diagnostic algorithm |
57 |
public :: calculate_grad, calculate_div, calculate_curl, calculate_perp, & |
58 |
& calculate_curl_2d, calculate_finite_element_divergence, & |
|
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
59 |
& calculate_finite_element_divergence_transpose, & |
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
60 |
& calculate_scalar_advection, calculate_scalar_laplacian, & |
61 |
& calculate_vector_advection, calculate_vector_laplacian |
|
62 |
||
63 |
contains
|
|
64 |
||
65 |
subroutine calculate_grad(state, v_field) |
|
66 |
type(state_type), intent(in) :: state |
|
67 |
type(vector_field), intent(inout) :: v_field |
|
68 |
||
69 |
type(scalar_field), pointer :: source_field |
|
70 |
type(vector_field), pointer :: positions |
|
71 |
||
72 |
source_field => scalar_source_field(state, v_field) |
|
73 |
||
74 |
positions => extract_vector_field(state, "Coordinate") |
|
75 |
||
|
1997
by skramer
Taking the derivative of a discontinuous field is not supported in most diagnostic algorithms. Checks for this however were missing and fluidity would happily provide you with the wrong answer. The following diagnostic algorithms do not work for discontinuous source fields and now give an appropriate error message: |
76 |
call check_source_mesh_derivative(source_field, "grad") |
77 |
||
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
78 |
call grad(source_field, positions, v_field) |
79 |
||
80 |
end subroutine calculate_grad |
|
81 |
||
82 |
subroutine calculate_div(state, s_field) |
|
83 |
type(state_type), intent(in) :: state |
|
84 |
type(scalar_field), intent(inout) :: s_field |
|
85 |
||
86 |
type(vector_field), pointer :: positions, source_field |
|
87 |
||
88 |
source_field => vector_source_field(state, s_field) |
|
89 |
||
90 |
positions => extract_vector_field(state, "Coordinate") |
|
91 |
||
|
1997
by skramer
Taking the derivative of a discontinuous field is not supported in most diagnostic algorithms. Checks for this however were missing and fluidity would happily provide you with the wrong answer. The following diagnostic algorithms do not work for discontinuous source fields and now give an appropriate error message: |
92 |
call check_source_mesh_derivative(source_field, "div") |
93 |
||
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
94 |
call div(source_field, positions, s_field) |
95 |
||
96 |
end subroutine calculate_div |
|
|
773
by maddison
Adding a finite_element_divergence diagnostic algorithm. Bugfix in div algorithm and removing some invalid options from the shallow water schema |
97 |
|
98 |
subroutine calculate_finite_element_divergence(state, s_field) |
|
99 |
type(state_type), intent(inout) :: state |
|
100 |
type(scalar_field), intent(inout) :: s_field |
|
101 |
||
102 |
character(len = OPTION_PATH_LEN) :: path |
|
|
1092
by maddison
Add some caching to the finite_element_divergence diagnostic |
103 |
type(block_csr_matrix) :: ct_m |
|
773
by maddison
Adding a finite_element_divergence diagnostic algorithm. Bugfix in div algorithm and removing some invalid options from the shallow water schema |
104 |
type(csr_sparsity), pointer :: divergence_sparsity |
105 |
type(csr_matrix), pointer :: mass |
|
106 |
type(scalar_field) :: ctfield, ct_rhs |
|
107 |
type(scalar_field), pointer :: masslump |
|
108 |
type(vector_field), pointer :: positions, source_field |
|
109 |
||
110 |
source_field => vector_source_field(state, s_field) |
|
111 |
path = trim(complete_field_path(s_field%option_path)) // "/algorithm" |
|
|
1092
by maddison
Add some caching to the finite_element_divergence diagnostic |
112 |
|
|
1093
by maddison
Applying usual tricks to make the divergence matrix caching work with mesh movement |
113 |
if(use_divergence_matrix_cache(state)) then |
114 |
ct_m = extract_block_csr_matrix(state, trim(s_field%name) // "DivergenceMatrix") |
|
115 |
call incref(ct_m) |
|
|
1092
by maddison
Add some caching to the finite_element_divergence diagnostic |
116 |
ct_rhs = extract_scalar_field(state, trim(s_field%name) // "DivergenceRHS") |
117 |
call incref(ct_rhs) |
|
118 |
else
|
|
119 |
positions => extract_vector_field(state, "Coordinate") |
|
120 |
||
121 |
divergence_sparsity => get_csr_sparsity_firstorder(state, s_field%mesh, source_field%mesh) |
|
122 |
call allocate(ct_m, divergence_sparsity, (/1, source_field%dim/), name = "DivergenceMatrix" ) |
|
123 |
call allocate(ct_rhs, s_field%mesh, name = "CTRHS") |
|
124 |
||
125 |
call assemble_divergence_matrix_cg(ct_m, state, ct_rhs = ct_rhs, & |
|
126 |
& test_mesh = s_field%mesh, field = source_field, & |
|
127 |
& option_path = path) |
|
128 |
call insert(state, ct_m, trim(s_field%name) // "DivergenceMatrix") |
|
129 |
call insert(state, ct_rhs, trim(s_field%name) // "DivergenceRHS") |
|
130 |
end if
|
|
|
773
by maddison
Adding a finite_element_divergence diagnostic algorithm. Bugfix in div algorithm and removing some invalid options from the shallow water schema |
131 |
|
|
1092
by maddison
Add some caching to the finite_element_divergence diagnostic |
132 |
call allocate(ctfield, s_field%mesh, name = "CTField") |
133 |
call mult(ctfield, ct_m, source_field) |
|
|
773
by maddison
Adding a finite_element_divergence diagnostic algorithm. Bugfix in div algorithm and removing some invalid options from the shallow water schema |
134 |
call addto(ctfield, ct_rhs, -1.0) |
135 |
call deallocate(ct_rhs) |
|
136 |
||
137 |
if(have_option(trim(path) // "/lump_mass")) then |
|
138 |
masslump => get_lumped_mass(state, s_field%mesh) |
|
139 |
s_field%val = ctfield%val / masslump%val |
|
140 |
else
|
|
141 |
mass => get_mass_matrix(state, s_field%mesh) |
|
142 |
call petsc_solve(s_field, mass, ctfield, option_path = path) |
|
143 |
end if
|
|
144 |
||
|
1092
by maddison
Add some caching to the finite_element_divergence diagnostic |
145 |
call deallocate(ct_m) |
|
773
by maddison
Adding a finite_element_divergence diagnostic algorithm. Bugfix in div algorithm and removing some invalid options from the shallow water schema |
146 |
call deallocate(ctfield) |
|
1093
by maddison
Applying usual tricks to make the divergence matrix caching work with mesh movement |
147 |
|
148 |
contains
|
|
149 |
||
150 |
function use_divergence_matrix_cache(state) result(use_cache) |
|
151 |
type(state_type), intent(in) :: state |
|
152 |
||
153 |
logical :: use_cache |
|
154 |
||
155 |
integer, save :: last_mesh_movement = -1 |
|
156 |
||
157 |
if(has_block_csr_matrix(state, trim(s_field%name) // "DivergenceMatrix")) then |
|
|
1132
by maddison
Adding some more diagnostics: vector sum, universal node numbering, and boundary condition option to finite element divergence transpose |
158 |
use_cache = (last_mesh_movement == eventcount(EVENT_MESH_MOVEMENT)) |
|
1093
by maddison
Applying usual tricks to make the divergence matrix caching work with mesh movement |
159 |
else
|
160 |
use_cache = .false. |
|
161 |
end if
|
|
|
1132
by maddison
Adding some more diagnostics: vector sum, universal node numbering, and boundary condition option to finite element divergence transpose |
162 |
|
163 |
if(use_cache) then |
|
164 |
ewrite(2, *) "Using cached divergence matrix" |
|
165 |
else
|
|
166 |
ewrite(2, *) "Not using cached divergence matrix" |
|
167 |
last_mesh_movement = eventcount(EVENT_MESH_MOVEMENT) |
|
168 |
end if
|
|
|
1093
by maddison
Applying usual tricks to make the divergence matrix caching work with mesh movement |
169 |
|
170 |
end function use_divergence_matrix_cache |
|
|
773
by maddison
Adding a finite_element_divergence diagnostic algorithm. Bugfix in div algorithm and removing some invalid options from the shallow water schema |
171 |
|
172 |
end subroutine calculate_finite_element_divergence |
|
|
825
by maddison
Adding finite_element_transpose diagnostic algorithm |
173 |
|
174 |
subroutine calculate_finite_element_divergence_transpose(state, v_field) |
|
175 |
type(state_type), intent(inout) :: state |
|
176 |
type(vector_field), intent(inout) :: v_field |
|
177 |
||
|
1132
by maddison
Adding some more diagnostics: vector sum, universal node numbering, and boundary condition option to finite element divergence transpose |
178 |
character(len = FIELD_NAME_LEN) :: bcfield_name |
|
825
by maddison
Adding finite_element_transpose diagnostic algorithm |
179 |
character(len = OPTION_PATH_LEN) :: path |
180 |
type(cmc_matrices) :: matrices |
|
181 |
type(scalar_field), pointer :: source_field |
|
|
1132
by maddison
Adding some more diagnostics: vector sum, universal node numbering, and boundary condition option to finite element divergence transpose |
182 |
type(vector_field), pointer :: bcfield |
|
825
by maddison
Adding finite_element_transpose diagnostic algorithm |
183 |
|
184 |
source_field => scalar_source_field(state, v_field) |
|
185 |
||
186 |
path = trim(complete_field_path(v_field%option_path)) // "/algorithm" |
|
|
1132
by maddison
Adding some more diagnostics: vector sum, universal node numbering, and boundary condition option to finite element divergence transpose |
187 |
if(have_option(trim(path) // "/bc_field")) then |
188 |
call get_option(trim(path) // "/bc_field/name", bcfield_name) |
|
189 |
bcfield => extract_vector_field(state, bcfield_name) |
|
190 |
call allocate(matrices, state, v_field, source_field, bcfield = bcfield, option_path = path, add_cmc = .false.) |
|
191 |
else
|
|
192 |
call allocate(matrices, state, v_field, source_field, option_path = path, add_cmc = .false.) |
|
193 |
end if
|
|
194 |
|
|
|
825
by maddison
Adding finite_element_transpose diagnostic algorithm |
195 |
call compute_conservative(matrices, v_field, source_field) |
196 |
call deallocate(matrices) |
|
197 |
||
198 |
end subroutine calculate_finite_element_divergence_transpose |
|
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
199 |
|
|
1488
by maddison
Perp diagnostic algorithm |
200 |
subroutine calculate_perp(state, v_field) |
201 |
type(state_type), intent(inout) :: state |
|
202 |
type(vector_field), intent(inout) :: v_field |
|
203 |
||
204 |
character(len = OPTION_PATH_LEN) :: path |
|
205 |
integer :: i |
|
206 |
type(csr_matrix), pointer :: mass |
|
207 |
type(scalar_field), pointer :: masslump |
|
208 |
type(scalar_field), pointer :: source_field |
|
209 |
type(vector_field) :: rhs |
|
210 |
type(vector_field), pointer :: positions |
|
211 |
||
212 |
write(1, *) "In calculate_perp" |
|
213 |
||
214 |
source_field => scalar_source_field(state, v_field) |
|
215 |
ewrite(2, *) "Calculating perp of field: " // trim(source_field%name) |
|
216 |
ewrite(2, *) "On mesh: " // trim(source_field%mesh%name) |
|
217 |
ewrite(2, *) "Diagnostic field: " // trim(v_field%name) |
|
218 |
ewrite(2, *) "On mesh: " // trim(v_field%mesh%name) |
|
219 |
||
220 |
if(v_field%dim /= 2) then |
|
|
1795
by hhiester
FLAbort/Exit updates |
221 |
FLExit("Can only calculate perp in 2D") |
|
1488
by maddison
Perp diagnostic algorithm |
222 |
end if
|
223 |
|
|
|
1997
by skramer
Taking the derivative of a discontinuous field is not supported in most diagnostic algorithms. Checks for this however were missing and fluidity would happily provide you with the wrong answer. The following diagnostic algorithms do not work for discontinuous source fields and now give an appropriate error message: |
224 |
call check_source_mesh_derivative(source_field, "perp") |
225 |
||
|
1488
by maddison
Perp diagnostic algorithm |
226 |
positions => extract_vector_field(state, "Coordinate") |
227 |
path = trim(complete_field_path(v_field%option_path)) // "/algorithm" |
|
228 |
||
229 |
if(have_option(trim(path) // "/lump_mass")) then |
|
230 |
masslump => get_lumped_mass(state, v_field%mesh) |
|
231 |
call zero(v_field) |
|
232 |
do i = 1, ele_count(v_field) |
|
233 |
call assemble_perp_ele(i, positions, source_field, v_field) |
|
234 |
end do
|
|
235 |
do i = 1, v_field%dim |
|
|
2445
by skramer
Changing the storage of vector_field from vfield%val(dim)%ptr(node) to vfield%val(dim,node). |
236 |
v_field%val(i,:) = v_field%val(i,:) / masslump%val |
|
1488
by maddison
Perp diagnostic algorithm |
237 |
end do
|
238 |
else
|
|
239 |
select case(continuity(v_field)) |
|
240 |
case(0) |
|
241 |
if(.not. have_option(trim(path) // "/solver")) then |
|
|
1531
by maddison
Very minor diagnostic output / FLAbort corrections |
242 |
FLExit("For continuous perp, must supply solver options when not lumping mass") |
|
1488
by maddison
Perp diagnostic algorithm |
243 |
end if
|
244 |
mass => get_mass_matrix(state, v_field%mesh) |
|
245 |
call allocate(rhs, v_field%dim, v_field%mesh, "PerpRHS") |
|
246 |
call zero(rhs) |
|
247 |
do i = 1, ele_count(rhs) |
|
248 |
call assemble_perp_ele(i, positions, source_field, rhs) |
|
249 |
end do
|
|
250 |
call petsc_solve(v_field, mass, rhs, option_path = path) |
|
251 |
call deallocate(rhs) |
|
252 |
case(-1) |
|
253 |
do i = 1, ele_count(v_field) |
|
254 |
call solve_perp_ele(i, positions, source_field, v_field) |
|
255 |
end do
|
|
256 |
case default |
|
257 |
ewrite(-1, *) "For mesh continuity: ", continuity(v_field) |
|
258 |
FLAbort("Unrecognised mesh continuity") |
|
259 |
end select
|
|
260 |
end if
|
|
261 |
|
|
262 |
contains
|
|
263 |
|
|
264 |
subroutine assemble_perp_ele(ele, positions, source, rhs) |
|
265 |
integer, intent(in) :: ele |
|
266 |
type(vector_field), intent(in) :: positions |
|
267 |
type(scalar_field), intent(in) :: source |
|
268 |
type(vector_field), intent(inout) :: rhs |
|
269 |
||
270 |
real, dimension(ele_ngi(positions, ele)) :: detwei |
|
271 |
real, dimension(ele_ngi(source, ele), positions%dim) :: grad_source_gi |
|
272 |
real, dimension(ele_loc(source, ele), ele_ngi(source, ele), positions%dim) :: dshape |
|
273 |
||
274 |
call transform_to_physical(positions, ele, ele_shape(source, ele), & |
|
275 |
& dshape = dshape, detwei = detwei) |
|
276 |
||
277 |
grad_source_gi = ele_grad_at_quad(source, ele, dshape) |
|
278 |
||
279 |
call addto(rhs, 1, ele_nodes(rhs, ele), & |
|
280 |
& shape_rhs(ele_shape(rhs, ele), grad_source_gi(:, 2) * detwei)) |
|
281 |
call addto(rhs, 2, ele_nodes(rhs, ele), & |
|
282 |
& shape_rhs(ele_shape(rhs, ele), -grad_source_gi(:, 1) * detwei)) |
|
283 |
||
284 |
end subroutine assemble_perp_ele |
|
285 |
||
286 |
subroutine solve_perp_ele(ele, positions, source, perp) |
|
287 |
integer, intent(in) :: ele |
|
288 |
type(vector_field), intent(in) :: positions |
|
289 |
type(scalar_field), intent(in) :: source |
|
290 |
type(vector_field), intent(inout) :: perp |
|
291 |
||
292 |
real, dimension(ele_ngi(positions, ele)) :: detwei |
|
293 |
real, dimension(ele_ngi(source, ele), positions%dim) :: grad_source_gi |
|
294 |
real, dimension(ele_loc(perp, ele), perp%dim) :: little_rhs |
|
295 |
real, dimension(ele_loc(perp, ele), ele_loc(perp, ele)) :: little_mass |
|
296 |
real, dimension(ele_loc(source, ele), ele_ngi(source, ele), positions%dim) :: dshape |
|
297 |
type(element_type), pointer :: shape |
|
298 |
|
|
299 |
call transform_to_physical(positions, ele, ele_shape(source, ele), & |
|
300 |
& dshape = dshape, detwei = detwei) |
|
301 |
||
302 |
grad_source_gi = ele_grad_at_quad(source, ele, dshape) |
|
303 |
||
304 |
shape => ele_shape(perp, ele) |
|
305 |
||
306 |
little_mass = shape_shape(shape, shape, detwei) |
|
307 |
little_rhs(:, 1) = shape_rhs(shape, grad_source_gi(:, 2) * detwei) |
|
308 |
little_rhs(:, 2) = shape_rhs(shape, -grad_source_gi(:, 1) * detwei) |
|
309 |
||
310 |
call solve(little_mass, little_rhs) |
|
311 |
||
312 |
call set(perp, ele_nodes(perp, ele), transpose(little_rhs)) |
|
313 |
||
314 |
end subroutine solve_perp_ele |
|
315 |
||
316 |
end subroutine calculate_perp |
|
317 |
||
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
318 |
subroutine calculate_curl_2d(state, s_field) |
319 |
type(state_type), intent(inout) :: state |
|
320 |
type(scalar_field), intent(inout) :: s_field |
|
321 |
||
322 |
character(len = OPTION_PATH_LEN) :: path |
|
323 |
integer :: i |
|
324 |
type(csr_matrix), pointer :: mass |
|
325 |
type(scalar_field), pointer :: masslump |
|
326 |
type(scalar_field) :: rhs |
|
327 |
type(vector_field), pointer :: positions, source_field |
|
328 |
||
329 |
ewrite(1, *) "In calculate_curl_2d" |
|
330 |
||
331 |
source_field => vector_source_field(state, s_field) |
|
332 |
ewrite(2, *) "Calculating curl of field: " // trim(source_field%name) |
|
333 |
ewrite(2, *) "On mesh: " // trim(source_field%mesh%name) |
|
334 |
ewrite(2, *) "Diagnostic field: " // trim(s_field%name) |
|
335 |
ewrite(2, *) "On mesh: " // trim(s_field%mesh%name) |
|
336 |
||
337 |
if(source_field%dim /= 2) then |
|
|
1795
by hhiester
FLAbort/Exit updates |
338 |
FLExit("Can only calculate 2D curl in 2D") |
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
339 |
end if
|
340 |
|
|
|
1997
by skramer
Taking the derivative of a discontinuous field is not supported in most diagnostic algorithms. Checks for this however were missing and fluidity would happily provide you with the wrong answer. The following diagnostic algorithms do not work for discontinuous source fields and now give an appropriate error message: |
341 |
call check_source_mesh_derivative(source_field, "curl_2d") |
342 |
||
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
343 |
positions => extract_vector_field(state, "Coordinate") |
344 |
path = trim(complete_field_path(s_field%option_path)) // "/algorithm" |
|
345 |
if(have_option(trim(path) // "/lump_mass")) then |
|
|
1071
by maddison
Removing sub-mesh lumping options for vorticity diagnostics |
346 |
masslump => get_lumped_mass(state, s_field%mesh) |
|
999
by maddison
Adding a .steady_state diagnostics file, for steady state diagnostics. Also marking some diagnostic fields as "legacy". |
347 |
call zero(s_field) |
|
1065
by maddison
Bugfix for lumped mass curl diagnostics |
348 |
do i = 1, ele_count(s_field) |
|
999
by maddison
Adding a .steady_state diagnostics file, for steady state diagnostics. Also marking some diagnostic fields as "legacy". |
349 |
call assemble_curl_ele(i, positions, source_field, s_field) |
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
350 |
end do
|
|
999
by maddison
Adding a .steady_state diagnostics file, for steady state diagnostics. Also marking some diagnostic fields as "legacy". |
351 |
s_field%val = s_field%val / masslump%val |
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
352 |
else
|
353 |
select case(continuity(s_field)) |
|
354 |
case(0) |
|
355 |
if(.not. have_option(trim(path) // "/solver")) then |
|
356 |
FLExit("For continuous curl, must supply solver options when not lumping mass") |
|
357 |
end if
|
|
358 |
mass => get_mass_matrix(state, s_field%mesh) |
|
359 |
call allocate(rhs, s_field%mesh, "CurlRHS") |
|
360 |
call zero(rhs) |
|
361 |
do i = 1, ele_count(rhs) |
|
362 |
call assemble_curl_ele(i, positions, source_field, rhs) |
|
363 |
end do
|
|
364 |
call petsc_solve(s_field, mass, rhs, option_path = path) |
|
365 |
call deallocate(rhs) |
|
366 |
case(-1) |
|
367 |
do i = 1, ele_count(s_field) |
|
368 |
call solve_curl_ele(i, positions, source_field, s_field) |
|
369 |
end do
|
|
370 |
case default |
|
371 |
ewrite(-1, *) "For mesh continuity: ", continuity(s_field) |
|
372 |
FLAbort("Unrecognised mesh continuity") |
|
373 |
end select
|
|
374 |
end if
|
|
375 |
|
|
376 |
ewrite_minmax(s_field%val) |
|
377 |
||
378 |
ewrite(1, *) "Exiting calculate_curl_2d" |
|
379 |
||
380 |
contains
|
|
381 |
|
|
382 |
subroutine assemble_curl_ele(ele, positions, source, rhs) |
|
383 |
integer, intent(in) :: ele |
|
384 |
type(vector_field), intent(in) :: positions |
|
385 |
type(vector_field), intent(in) :: source |
|
386 |
type(scalar_field), intent(inout) :: rhs |
|
387 |
||
388 |
real, dimension(ele_ngi(positions, ele)) :: detwei |
|
389 |
real, dimension(ele_loc(source, ele), ele_ngi(source, ele), positions%dim) :: dshape |
|
390 |
||
391 |
call transform_to_physical(positions, ele, ele_shape(source, ele), & |
|
392 |
& dshape = dshape, detwei = detwei) |
|
393 |
||
394 |
call addto(rhs, ele_nodes(rhs, ele), & |
|
395 |
& shape_rhs(ele_shape(rhs, ele), & |
|
396 |
& ele_2d_curl_at_quad(source, ele, dshape) * detwei)) |
|
397 |
||
398 |
end subroutine assemble_curl_ele |
|
399 |
||
400 |
subroutine solve_curl_ele(ele, positions, source, curl) |
|
401 |
integer, intent(in) :: ele |
|
402 |
type(vector_field), intent(in) :: positions |
|
403 |
type(vector_field), intent(in) :: source |
|
404 |
type(scalar_field), intent(inout) :: curl |
|
405 |
||
406 |
real, dimension(ele_ngi(positions, ele)) :: detwei |
|
407 |
real, dimension(ele_loc(curl, ele)) :: little_rhs |
|
408 |
real, dimension(ele_loc(curl, ele), ele_loc(curl, ele)) :: little_mass |
|
409 |
real, dimension(ele_loc(source, ele), ele_ngi(source, ele), positions%dim) :: dshape |
|
410 |
type(element_type), pointer :: shape |
|
411 |
|
|
412 |
call transform_to_physical(positions, ele, ele_shape(source, ele), & |
|
413 |
& dshape = dshape, detwei = detwei) |
|
414 |
||
415 |
shape => ele_shape(curl, ele) |
|
416 |
||
417 |
little_mass = shape_shape(shape, shape, detwei) |
|
418 |
little_rhs = shape_rhs(shape, ele_2d_curl_at_quad(source, ele, dshape) * detwei) |
|
419 |
||
420 |
call solve(little_mass, little_rhs) |
|
421 |
||
422 |
call set(curl, ele_nodes(curl, ele), little_rhs) |
|
423 |
||
424 |
end subroutine solve_curl_ele |
|
425 |
||
426 |
end subroutine calculate_curl_2d |
|
427 |
||
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
428 |
subroutine calculate_curl(state, v_field) |
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
429 |
type(state_type), intent(inout) :: state |
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
430 |
type(vector_field), intent(inout) :: v_field |
431 |
||
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
432 |
character(len = OPTION_PATH_LEN) :: path |
433 |
integer :: i |
|
434 |
type(csr_matrix), pointer :: mass |
|
435 |
type(scalar_field), pointer :: masslump |
|
436 |
type(vector_field) :: rhs |
|
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
437 |
type(vector_field), pointer :: positions, source_field |
438 |
||
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
439 |
ewrite(1, *) "In calculate_curl" |
440 |
||
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
441 |
source_field => vector_source_field(state, v_field) |
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
442 |
ewrite(2, *) "Calculating curl of field: " // trim(source_field%name) |
443 |
ewrite(2, *) "On mesh: " // trim(source_field%mesh%name) |
|
444 |
ewrite(2, *) "Diagnostic field: " // trim(v_field%name) |
|
445 |
ewrite(2, *) "On mesh: " // trim(v_field%mesh%name) |
|
446 |
||
447 |
if(source_field%dim /= 3) then |
|
|
1795
by hhiester
FLAbort/Exit updates |
448 |
FLExit("Can only calculate curl in 3D") |
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
449 |
end if
|
450 |
|
|
|
1997
by skramer
Taking the derivative of a discontinuous field is not supported in most diagnostic algorithms. Checks for this however were missing and fluidity would happily provide you with the wrong answer. The following diagnostic algorithms do not work for discontinuous source fields and now give an appropriate error message: |
451 |
call check_source_mesh_derivative(source_field, "curl") |
452 |
||
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
453 |
positions => extract_vector_field(state, "Coordinate") |
454 |
path = trim(complete_field_path(v_field%option_path)) // "/algorithm" |
|
455 |
if(have_option(trim(path) // "/lump_mass")) then |
|
|
1071
by maddison
Removing sub-mesh lumping options for vorticity diagnostics |
456 |
masslump => get_lumped_mass(state, v_field%mesh) |
|
999
by maddison
Adding a .steady_state diagnostics file, for steady state diagnostics. Also marking some diagnostic fields as "legacy". |
457 |
call zero(v_field) |
|
1065
by maddison
Bugfix for lumped mass curl diagnostics |
458 |
do i = 1, ele_count(v_field) |
|
999
by maddison
Adding a .steady_state diagnostics file, for steady state diagnostics. Also marking some diagnostic fields as "legacy". |
459 |
call assemble_curl_ele(i, positions, source_field, v_field) |
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
460 |
end do
|
461 |
do i = 1, v_field%dim |
|
|
2445
by skramer
Changing the storage of vector_field from vfield%val(dim)%ptr(node) to vfield%val(dim,node). |
462 |
v_field%val(i,:) = v_field%val(i,:) / masslump%val |
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
463 |
end do
|
464 |
else
|
|
465 |
select case(continuity(v_field)) |
|
466 |
case(0) |
|
467 |
if(.not. have_option(trim(path) // "/solver")) then |
|
468 |
FLExit("For continuous curl, must supply solver options when not lumping mass") |
|
469 |
end if
|
|
470 |
mass => get_mass_matrix(state, v_field%mesh) |
|
471 |
call allocate(rhs, v_field%dim, v_field%mesh, "CurlRHS") |
|
472 |
call zero(rhs) |
|
473 |
do i = 1, ele_count(rhs) |
|
474 |
call assemble_curl_ele(i, positions, source_field, rhs) |
|
475 |
end do
|
|
476 |
call petsc_solve(v_field, mass, rhs, option_path = path) |
|
477 |
call deallocate(rhs) |
|
478 |
case(-1) |
|
479 |
do i = 1, ele_count(v_field) |
|
480 |
call solve_curl_ele(i, positions, source_field, v_field) |
|
481 |
end do
|
|
482 |
case default |
|
483 |
ewrite(-1, *) "For mesh continuity: ", continuity(v_field) |
|
484 |
FLAbort("Unrecognised mesh continuity") |
|
485 |
end select
|
|
486 |
end if
|
|
487 |
|
|
488 |
do i = 1, v_field%dim |
|
|
2445
by skramer
Changing the storage of vector_field from vfield%val(dim)%ptr(node) to vfield%val(dim,node). |
489 |
ewrite_minmax(v_field%val(i,:)) |
|
993
by maddison
Adding some more options to the Vorticity and Vorticity2D diagnostics |
490 |
end do
|
491 |
|
|
492 |
ewrite(1, *) "Exiting calculate_curl" |
|
493 |
||
494 |
contains
|
|
495 |
|
|
496 |
subroutine assemble_curl_ele(ele, positions, source, rhs) |
|
497 |
integer, intent(in) :: ele |
|
498 |
type(vector_field), intent(in) :: positions |
|
499 |
type(vector_field), intent(in) :: source |
|
500 |
type(vector_field), intent(inout) :: rhs |
|
501 |
||
502 |
real, dimension(ele_ngi(positions, ele)) :: detwei |
|
503 |
real, dimension(ele_loc(source, ele), ele_ngi(source, ele), positions%dim) :: dshape |
|
504 |
||
505 |
call transform_to_physical(positions, ele, ele_shape(source, ele), & |
|
506 |
& dshape = dshape, detwei = detwei) |
|
507 |
||
508 |
call addto(rhs, ele_nodes(rhs, ele), & |
|
509 |
& shape_vector_rhs(ele_shape(rhs, ele), & |
|
510 |
& ele_curl_at_quad(source, ele, dshape), detwei)) |
|
511 |
||
512 |
end subroutine assemble_curl_ele |
|
513 |
||
514 |
subroutine solve_curl_ele(ele, positions, source, curl) |
|
515 |
integer, intent(in) :: ele |
|
516 |
type(vector_field), intent(in) :: positions |
|
517 |
type(vector_field), intent(in) :: source |
|
518 |
type(vector_field), intent(inout) :: curl |
|
519 |
||
520 |
real, dimension(ele_ngi(positions, ele)) :: detwei |
|
521 |
real, dimension(ele_loc(curl, ele), curl%dim) :: little_rhs |
|
522 |
real, dimension(ele_loc(curl, ele), ele_loc(curl, ele)) :: little_mass |
|
523 |
real, dimension(ele_loc(source, ele), ele_ngi(source, ele), positions%dim) :: dshape |
|
524 |
type(element_type), pointer :: shape |
|
525 |
|
|
526 |
call transform_to_physical(positions, ele, ele_shape(source, ele), & |
|
527 |
& dshape = dshape, detwei = detwei) |
|
528 |
||
529 |
shape => ele_shape(curl, ele) |
|
530 |
||
531 |
little_mass = shape_shape(shape, shape, detwei) |
|
532 |
little_rhs = transpose(shape_vector_rhs(shape, & |
|
533 |
& ele_curl_at_quad(source, ele, dshape), detwei)) |
|
534 |
||
535 |
call solve(little_mass, little_rhs) |
|
536 |
||
537 |
call set(curl, ele_nodes(curl, ele), transpose(little_rhs)) |
|
538 |
||
539 |
end subroutine solve_curl_ele |
|
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
540 |
|
541 |
end subroutine calculate_curl |
|
542 |
||
543 |
subroutine calculate_scalar_advection(state, s_field) |
|
544 |
type(state_type), intent(in) :: state |
|
545 |
type(scalar_field), intent(inout) :: s_field |
|
546 |
||
547 |
integer :: stat |
|
548 |
type(scalar_field), pointer :: source_field |
|
549 |
type(vector_field), pointer :: positions, velocity |
|
550 |
||
551 |
source_field => scalar_source_field(state, s_field) |
|
552 |
||
553 |
velocity => extract_vector_field(state, "Velocity", stat = stat) |
|
554 |
if(stat /= 0) then |
|
555 |
ewrite(0, *) "For field " // trim(s_field%name) |
|
556 |
ewrite(0, *) "Warning: Calculating advection with no velocity" |
|
557 |
call zero(s_field) |
|
558 |
return
|
|
559 |
end if
|
|
560 |
|
|
|
1997
by skramer
Taking the derivative of a discontinuous field is not supported in most diagnostic algorithms. Checks for this however were missing and fluidity would happily provide you with the wrong answer. The following diagnostic algorithms do not work for discontinuous source fields and now give an appropriate error message: |
561 |
call check_source_mesh_derivative(source_field, "scalar_advection") |
562 |
||
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
563 |
positions => extract_vector_field(state, "Coordinate") |
564 |
||
565 |
call u_dot_nabla(velocity, source_field, positions, s_field) |
|
566 |
||
567 |
end subroutine calculate_scalar_advection |
|
568 |
||
569 |
subroutine calculate_vector_advection(state, v_field) |
|
570 |
type(state_type), intent(in) :: state |
|
571 |
type(vector_field), intent(inout) :: v_field |
|
572 |
||
573 |
integer :: stat |
|
574 |
type(vector_field), pointer :: positions, source_field, velocity |
|
575 |
||
576 |
source_field => vector_source_field(state, v_field) |
|
577 |
||
578 |
velocity => extract_vector_field(state, "Velocity", stat = stat) |
|
579 |
if(stat /= 0) then |
|
580 |
ewrite(0, *) "For field " // trim(v_field%name) |
|
581 |
ewrite(0, *) "Warning: Calculating advection with no velocity" |
|
582 |
call zero(v_field) |
|
583 |
return
|
|
584 |
end if
|
|
585 |
|
|
|
1997
by skramer
Taking the derivative of a discontinuous field is not supported in most diagnostic algorithms. Checks for this however were missing and fluidity would happily provide you with the wrong answer. The following diagnostic algorithms do not work for discontinuous source fields and now give an appropriate error message: |
586 |
call check_source_mesh_derivative(source_field, "vector_advection") |
587 |
||
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
588 |
positions => extract_vector_field(state, "Coordinate") |
589 |
||
590 |
call u_dot_nabla(velocity, source_field, positions, v_field) |
|
591 |
||
592 |
end subroutine calculate_vector_advection |
|
593 |
||
594 |
subroutine calculate_vector_laplacian(state, v_field) |
|
595 |
type(state_type), intent(inout) :: state |
|
596 |
type(vector_field), intent(inout) :: v_field |
|
597 |
||
|
1555
by maddison
Add on the boundary term for scalar_laplacian and vector_laplacian algorithms |
598 |
integer :: i, j |
599 |
type(scalar_field) :: source_field_comp, v_field_comp |
|
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
600 |
type(scalar_field), pointer :: masslump |
601 |
type(vector_field), pointer :: positions, source_field |
|
602 |
||
603 |
ewrite(1, *) "In calculate_vector_laplacian" |
|
604 |
ewrite(2, *) "Computing laplacian for field " // trim(v_field%name) |
|
605 |
||
606 |
source_field => vector_source_field(state, v_field) |
|
607 |
assert(ele_count(source_field) == ele_count(v_field)) |
|
608 |
assert(source_field%dim == v_field%dim) |
|
609 |
do i = 1, source_field%dim |
|
|
2445
by skramer
Changing the storage of vector_field from vfield%val(dim)%ptr(node) to vfield%val(dim,node). |
610 |
ewrite_minmax(source_field%val(i,:)) |
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
611 |
end do
|
|
1997
by skramer
Taking the derivative of a discontinuous field is not supported in most diagnostic algorithms. Checks for this however were missing and fluidity would happily provide you with the wrong answer. The following diagnostic algorithms do not work for discontinuous source fields and now give an appropriate error message: |
612 |
|
613 |
call check_source_mesh_derivative(source_field, "vector_laplacian") |
|
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
614 |
|
615 |
positions => extract_vector_field(state, "Coordinate") |
|
616 |
assert(ele_count(positions) == ele_count(v_field)) |
|
617 |
||
618 |
call zero(v_field) |
|
619 |
do i = 1, ele_count(v_field) |
|
620 |
call assemble_vector_laplacian_ele(i, v_field, positions, source_field) |
|
621 |
end do
|
|
622 |
do i = 1, v_field%dim |
|
|
1555
by maddison
Add on the boundary term for scalar_laplacian and vector_laplacian algorithms |
623 |
v_field_comp = extract_scalar_field(v_field, i) |
624 |
source_field_comp = extract_scalar_field(source_field, i) |
|
625 |
do j = 1, surface_element_count(v_field) |
|
626 |
! This could be made more efficent by assembling all components at the
|
|
627 |
! same time
|
|
628 |
call assemble_scalar_laplacian_face(j, v_field_comp, positions, source_field_comp) |
|
629 |
end do
|
|
630 |
end do
|
|
631 |
do i = 1, v_field%dim |
|
|
2445
by skramer
Changing the storage of vector_field from vfield%val(dim)%ptr(node) to vfield%val(dim,node). |
632 |
ewrite_minmax(v_field%val(i,:)) |
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
633 |
end do
|
634 |
|
|
635 |
masslump => get_lumped_mass(state, v_field%mesh) |
|
636 |
||
637 |
do i = 1, v_field%dim |
|
|
2445
by skramer
Changing the storage of vector_field from vfield%val(dim)%ptr(node) to vfield%val(dim,node). |
638 |
v_field%val(i,:) = v_field%val(i,:) / masslump%val |
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
639 |
end do
|
640 |
do i = 1, v_field%dim |
|
|
2445
by skramer
Changing the storage of vector_field from vfield%val(dim)%ptr(node) to vfield%val(dim,node). |
641 |
ewrite_minmax(v_field%val(i,:)) |
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
642 |
end do
|
643 |
|
|
644 |
ewrite(1, *) "Exiting calculate_vector_laplacian" |
|
645 |
||
646 |
end subroutine calculate_vector_laplacian |
|
647 |
||
648 |
subroutine assemble_vector_laplacian_ele(ele, v_field, positions, source_field) |
|
649 |
integer, intent(in) :: ele |
|
650 |
type(vector_field), intent(inout) :: v_field |
|
651 |
type(vector_field), intent(in) :: positions |
|
652 |
type(vector_field), intent(in) :: source_field |
|
653 |
||
654 |
integer :: i, j |
|
655 |
integer, dimension(:), pointer :: element_nodes |
|
656 |
real, dimension(ele_ngi(v_field, ele)) :: detwei |
|
657 |
real, dimension(source_field%dim, ele_ngi(source_field, ele)) :: grad_gi |
|
658 |
real, dimension(ele_loc(source_field, ele), ele_ngi(source_field, ele), source_field%dim) :: dn_t |
|
659 |
||
660 |
assert(ele_ngi(positions, ele) == ele_ngi(v_field, ele)) |
|
661 |
assert(ele_ngi(source_field, ele) == ele_ngi(v_field, ele)) |
|
662 |
||
663 |
call transform_to_physical(positions, ele, ele_shape(source_field, ele), & |
|
664 |
& dshape = dn_t, detwei = detwei) |
|
665 |
||
666 |
element_nodes => ele_nodes(v_field, ele) |
|
667 |
||
668 |
do i = 1, source_field%dim |
|
669 |
do j = 1, source_field%dim |
|
670 |
grad_gi(j, :) = matmul(ele_val(source_field, i, ele), dn_t(:, :, j)) |
|
671 |
end do
|
|
672 |
|
|
673 |
call addto(v_field, i, element_nodes, -dshape_dot_vector_rhs(dn_t, grad_gi, detwei)) |
|
674 |
end do
|
|
675 |
|
|
676 |
end subroutine assemble_vector_laplacian_ele |
|
677 |
||
|
1555
by maddison
Add on the boundary term for scalar_laplacian and vector_laplacian algorithms |
678 |
subroutine assemble_scalar_laplacian_face(face, s_field, positions, source_field) |
679 |
integer, intent(in) :: face |
|
680 |
type(scalar_field), intent(inout) :: s_field |
|
681 |
type(vector_field), intent(in) :: positions |
|
682 |
type(scalar_field), intent(in) :: source_field |
|
683 |
||
684 |
integer :: ele, lface, i, j |
|
685 |
real, dimension(face_ngi(s_field, face)) :: detwei, grad_sgi_n |
|
686 |
real, dimension(positions%dim, face_ngi(s_field, face)) :: grad_sgi, normal |
|
687 |
real, dimension(positions%dim, positions%dim, ele_ngi(s_field, face_ele(s_field, face))) :: invj |
|
688 |
real, dimension(positions%dim, positions%dim, face_ngi(s_field, face)) :: invj_face |
|
689 |
real, dimension(ele_loc(s_field, face_ele(source_field, face)), face_ngi(s_field, face), positions%dim) :: dshape_face |
|
690 |
real, dimension(ele_loc(s_field, face_ele(source_field, face))) :: s_ele_val |
|
691 |
type(element_type) :: augmented_shape |
|
692 |
type(element_type), pointer :: fshape, source_shape, positions_shape |
|
693 |
||
694 |
ele = face_ele(s_field, face) |
|
695 |
lface = local_face_number(s_field, face) |
|
696 |
||
697 |
positions_shape => ele_shape(positions, ele) |
|
698 |
fshape => face_shape(s_field, face) |
|
699 |
||
700 |
source_shape => ele_shape(source_field, ele) |
|
701 |
if(associated(source_shape%dn_s)) then |
|
702 |
augmented_shape = source_shape |
|
703 |
call incref(augmented_shape) |
|
704 |
else
|
|
705 |
augmented_shape = make_element_shape(positions_shape%loc, source_shape%dim, & |
|
706 |
& source_shape%degree, source_shape%quadrature, & |
|
707 |
& quad_s = fshape%quadrature) |
|
708 |
end if
|
|
709 |
|
|
710 |
call transform_facet_to_physical(positions, face, & |
|
711 |
& detwei_f = detwei, normal = normal) |
|
712 |
call compute_inverse_jacobian( & |
|
713 |
& ele_val(positions, ele), positions_shape, & |
|
714 |
& invj = invj) |
|
715 |
assert(positions_shape%degree == 1) |
|
716 |
assert(ele_numbering_family(positions_shape) == FAMILY_SIMPLEX) |
|
717 |
invj_face = spread(invj(:, :, 1), 3, size(invj_face, 3)) |
|
718 |
||
719 |
dshape_face = eval_volume_dshape_at_face_quad(augmented_shape, lface, invj_face) |
|
720 |
call deallocate(augmented_shape) |
|
721 |
||
722 |
s_ele_val = ele_val(source_field, ele) |
|
723 |
forall(i = 1:positions%dim, j = 1:face_ngi(s_field, face)) |
|
724 |
grad_sgi(i, j) = dot_product(s_ele_val, dshape_face(:, j, i)) |
|
725 |
end forall
|
|
726 |
|
|
727 |
do i = 1, face_ngi(s_field, face) |
|
728 |
grad_sgi_n = dot_product(grad_sgi(:, i), normal(:, i)) |
|
729 |
end do
|
|
730 |
|
|
731 |
call addto(s_field, face_global_nodes(s_field, face), shape_rhs(fshape, detwei * grad_sgi_n)) |
|
732 |
||
733 |
end subroutine assemble_scalar_laplacian_face |
|
734 |
||
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
735 |
subroutine calculate_scalar_laplacian(state, s_field) |
736 |
type(state_type), intent(inout) :: state |
|
737 |
type(scalar_field), intent(inout) :: s_field |
|
738 |
||
739 |
integer :: i |
|
740 |
type(scalar_field), pointer :: masslump, source_field |
|
741 |
type(vector_field), pointer :: positions |
|
742 |
||
743 |
ewrite(1, *) "In calculate_scalar_laplacian" |
|
744 |
ewrite(2, *) "Computing laplacian for field " // trim(s_field%name) |
|
745 |
||
746 |
source_field => scalar_source_field(state, s_field) |
|
747 |
assert(ele_count(source_field) == ele_count(s_field)) |
|
748 |
assert(mesh_dim(source_field) == mesh_dim(s_field)) |
|
749 |
ewrite_minmax(source_field%val) |
|
750 |
||
|
1997
by skramer
Taking the derivative of a discontinuous field is not supported in most diagnostic algorithms. Checks for this however were missing and fluidity would happily provide you with the wrong answer. The following diagnostic algorithms do not work for discontinuous source fields and now give an appropriate error message: |
751 |
call check_source_mesh_derivative(source_field, "scalar_laplacian") |
752 |
||
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
753 |
positions => extract_vector_field(state, "Coordinate") |
754 |
assert(ele_count(positions) == ele_count(s_field)) |
|
755 |
||
756 |
call zero(s_field) |
|
757 |
do i = 1, ele_count(s_field) |
|
758 |
call assemble_scalar_laplacian_ele(i, s_field, positions, source_field) |
|
759 |
end do
|
|
|
1555
by maddison
Add on the boundary term for scalar_laplacian and vector_laplacian algorithms |
760 |
do i = 1, surface_element_count(s_field) |
761 |
call assemble_scalar_laplacian_face(i, s_field, positions, source_field) |
|
762 |
end do
|
|
|
1
by hhiester
deleting initialisation as none of tools are used any longer. |
763 |
ewrite_minmax(s_field%val) |
764 |
||
765 |
masslump => get_lumped_mass(state, s_field%mesh) |
|
766 |
||
767 |
s_field%val = s_field%val / masslump%val |
|
768 |
ewrite_minmax(s_field%val) |
|
769 |
||
770 |
ewrite(1, *) "Exiting calculate_scalar_laplacian" |
|
771 |
||
772 |
end subroutine calculate_scalar_laplacian |
|
773 |
||
774 |
subroutine assemble_scalar_laplacian_ele(ele, s_field, positions, source_field) |
|
775 |
integer, intent(in) :: ele |
|
776 |
type(scalar_field), intent(inout) :: s_field |
|
777 |
type(vector_field), intent(in) :: positions |
|
778 |
type(scalar_field), intent(in) :: source_field |
|
779 |
||
780 |
integer, dimension(:), pointer :: element_nodes |
|
781 |
real, dimension(ele_ngi(s_field, ele)) :: detwei |
|
782 |
real, dimension(ele_loc(source_field, ele), ele_ngi(source_field, ele), mesh_dim(source_field)) :: dn_t |
|
783 |
||
784 |
assert(ele_ngi(positions, ele) == ele_ngi(s_field, ele)) |
|
785 |
assert(ele_ngi(source_field, ele) == ele_ngi(s_field, ele)) |
|
786 |
||
787 |
call transform_to_physical(positions, ele, ele_shape(source_field, ele), & |
|
788 |
& dshape = dn_t, detwei = detwei) |
|
789 |
||
790 |
element_nodes => ele_nodes(s_field, ele) |
|
791 |
||
792 |
call addto(s_field, element_nodes, -dshape_dot_vector_rhs(dn_t, transpose(ele_grad_at_quad(source_field, ele, dn_t)), detwei)) |
|
793 |
||
794 |
end subroutine assemble_scalar_laplacian_ele |
|
795 |
||
796 |
end module differential_operator_diagnostics |