Merge lp:~btollit/fluidity/pressure_CV_weak_pressure_BC into lp:fluidity
- pressure_CV_weak_pressure_BC
- Merge into dev-trunk
Status: | Merged |
---|---|
Approved by: | Cian Wilson |
Approved revision: | 3873 |
Merged at revision: | 3882 |
Proposed branch: | lp:~btollit/fluidity/pressure_CV_weak_pressure_BC |
Merge into: | lp:fluidity |
Diff against target: |
7063 lines (+6824/-2) 40 files modified
assemble/Makefile.dependencies (+17/-1) assemble/Makefile.in (+1/-1) assemble/Momentum_Equation.F90 (+6/-0) assemble/Pressure_Dirichlet_BCS_CV.F90 (+280/-0) tests/darcy_p0p1_test_cty_cv_pressBCinlet/Makefile (+32/-0) tests/darcy_p0p1_test_cty_cv_pressBCinlet/cube.geo (+39/-0) tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet.xml (+118/-0) tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_1d.flml (+411/-0) tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_2d.flml (+427/-0) tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_3d.flml (+442/-0) tests/darcy_p0p1_test_cty_cv_pressBCinlet/square.geo (+22/-0) tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/Makefile (+23/-0) tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain.xml (+46/-0) tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d.flml (+372/-0) tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/square_rotated.geo (+38/-0) tests/darcy_p0p1_test_cty_cv_velBCinlet/Makefile (+32/-0) tests/darcy_p0p1_test_cty_cv_velBCinlet/cube.geo (+39/-0) tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet.xml (+94/-0) tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_1d.flml (+391/-0) tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_2d.flml (+407/-0) tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_3d.flml (+404/-0) tests/darcy_p0p1_test_cty_cv_velBCinlet/square.geo (+22/-0) tests/darcy_p0p1cv_pressBCinlet/Makefile (+32/-0) tests/darcy_p0p1cv_pressBCinlet/cube.geo (+39/-0) tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet.xml (+94/-0) tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_1d.flml (+384/-0) tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_2d.flml (+392/-0) tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_3d.flml (+398/-0) tests/darcy_p0p1cv_pressBCinlet/square.geo (+22/-0) tests/darcy_p0p1cv_pressBCinlet_1solidphase/Makefile (+32/-0) tests/darcy_p0p1cv_pressBCinlet_1solidphase/cube.geo (+39/-0) tests/darcy_p0p1cv_pressBCinlet_1solidphase/darcy_p0p1cv_pressBCinlet_1solidphase.xml (+94/-0) tests/darcy_p0p1cv_pressBCinlet_1solidphase/darcy_p0p1cv_pressBCinlet_1solidphase_1d.flml (+363/-0) tests/darcy_p0p1cv_pressBCinlet_1solidphase/darcy_p0p1cv_pressBCinlet_1solidphase_2d.flml (+379/-0) tests/darcy_p0p1cv_pressBCinlet_1solidphase/darcy_p0p1cv_pressBCinlet_1solidphase_3d.flml (+394/-0) tests/darcy_p0p1cv_pressBCinlet_1solidphase/square.geo (+22/-0) tests/darcy_p0p1cv_pressBCinlet_rotated_domain/Makefile (+23/-0) tests/darcy_p0p1cv_pressBCinlet_rotated_domain/darcy_p0p1cv_pressBCinlet_rotated_domain.xml (+46/-0) tests/darcy_p0p1cv_pressBCinlet_rotated_domain/darcy_p0p1cv_pressBCinlet_rotated_domain_2d.flml (+370/-0) tests/darcy_p0p1cv_pressBCinlet_rotated_domain/square_rotated.geo (+38/-0) |
To merge this branch: | bzr merge lp:~btollit/fluidity/pressure_CV_weak_pressure_BC |
Related bugs: |
Reviewer | Review Type | Date Requested | Status |
---|---|---|---|
Cian Wilson | Approve | ||
Review via email: mp+83558@code.launchpad.net |
Commit message
Description of the change
This merge will add a procedure to be able to have non zero pressure Dirichlet boundary conditions when using a CV discretisation for pressure. As for the CG pressure BCs the integrals are included for both weak and strong Dirichlet. For multiphase the phase volume fraction is included via a FE spatial interpolation, akin to their inclusion in interior domain terms.
Code Additions:
- assemble/
Code Changes:
- assemble/
- assemble/
- assemble/
Test Cases Removed:
- tests/mphase_
Test Cases Added:
- tests/darcy_
- tests/darcy_
- tests/darcy_
- tests/darcy_
- tests/darcy_
- tests/darcy_
Buildbot green (on Nov 25 18:48)
- 3872. By Brendan Tollit
-
Very small change to a comment.
- 3873. By Brendan Tollit
-
Merge trunk into branch
Brendan Tollit (btollit) wrote : | # |
Cian Wilson (cwilson) wrote : | # |
Looks good. Sorry for the delay in reviewing.
One comment/suggestion:
In Makefile.in:
Pressure_
Clearly, as it still compiles this doesn't make much difference!
- 3874. By Brendan Tollit
-
Change .F90 to .o spotted by Cian in review.
- 3875. By Brendan Tollit
-
Merge trunk into branch
Preview Diff
1 | === modified file 'assemble/Makefile.dependencies' |
2 | --- assemble/Makefile.dependencies 2011-10-27 23:11:25 +0000 |
3 | +++ assemble/Makefile.dependencies 2011-12-23 12:05:25 +0000 |
4 | @@ -579,7 +579,8 @@ |
5 | ../include/momentum_cg.mod ../include/momentum_dg.mod \ |
6 | ../include/momentum_diagnostic_fields.mod ../include/multiphase_module.mod \ |
7 | ../include/oceansurfaceforcing.mod ../include/parallel_tools.mod \ |
8 | - ../include/petsc_solve_state_module.mod ../include/profiler.mod \ |
9 | + ../include/petsc_solve_state_module.mod \ |
10 | + ../include/pressure_dirichlet_bcs_cv.mod ../include/profiler.mod \ |
11 | ../include/reduced_model_runtime.mod \ |
12 | ../include/rotated_boundary_conditions.mod ../include/slope_limiters_dg.mod \ |
13 | ../include/solvers.mod ../include/sparse_matrices_fields.mod \ |
14 | @@ -645,6 +646,21 @@ |
15 | ../include/fdebug.h ../include/fields.mod ../include/fldebug.mod \ |
16 | ../include/spud.mod ../include/state_module.mod |
17 | |
18 | +../include/pressure_dirichlet_bcs_cv.mod: Pressure_Dirichlet_BCS_CV.o |
19 | + @true |
20 | + |
21 | +Pressure_Dirichlet_BCS_CV.o ../include/pressure_dirichlet_bcs_cv.mod: \ |
22 | + Pressure_Dirichlet_BCS_CV.F90 ../include/boundary_conditions.mod \ |
23 | + ../include/cv_face_values.mod ../include/cv_faces.mod \ |
24 | + ../include/cv_fields.mod ../include/cv_options.mod \ |
25 | + ../include/cv_shape_functions.mod ../include/cv_upwind_values.mod \ |
26 | + ../include/cvtools.mod ../include/diagnostic_fields.mod ../include/fdebug.h \ |
27 | + ../include/fetools.mod ../include/field_derivatives.mod \ |
28 | + ../include/field_options.mod ../include/fields.mod ../include/futils.mod \ |
29 | + ../include/global_parameters.mod ../include/multiphase_module.mod \ |
30 | + ../include/quadrature.mod ../include/sparsity_patterns_meshes.mod \ |
31 | + ../include/spud.mod ../include/state_module.mod |
32 | + |
33 | ../include/pseudo_supermesh.mod: Pseudo_supermesh.o |
34 | @true |
35 | |
36 | |
37 | === modified file 'assemble/Makefile.in' |
38 | --- assemble/Makefile.in 2011-10-21 10:05:25 +0000 |
39 | +++ assemble/Makefile.in 2011-12-23 12:05:25 +0000 |
40 | @@ -76,7 +76,7 @@ |
41 | Linear_Shallow_Water.o \ |
42 | Zoltan_global_variables.o Zoltan_integration.o Zoltan_callbacks.o Zoltan_detectors.o \ |
43 | Diagnostic_Children.o Reduced_Model_Runtime.o Turbine.o Implicit_Solids.o \ |
44 | - Manifold_Projections.o Hybridized_Helmholtz.o Burgers_Assembly.o |
45 | + Manifold_Projections.o Hybridized_Helmholtz.o Burgers_Assembly.o Pressure_Dirichlet_BCS_CV.o |
46 | |
47 | .SUFFIXES: .F90 .c .cpp .o .a |
48 | |
49 | |
50 | === modified file 'assemble/Momentum_Equation.F90' |
51 | --- assemble/Momentum_Equation.F90 2011-11-23 17:00:55 +0000 |
52 | +++ assemble/Momentum_Equation.F90 2011-12-23 12:05:25 +0000 |
53 | @@ -74,6 +74,7 @@ |
54 | use slope_limiters_dg |
55 | use implicit_solids |
56 | use multiphase_module |
57 | + use pressure_dirichlet_bcs_cv |
58 | |
59 | implicit none |
60 | |
61 | @@ -589,6 +590,11 @@ |
62 | include_pressure_and_continuity_bcs=.not. cv_pressure) |
63 | end if |
64 | |
65 | + ! If CV pressure then add in any dirichlet pressure BC integrals to the mom_rhs. |
66 | + if (cv_pressure) then |
67 | + call add_pressure_dirichlet_bcs_cv(mom_rhs(istate), u, p, state(istate)) |
68 | + end if |
69 | + |
70 | ! Add in multiphase interactions (e.g. fluid-particle drag) if necessary |
71 | ! Note: this is done outside of construct_momentum_cg/dg to keep things |
72 | ! neater in Momentum_CG/DG.F90, since we would need to pass around multiple phases |
73 | |
74 | === added file 'assemble/Pressure_Dirichlet_BCS_CV.F90' |
75 | --- assemble/Pressure_Dirichlet_BCS_CV.F90 1970-01-01 00:00:00 +0000 |
76 | +++ assemble/Pressure_Dirichlet_BCS_CV.F90 2011-12-23 12:05:25 +0000 |
77 | @@ -0,0 +1,280 @@ |
78 | +! Copyright (C) 2006 Imperial College London and others. |
79 | +! |
80 | +! Please see the AUTHORS file in the main source directory for a full list |
81 | +! of copyright holders. |
82 | +! |
83 | +! Prof. C Pain |
84 | +! Applied Modelling and Computation Group |
85 | +! Department of Earth Science and Engineering |
86 | +! Imperial College London |
87 | +! |
88 | +! amcgsoftware@imperial.ac.uk |
89 | +! |
90 | +! This library is free software; you can redistribute it and/or |
91 | +! modify it under the terms of the GNU Lesser General Public |
92 | +! License as published by the Free Software Foundation, |
93 | +! version 2.1 of the License. |
94 | +! |
95 | +! This library is distributed in the hope that it will be useful, |
96 | +! but WITHOUT ANY WARRANTY; without even the implied warranty of |
97 | +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
98 | +! Lesser General Public License for more details. |
99 | +! |
100 | +! You should have received a copy of the GNU Lesser General Public |
101 | +! License along with this library; if not, write to the Free Software |
102 | +! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 |
103 | +! USA |
104 | +#include "fdebug.h" |
105 | + |
106 | +module pressure_dirichlet_bcs_cv |
107 | + |
108 | + use quadrature |
109 | + use fields |
110 | + use field_derivatives |
111 | + use state_module |
112 | + use futils |
113 | + use fetools |
114 | + use spud |
115 | + use cv_faces |
116 | + use cv_shape_functions |
117 | + use cvtools |
118 | + use cv_fields |
119 | + use cv_upwind_values |
120 | + use cv_face_values |
121 | + use cv_options |
122 | + use diagnostic_fields, only: calculate_diagnostic_variable |
123 | + use boundary_conditions |
124 | + use global_parameters, only: OPTION_PATH_LEN |
125 | + use field_options, only: get_coordinate_field |
126 | + use sparsity_patterns_meshes |
127 | + use multiphase_module |
128 | + |
129 | + implicit none |
130 | + |
131 | + private |
132 | + |
133 | + public :: add_pressure_dirichlet_bcs_cv |
134 | + |
135 | +contains |
136 | + |
137 | + ! -------------------------------------------------------------------------------------- |
138 | + |
139 | + subroutine add_pressure_dirichlet_bcs_cv(mom_rhs, u, p, state) |
140 | + |
141 | + !!< Add any CV pressure dirichlet BC integrals to the mom_rhs field if required. |
142 | + !!< If this is a multiphase simulation then the phase volume fraction spatial |
143 | + !!< interpolation uses finite elements. This isnt striclty correct if the |
144 | + !!< phase volume fractions were solved with a control volume discretisation |
145 | + !!< and also whether of not they themselves have a weak BC is ignored. If the |
146 | + !!< latter was to be considered then the upwind value should be used. |
147 | + |
148 | + ! Note the pressure BC values come from the end of time step Pressure field |
149 | + ! rather than using a theta average. The latter was tried but didnt work |
150 | + ! because OldPressure had no BC info. Taking the end of time step value |
151 | + ! for the BC is also what the pressure CG routine(s) do. |
152 | + |
153 | + type(vector_field), intent(inout) :: mom_rhs |
154 | + type(vector_field), intent(in) :: u |
155 | + type(scalar_field), intent(in) :: p |
156 | + type(state_type), intent(inout) :: state ! required for phase volume fraction |
157 | + |
158 | + ! local variables |
159 | + |
160 | + ! degree of quadrature over cv faces |
161 | + integer :: quaddegree |
162 | + |
163 | + real, dimension(:,:), allocatable :: x_ele, x_ele_bdy |
164 | + real, dimension(:,:), allocatable :: normal_bdy |
165 | + real, dimension(:), allocatable :: detwei_bdy |
166 | + integer, dimension(:), allocatable :: u_nodes_bdy |
167 | + |
168 | + integer :: ele, sele, iloc, jloc, face, gi, ggi, dim |
169 | + |
170 | + ! information about cv faces |
171 | + type(cv_faces_type) :: cvfaces |
172 | + ! shape functions for surface |
173 | + type(element_type) :: x_cvbdyshape |
174 | + type(element_type) :: u_cvbdyshape |
175 | + |
176 | + ! pointer to coordinates |
177 | + type(vector_field), pointer :: x |
178 | + |
179 | + ! pressure BC info |
180 | + integer, dimension(:), allocatable :: pressure_bc_type |
181 | + real, dimension(:), allocatable :: pressure_bc_val |
182 | + type(scalar_field) :: pressure_bc |
183 | + |
184 | + ! local ct matrix |
185 | + real, dimension(:,:,:), allocatable :: ct_mat_local_bdy |
186 | + |
187 | + ! Multiphase variables |
188 | + logical :: multiphase |
189 | + ! Volume fraction fields |
190 | + type(scalar_field), pointer :: vfrac |
191 | + type(scalar_field) :: nvfrac |
192 | + |
193 | + ! Variables used for FE interpolation of vfrac on surface |
194 | + ! Values of nvfrac at the faces |
195 | + real, dimension(:), allocatable :: nvfrac_gi_f |
196 | + ! CV shape functions for nvfrac (computed only if the Coordinate and |
197 | + ! PhaseVolumeFraction meshes are different, otherwise it is |
198 | + ! assigned x_cvbdyshape) |
199 | + type(element_type) :: nvfrac_cvbdyshape |
200 | + |
201 | + ewrite(1,*) 'In add_cv_pressure_weak_dirichlet_bcs' |
202 | + |
203 | + x => extract_vector_field(state, "Coordinate") |
204 | + |
205 | + allocate(x_ele(x%dim, x%mesh%shape%loc)) |
206 | + |
207 | + call get_option("/geometry/quadrature/controlvolume_surface_degree", quaddegree, default=1) |
208 | + |
209 | + cvfaces = find_cv_faces(vertices = ele_vertices(p, 1), & |
210 | + dimension = mesh_dim(p), & |
211 | + polydegree = p%mesh%shape%degree, & |
212 | + quaddegree = quaddegree) |
213 | + |
214 | + x_cvbdyshape = make_cvbdy_element_shape(cvfaces, x%mesh%faces%shape) |
215 | + u_cvbdyshape = make_cvbdy_element_shape(cvfaces, u%mesh%faces%shape) |
216 | + |
217 | + assert(surface_element_count(p) == surface_element_count(u)) |
218 | + |
219 | + ! get the pressure BC field |
220 | + allocate(pressure_bc_type(surface_element_count(p))) |
221 | + |
222 | + call get_entire_boundary_condition(p, (/"weakdirichlet", & |
223 | + "dirichlet "/), pressure_bc, pressure_bc_type) |
224 | + |
225 | + allocate(x_ele_bdy(x%dim,x%mesh%faces%shape%loc), & |
226 | + detwei_bdy(x_cvbdyshape%ngi), & |
227 | + normal_bdy(x%dim, x_cvbdyshape%ngi), & |
228 | + pressure_bc_val(pressure_bc%mesh%shape%loc)) |
229 | + |
230 | + allocate(u_nodes_bdy(u%mesh%faces%shape%loc)) |
231 | + |
232 | + allocate(ct_mat_local_bdy(x%dim, p%mesh%faces%shape%loc, u%mesh%faces%shape%loc)) |
233 | + |
234 | + ! Check if we need to multiply through by the non-linear volume fraction |
235 | + if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then |
236 | + multiphase = .true. |
237 | + |
238 | + vfrac => extract_scalar_field(state, "PhaseVolumeFraction") |
239 | + call allocate(nvfrac, vfrac%mesh, "NonlinearPhaseVolumeFraction") |
240 | + call zero(nvfrac) |
241 | + call get_nonlinear_volume_fraction(state, nvfrac) |
242 | + |
243 | + assert(surface_element_count(p) == surface_element_count(vfrac)) |
244 | + |
245 | + ewrite_minmax(nvfrac) |
246 | + |
247 | + ! If the Coordinate and PhaseVolumeFraction meshes are different, then we need to |
248 | + ! generate the PhaseVolumeFraction CV shape functions. |
249 | + if(.not.(nvfrac%mesh == x%mesh)) then |
250 | + nvfrac_cvbdyshape = make_cvbdy_element_shape(cvfaces, nvfrac%mesh%faces%shape) |
251 | + else |
252 | + nvfrac_cvbdyshape = x_cvbdyshape |
253 | + ! incease reference for deallocates below |
254 | + call incref(nvfrac_cvbdyshape) |
255 | + end if |
256 | + |
257 | + allocate(nvfrac_gi_f(nvfrac_cvbdyshape%ngi)) |
258 | + |
259 | + else |
260 | + multiphase = .false. |
261 | + nullify(vfrac) |
262 | + end if |
263 | + |
264 | + surface_element_loop: do sele = 1, surface_element_count(p) |
265 | + |
266 | + ! cycle if this not a dirichlet pressure BC. |
267 | + if (pressure_bc_type(sele) == 0) cycle |
268 | + |
269 | + ele = face_ele(x, sele) |
270 | + x_ele = ele_val(x, ele) |
271 | + x_ele_bdy = face_val(x, sele) |
272 | + u_nodes_bdy = face_global_nodes(u, sele) |
273 | + |
274 | + ! Get the phase volume fraction face value |
275 | + if(multiphase) then |
276 | + nvfrac_gi_f = face_val_at_quad(nvfrac, sele, nvfrac_cvbdyshape) |
277 | + end if |
278 | + |
279 | + call transform_cvsurf_facet_to_physical(x_ele, x_ele_bdy, & |
280 | + x_cvbdyshape, normal_bdy, detwei_bdy) |
281 | + |
282 | + ct_mat_local_bdy = 0.0 |
283 | + |
284 | + ! calculate the ct local matrix |
285 | + surface_nodal_loop_i: do iloc = 1, p%mesh%faces%shape%loc |
286 | + |
287 | + surface_face_loop: do face = 1, cvfaces%sfaces |
288 | + |
289 | + if(cvfaces%sneiloc(iloc,face) /= 0) then |
290 | + |
291 | + surface_quadrature_loop: do gi = 1, cvfaces%shape%ngi |
292 | + |
293 | + ggi = (face-1)*cvfaces%shape%ngi + gi |
294 | + |
295 | + surface_nodal_loop_j: do jloc = 1, u%mesh%faces%shape%loc |
296 | + |
297 | + surface_inner_dimension_loop: do dim = 1, size(normal_bdy,1) |
298 | + |
299 | + if(multiphase) then |
300 | + ct_mat_local_bdy(dim, iloc, jloc) = ct_mat_local_bdy(dim, iloc, jloc) + & |
301 | + u_cvbdyshape%n(jloc,ggi)*detwei_bdy(ggi)*nvfrac_gi_f(ggi)*normal_bdy(dim, ggi) |
302 | + else |
303 | + ct_mat_local_bdy(dim, iloc, jloc) = ct_mat_local_bdy(dim, iloc, jloc) + & |
304 | + u_cvbdyshape%n(jloc,ggi)*detwei_bdy(ggi)*normal_bdy(dim, ggi) |
305 | + end if |
306 | + |
307 | + end do surface_inner_dimension_loop |
308 | + |
309 | + end do surface_nodal_loop_j |
310 | + |
311 | + end do surface_quadrature_loop |
312 | + |
313 | + end if |
314 | + |
315 | + end do surface_face_loop |
316 | + |
317 | + end do surface_nodal_loop_i |
318 | + |
319 | + ! pressure dirichlet BC is -c*press_bc_val = -press_bc_val*ct, integrated over surface elements appropriate |
320 | + surface_outer_dimension_loop: do dim = 1, size(normal_bdy,1) |
321 | + |
322 | + ! for weak and strong pressure dirichlet bcs: |
323 | + ! / |
324 | + ! add -| N_i M_j \vec n p_j, where p_j are the prescribed bc values |
325 | + ! / |
326 | + |
327 | + call addto(mom_rhs, dim, u_nodes_bdy, -matmul(ele_val(pressure_bc, sele), ct_mat_local_bdy(dim,:,:))) |
328 | + |
329 | + end do surface_outer_dimension_loop |
330 | + |
331 | + end do surface_element_loop |
332 | + |
333 | + call deallocate(pressure_bc) |
334 | + deallocate(pressure_bc_type) |
335 | + deallocate(pressure_bc_val) |
336 | + call deallocate(x_cvbdyshape) |
337 | + call deallocate(u_cvbdyshape) |
338 | + deallocate(x_ele_bdy, detwei_bdy, normal_bdy) |
339 | + |
340 | + deallocate(u_nodes_bdy) |
341 | + deallocate(ct_mat_local_bdy) |
342 | + |
343 | + call deallocate(cvfaces) |
344 | + deallocate(x_ele) |
345 | + |
346 | + if(multiphase) then |
347 | + call deallocate(nvfrac) |
348 | + deallocate(nvfrac_gi_f) |
349 | + call deallocate(nvfrac_cvbdyshape) |
350 | + end if |
351 | + |
352 | + end subroutine add_pressure_dirichlet_bcs_cv |
353 | + |
354 | + ! -------------------------------------------------------------------------------------- |
355 | + |
356 | +end module pressure_dirichlet_bcs_cv |
357 | + |
358 | |
359 | === added directory 'tests/darcy_p0p1_test_cty_cv_pressBCinlet' |
360 | === added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/Makefile' |
361 | --- tests/darcy_p0p1_test_cty_cv_pressBCinlet/Makefile 1970-01-01 00:00:00 +0000 |
362 | +++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/Makefile 2011-12-23 12:05:25 +0000 |
363 | @@ -0,0 +1,32 @@ |
364 | +SHELL = sh |
365 | + |
366 | +SIM = darcy_p0p1_test_cty_cv_pressBCinlet |
367 | + |
368 | +# options for 1d mesh interval script |
369 | +MESH_SIZE = 100.0 |
370 | +LEFT_X = 0.0 |
371 | +RIGHT_X = 300.0 |
372 | + |
373 | +input: clean |
374 | + interval --dx=$(MESH_SIZE) $(LEFT_X) $(RIGHT_X) $(SIM)_1d |
375 | + gmsh -2 square.geo -o square.msh |
376 | + gmsh -3 cube.geo -o cube.msh |
377 | + |
378 | +clean: |
379 | + rm -f *.ele |
380 | + rm -f *.node |
381 | + rm -f *.bound |
382 | + rm -f *.edge |
383 | + rm -f *.face |
384 | + rm -f *.vtu |
385 | + rm -f *.pvtu |
386 | + rm -f *.s |
387 | + rm -f *.stat |
388 | + rm -f *.log-0 |
389 | + rm -f *.err-0 |
390 | + rm -f *.msh |
391 | + rm -f *.halo |
392 | + rm -f fluidity.err* |
393 | + rm -f fluidity.log* |
394 | + rm -f matrixdump* |
395 | + rm -f first_timestep_adapted_mesh* |
396 | |
397 | === added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/cube.geo' |
398 | --- tests/darcy_p0p1_test_cty_cv_pressBCinlet/cube.geo 1970-01-01 00:00:00 +0000 |
399 | +++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/cube.geo 2011-12-23 12:05:25 +0000 |
400 | @@ -0,0 +1,39 @@ |
401 | +// define a layer variable |
402 | +lay = 5; |
403 | + |
404 | +// define a len variable |
405 | +len = 300.0; |
406 | + |
407 | +Point(1) = {0.0, 0.0, 0.0, 1.0}; |
408 | + |
409 | +Extrude {len, 0.0, 0.0} { |
410 | + Point{1}; Layers{lay}; |
411 | +} |
412 | + |
413 | +Extrude {0.0, len, 0.0} { |
414 | + Line{1}; Layers{lay}; |
415 | +} |
416 | + |
417 | +Extrude {0.0, 0.0, len} { |
418 | + Surface{5}; Layers{lay}; |
419 | +} |
420 | + |
421 | +// bottom |
422 | +Physical Surface(1) = {5}; |
423 | + |
424 | +// top |
425 | +Physical Surface(2) = {27}; |
426 | + |
427 | +// left |
428 | +Physical Surface(3) = {26}; |
429 | + |
430 | +// right |
431 | +Physical Surface(4) = {18}; |
432 | + |
433 | +// front |
434 | +Physical Surface(5) = {14}; |
435 | + |
436 | +// back |
437 | +Physical Surface(6) = {22}; |
438 | + |
439 | +Physical Volume(1) = {1}; |
440 | |
441 | === added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet.xml' |
442 | --- tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet.xml 1970-01-01 00:00:00 +0000 |
443 | +++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet.xml 2011-12-23 12:05:25 +0000 |
444 | @@ -0,0 +1,118 @@ |
445 | +<?xml version="1.0" encoding="UTF-8" ?> |
446 | +<!DOCTYPE testproblem SYSTEM "regressiontest.dtd"> |
447 | + |
448 | +<testproblem> |
449 | + <name>darcy_p0p1_test_cty_cv_pressBCinlet</name> |
450 | + <owner userid="btollit"/> |
451 | + <tags>flml darcy</tags> |
452 | + <problem_definition length="short" nprocs="1"> |
453 | + <command_line> |
454 | +../../bin/fluidity darcy_p0p1_test_cty_cv_pressBCinlet_1d.flml |
455 | +../../bin/fluidity darcy_p0p1_test_cty_cv_pressBCinlet_2d.flml |
456 | +../../bin/fluidity darcy_p0p1_test_cty_cv_pressBCinlet_3d.flml |
457 | + </command_line> |
458 | + <!-- One/two/three dimensional problem for darcy flow with one region and one material using p0p1 element type testing the pressure gradient against analytic and the interstitial velocity. This tests using a diagnostic velocity absorption field associated with a mesh and momentum dg for darcy flows.--> |
459 | + </problem_definition> |
460 | + <variables> |
461 | + <variable name="pressure_1d" language="python"> |
462 | +import vtktools |
463 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_1d_1.vtu") |
464 | +pressure_1d = v.GetScalarRange("Pressure") |
465 | + </variable> |
466 | + <variable name="inter_vel_1d" language="python"> |
467 | +import vtktools |
468 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_1d_1.vtu") |
469 | +inter_vel_1d = max(v.GetScalarRange("interstitial_velocity")) |
470 | + </variable> |
471 | + <variable name="max_div_vel_1d" language="python"> |
472 | +import vtktools |
473 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_1d_1.vtu") |
474 | +max_div_vel_1d = max(v.GetScalarRange("ControlVolumeDivergence")) |
475 | + </variable> |
476 | + <variable name="pressure_2d" language="python"> |
477 | +import vtktools |
478 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_2d_1.vtu") |
479 | +pressure_2d = v.GetScalarRange("Pressure") |
480 | + </variable> |
481 | + <variable name="inter_vel_2d" language="python"> |
482 | +import vtktools |
483 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_2d_1.vtu") |
484 | +inter_vel_2d = max(v.GetScalarRange("interstitial_velocity")) |
485 | + </variable> |
486 | + <variable name="max_div_vel_2d" language="python"> |
487 | +import vtktools |
488 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_2d_1.vtu") |
489 | +max_div_vel_2d = max(v.GetScalarRange("ControlVolumeDivergence")) |
490 | + </variable> |
491 | + <variable name="pressure_3d" language="python"> |
492 | +import vtktools |
493 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_3d_1.vtu") |
494 | +pressure_3d = v.GetScalarRange("Pressure") |
495 | + </variable> |
496 | + <variable name="inter_vel_3d" language="python"> |
497 | +import vtktools |
498 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_3d_1.vtu") |
499 | +inter_vel_3d = max(v.GetScalarRange("interstitial_velocity")) |
500 | + </variable> |
501 | + <variable name="max_div_vel_3d" language="python"> |
502 | +import vtktools |
503 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_3d_1.vtu") |
504 | +max_div_vel_3d = max(v.GetScalarRange("ControlVolumeDivergence")) |
505 | + </variable> |
506 | + </variables> |
507 | + <pass_tests> |
508 | + <test name="change_P for 1d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python"> |
509 | +change_P = abs(max(pressure_1d) - min(pressure_1d)) |
510 | +visc = 1.0e-04 |
511 | +darcy_vel_BC = 5.0 |
512 | +perm = 1.0e-10 |
513 | +domain_length = 300.0 |
514 | +print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) |
515 | +assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) < 1.0e-09 |
516 | + </test> |
517 | + <test name="interstitial velocity 1d should equal darcy_vel/porosity), check relative difference to be under tolerance 1.0e-09" language="python"> |
518 | +analytic_vel = 10.0 |
519 | +print 'Relative error of interstitial velocity: ',abs((inter_vel_1d - analytic_vel)/analytic_vel) |
520 | +assert abs((inter_vel_1d - analytic_vel)/analytic_vel) < 1.0e-09 |
521 | + </test> |
522 | + <test name="max_div_vel_1d under tolerance 1.0e-09" language="python"> |
523 | +assert max_div_vel_1d < 1.0e-09 |
524 | + </test> |
525 | + <test name="change_P for 2d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python"> |
526 | +change_P = abs(max(pressure_2d) - min(pressure_2d)) |
527 | +visc = 1.0e-04 |
528 | +darcy_vel_BC = 5.0 |
529 | +perm = 1.0e-10 |
530 | +domain_length = 300.0 |
531 | +print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) |
532 | +assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) < 1.0e-09 |
533 | + </test> |
534 | + <test name="interstitial velocity 2d should equal darcy_vel/porosity), check relative difference to be under tolerance 1.0e-09" language="python"> |
535 | +analytic_vel = 10.0 |
536 | +print 'Relative error of interstitial velocity: ',abs((inter_vel_2d - analytic_vel)/analytic_vel) |
537 | +assert abs((inter_vel_2d - analytic_vel)/analytic_vel) < 1.0e-09 |
538 | + </test> |
539 | + <test name="max_div_vel_2d under tolerance 1.0e-09" language="python"> |
540 | +assert max_div_vel_2d < 1.0e-09 |
541 | + </test> |
542 | + <test name="change_P for 3d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python"> |
543 | +change_P = abs(max(pressure_3d) - min(pressure_3d)) |
544 | +visc = 1.0e-04 |
545 | +darcy_vel_BC = 5.0 |
546 | +perm = 1.0e-10 |
547 | +domain_length = 300.0 |
548 | +print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) |
549 | +assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) < 1.0e-09 |
550 | + </test> |
551 | + <test name="interstitial velocity 3d should equal darcy_vel/porosity), check relative difference to be under tolerance 1.0e-09" language="python"> |
552 | +analytic_vel = 10.0 |
553 | +print 'Relative error of interstitial velocity: ',abs((inter_vel_3d - analytic_vel)/analytic_vel) |
554 | +assert abs((inter_vel_3d - analytic_vel)/analytic_vel) < 1.0e-09 |
555 | + </test> |
556 | + <test name="max_div_vel_3d under tolerance 1.0e-09" language="python"> |
557 | +assert max_div_vel_3d < 1.0e-09 |
558 | + </test> |
559 | + </pass_tests> |
560 | + <warn_tests> |
561 | + </warn_tests> |
562 | +</testproblem> |
563 | |
564 | === added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_1d.flml' |
565 | --- tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_1d.flml 1970-01-01 00:00:00 +0000 |
566 | +++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_1d.flml 2011-12-23 12:05:25 +0000 |
567 | @@ -0,0 +1,411 @@ |
568 | +<?xml version='1.0' encoding='utf-8'?> |
569 | +<fluidity_options> |
570 | + <simulation_name> |
571 | + <string_value lines="1">darcy_p0p1_test_cty_cv_pressBCinlet_1d</string_value> |
572 | + <comment>a simple short test case for darcy flow using the element type p0p1 for velocity-pressure |
573 | + |
574 | +with the continuity equation tested with the control volume dual mesh |
575 | + |
576 | +the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution) |
577 | + |
578 | +it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct |
579 | + |
580 | +the darcy flow equation is |
581 | + |
582 | +sigma*darcy_vel = - grad P |
583 | + |
584 | +sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options) |
585 | + |
586 | +python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity |
587 | + |
588 | +this case has a 1 region 1 material with pressure boundary condition inflow. |
589 | + |
590 | +to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input. |
591 | + |
592 | +the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) |
593 | + |
594 | +the geometry is 1d and one time step is done (as nothing depends on time) |
595 | + |
596 | +the left and right boundary conditions of pressure are applied strongly. |
597 | + |
598 | +For 2d and 3d a second approach is also done via solving a pressure poisson equation then calculating the |
599 | +velocity diagnostically from this. This uses generic scalar and vector fields and python diagnostics. |
600 | +The latter didnt work in 1d. The same vel-press element pair is used as the mixed projection formulation.</comment> |
601 | + </simulation_name> |
602 | + <problem_type> |
603 | + <string_value lines="1">fluids</string_value> |
604 | + </problem_type> |
605 | + <geometry> |
606 | + <dimension> |
607 | + <integer_value rank="0">1</integer_value> |
608 | + </dimension> |
609 | + <mesh name="CoordinateMesh"> |
610 | + <from_file file_name="darcy_p0p1_test_cty_cv_pressBCinlet_1d"> |
611 | + <format name="triangle"/> |
612 | + <stat> |
613 | + <include_in_stat/> |
614 | + </stat> |
615 | + </from_file> |
616 | + </mesh> |
617 | + <mesh name="VelocityMesh"> |
618 | + <from_mesh> |
619 | + <mesh name="CoordinateMesh"/> |
620 | + <mesh_shape> |
621 | + <polynomial_degree> |
622 | + <integer_value rank="0">0</integer_value> |
623 | + </polynomial_degree> |
624 | + </mesh_shape> |
625 | + <mesh_continuity> |
626 | + <string_value>discontinuous</string_value> |
627 | + </mesh_continuity> |
628 | + <stat> |
629 | + <exclude_from_stat/> |
630 | + </stat> |
631 | + </from_mesh> |
632 | + </mesh> |
633 | + <mesh name="PressureMesh"> |
634 | + <from_mesh> |
635 | + <mesh name="CoordinateMesh"/> |
636 | + <mesh_shape> |
637 | + <polynomial_degree> |
638 | + <integer_value rank="0">1</integer_value> |
639 | + </polynomial_degree> |
640 | + </mesh_shape> |
641 | + <stat> |
642 | + <exclude_from_stat/> |
643 | + </stat> |
644 | + </from_mesh> |
645 | + </mesh> |
646 | + <mesh name="DGP0"> |
647 | + <from_mesh> |
648 | + <mesh name="CoordinateMesh"/> |
649 | + <mesh_shape> |
650 | + <polynomial_degree> |
651 | + <integer_value rank="0">0</integer_value> |
652 | + </polynomial_degree> |
653 | + </mesh_shape> |
654 | + <mesh_continuity> |
655 | + <string_value>discontinuous</string_value> |
656 | + </mesh_continuity> |
657 | + <stat> |
658 | + <exclude_from_stat/> |
659 | + </stat> |
660 | + </from_mesh> |
661 | + </mesh> |
662 | + <mesh name="OutputMesh"> |
663 | + <from_mesh> |
664 | + <mesh name="CoordinateMesh"/> |
665 | + <mesh_continuity> |
666 | + <string_value>discontinuous</string_value> |
667 | + </mesh_continuity> |
668 | + <stat> |
669 | + <exclude_from_stat/> |
670 | + </stat> |
671 | + </from_mesh> |
672 | + </mesh> |
673 | + <quadrature> |
674 | + <degree> |
675 | + <integer_value rank="0">3</integer_value> |
676 | + </degree> |
677 | + </quadrature> |
678 | + </geometry> |
679 | + <io> |
680 | + <dump_format> |
681 | + <string_value>vtk</string_value> |
682 | + </dump_format> |
683 | + <dump_period_in_timesteps> |
684 | + <constant> |
685 | + <integer_value rank="0">0</integer_value> |
686 | + </constant> |
687 | + </dump_period_in_timesteps> |
688 | + <output_mesh name="OutputMesh"/> |
689 | + <stat/> |
690 | + </io> |
691 | + <timestepping> |
692 | + <current_time> |
693 | + <real_value rank="0">0.0</real_value> |
694 | + </current_time> |
695 | + <timestep> |
696 | + <real_value rank="0">1.0</real_value> |
697 | + </timestep> |
698 | + <finish_time> |
699 | + <real_value rank="0">1.0</real_value> |
700 | + </finish_time> |
701 | + </timestepping> |
702 | + <material_phase name="fluid"> |
703 | + <scalar_field name="Pressure" rank="0"> |
704 | + <prognostic> |
705 | + <mesh name="PressureMesh"/> |
706 | + <spatial_discretisation> |
707 | + <continuous_galerkin> |
708 | + <test_continuity_with_cv_dual/> |
709 | + </continuous_galerkin> |
710 | + </spatial_discretisation> |
711 | + <scheme> |
712 | + <poisson_pressure_solution> |
713 | + <string_value lines="1">never</string_value> |
714 | + </poisson_pressure_solution> |
715 | + <use_projection_method/> |
716 | + </scheme> |
717 | + <solver> |
718 | + <iterative_method name="cg"/> |
719 | + <preconditioner name="sor"/> |
720 | + <relative_error> |
721 | + <real_value rank="0">1.0e-10</real_value> |
722 | + </relative_error> |
723 | + <max_iterations> |
724 | + <integer_value rank="0">1000</integer_value> |
725 | + </max_iterations> |
726 | + <never_ignore_solver_failures/> |
727 | + <diagnostics> |
728 | + <monitors/> |
729 | + </diagnostics> |
730 | + </solver> |
731 | + <boundary_conditions name="right_outflow"> |
732 | + <surface_ids> |
733 | + <integer_value shape="1" rank="1">2</integer_value> |
734 | + </surface_ids> |
735 | + <type name="dirichlet"> |
736 | + <constant> |
737 | + <real_value rank="0">0.0</real_value> |
738 | + </constant> |
739 | + </type> |
740 | + </boundary_conditions> |
741 | + <boundary_conditions name="left_inflow"> |
742 | + <surface_ids> |
743 | + <integer_value shape="1" rank="1">1</integer_value> |
744 | + </surface_ids> |
745 | + <type name="dirichlet"> |
746 | + <constant> |
747 | + <real_value rank="0">1.5e+09</real_value> |
748 | + </constant> |
749 | + </type> |
750 | + </boundary_conditions> |
751 | + <output/> |
752 | + <stat/> |
753 | + <convergence> |
754 | + <include_in_convergence/> |
755 | + </convergence> |
756 | + <detectors> |
757 | + <exclude_from_detectors/> |
758 | + </detectors> |
759 | + <steady_state> |
760 | + <include_in_steady_state/> |
761 | + </steady_state> |
762 | + <no_interpolation/> |
763 | + </prognostic> |
764 | + </scalar_field> |
765 | + <vector_field name="Velocity" rank="1"> |
766 | + <prognostic> |
767 | + <mesh name="VelocityMesh"/> |
768 | + <equation name="Boussinesq"/> |
769 | + <spatial_discretisation> |
770 | + <discontinuous_galerkin> |
771 | + <mass_terms> |
772 | + <exclude_mass_terms/> |
773 | + </mass_terms> |
774 | + <viscosity_scheme> |
775 | + <compact_discontinuous_galerkin/> |
776 | + </viscosity_scheme> |
777 | + <advection_scheme> |
778 | + <none/> |
779 | + <integrate_advection_by_parts> |
780 | + <twice/> |
781 | + </integrate_advection_by_parts> |
782 | + </advection_scheme> |
783 | + </discontinuous_galerkin> |
784 | + <conservative_advection> |
785 | + <real_value rank="0">0.0</real_value> |
786 | + </conservative_advection> |
787 | + </spatial_discretisation> |
788 | + <temporal_discretisation> |
789 | + <theta> |
790 | + <real_value rank="0">1.0</real_value> |
791 | + </theta> |
792 | + <relaxation> |
793 | + <real_value rank="0">1.0</real_value> |
794 | + </relaxation> |
795 | + </temporal_discretisation> |
796 | + <solver> |
797 | + <iterative_method name="gmres"> |
798 | + <restart> |
799 | + <integer_value rank="0">30</integer_value> |
800 | + </restart> |
801 | + </iterative_method> |
802 | + <preconditioner name="sor"/> |
803 | + <relative_error> |
804 | + <real_value rank="0">1.0e-10</real_value> |
805 | + </relative_error> |
806 | + <max_iterations> |
807 | + <integer_value rank="0">1000</integer_value> |
808 | + </max_iterations> |
809 | + <never_ignore_solver_failures/> |
810 | + <diagnostics> |
811 | + <monitors/> |
812 | + </diagnostics> |
813 | + </solver> |
814 | + <initial_condition name="WholeMesh"> |
815 | + <constant> |
816 | + <real_value shape="1" dim1="dim" rank="1">0.0</real_value> |
817 | + </constant> |
818 | + </initial_condition> |
819 | + <vector_field name="Absorption" rank="1"> |
820 | + <diagnostic> |
821 | + <mesh name="DGP0"/> |
822 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
823 | + <string_value lines="20" type="python"># get the prescribed fields |
824 | +perm = states["fluid"].scalar_fields["Permeability"] |
825 | +visc = states["fluid"].scalar_fields["Viscosity"] |
826 | + |
827 | +# loop the DMmesh_p0 nodes (which are element centred) |
828 | +for n in range(field.node_count): |
829 | + perm_n = perm.node_val(n) |
830 | + visc_n = visc.node_val(n) |
831 | + |
832 | + # calc the absorption term |
833 | + sigma_n = visc_n/perm_n |
834 | + |
835 | + # set the absorption term |
836 | + field.set(n,sigma_n)</string_value> |
837 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
838 | + |
839 | +also assume that the velocity is dg</comment> |
840 | + </algorithm> |
841 | + <output/> |
842 | + <stat> |
843 | + <include_in_stat/> |
844 | + </stat> |
845 | + <convergence> |
846 | + <include_in_convergence/> |
847 | + </convergence> |
848 | + <detectors> |
849 | + <include_in_detectors/> |
850 | + </detectors> |
851 | + <steady_state> |
852 | + <include_in_steady_state/> |
853 | + </steady_state> |
854 | + </diagnostic> |
855 | + <include_pressure_correction/> |
856 | + </vector_field> |
857 | + <output/> |
858 | + <stat> |
859 | + <include_in_stat/> |
860 | + <previous_time_step> |
861 | + <exclude_from_stat/> |
862 | + </previous_time_step> |
863 | + <nonlinear_field> |
864 | + <exclude_from_stat/> |
865 | + </nonlinear_field> |
866 | + </stat> |
867 | + <convergence> |
868 | + <include_in_convergence/> |
869 | + </convergence> |
870 | + <detectors> |
871 | + <include_in_detectors/> |
872 | + </detectors> |
873 | + <steady_state> |
874 | + <include_in_steady_state/> |
875 | + </steady_state> |
876 | + <consistent_interpolation/> |
877 | + </prognostic> |
878 | + </vector_field> |
879 | + <scalar_field name="Porosity" rank="0"> |
880 | + <prescribed> |
881 | + <mesh name="DGP0"/> |
882 | + <value name="WholeMesh"> |
883 | + <constant> |
884 | + <real_value rank="0">0.5</real_value> |
885 | + </constant> |
886 | + </value> |
887 | + <output/> |
888 | + <stat/> |
889 | + <detectors> |
890 | + <exclude_from_detectors/> |
891 | + </detectors> |
892 | + </prescribed> |
893 | + </scalar_field> |
894 | + <scalar_field name="Permeability" rank="0"> |
895 | + <prescribed> |
896 | + <mesh name="DGP0"/> |
897 | + <value name="WholeMesh"> |
898 | + <constant> |
899 | + <real_value rank="0">1.0e-10</real_value> |
900 | + </constant> |
901 | + </value> |
902 | + <output/> |
903 | + <stat/> |
904 | + <detectors> |
905 | + <exclude_from_detectors/> |
906 | + </detectors> |
907 | + </prescribed> |
908 | + </scalar_field> |
909 | + <scalar_field name="Viscosity" rank="0"> |
910 | + <prescribed> |
911 | + <mesh name="DGP0"/> |
912 | + <value name="WholeMesh"> |
913 | + <constant> |
914 | + <real_value rank="0">1.0e-04</real_value> |
915 | + </constant> |
916 | + </value> |
917 | + <output/> |
918 | + <stat/> |
919 | + <detectors> |
920 | + <exclude_from_detectors/> |
921 | + </detectors> |
922 | + </prescribed> |
923 | + </scalar_field> |
924 | + <scalar_field name="ControlVolumeDivergence" rank="0"> |
925 | + <diagnostic field_name="Velocity"> |
926 | + <algorithm name="Internal" material_phase_support="multiple"/> |
927 | + <mesh name="PressureMesh"/> |
928 | + <output/> |
929 | + <stat/> |
930 | + <convergence> |
931 | + <include_in_convergence/> |
932 | + </convergence> |
933 | + <detectors> |
934 | + <include_in_detectors/> |
935 | + </detectors> |
936 | + <steady_state> |
937 | + <include_in_steady_state/> |
938 | + </steady_state> |
939 | + </diagnostic> |
940 | + </scalar_field> |
941 | + <vector_field name="interstitial_velocity" rank="1"> |
942 | + <diagnostic> |
943 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
944 | + <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"] |
945 | +darcy_vel = states["fluid"].vector_fields["Velocity"] |
946 | + |
947 | +for ele in range(field.element_count): |
948 | + phi_ele = phi.node_val(ele) |
949 | + |
950 | + factor_ele = 1.0/max(phi_ele,1.0e-08) |
951 | + |
952 | + vel_ele = [] |
953 | + for i in range(len(darcy_vel.ele_val(ele))): |
954 | + vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i]) |
955 | + |
956 | + field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value> |
957 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
958 | + |
959 | +also assume that the velocity is dg</comment> |
960 | + </algorithm> |
961 | + <mesh name="VelocityMesh"/> |
962 | + <output/> |
963 | + <stat> |
964 | + <include_in_stat/> |
965 | + </stat> |
966 | + <convergence> |
967 | + <include_in_convergence/> |
968 | + </convergence> |
969 | + <detectors> |
970 | + <include_in_detectors/> |
971 | + </detectors> |
972 | + <steady_state> |
973 | + <include_in_steady_state/> |
974 | + </steady_state> |
975 | + </diagnostic> |
976 | + </vector_field> |
977 | + </material_phase> |
978 | +</fluidity_options> |
979 | |
980 | === added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_2d.flml' |
981 | --- tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_2d.flml 1970-01-01 00:00:00 +0000 |
982 | +++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_2d.flml 2011-12-23 12:05:25 +0000 |
983 | @@ -0,0 +1,427 @@ |
984 | +<?xml version='1.0' encoding='utf-8'?> |
985 | +<fluidity_options> |
986 | + <simulation_name> |
987 | + <string_value lines="1">darcy_p0p1_test_cty_cv_pressBCinlet_2d</string_value> |
988 | + <comment>a simple short test case for darcy flow using the element type p0p1 for velocity-pressure |
989 | + |
990 | +with the continuity equation tested with the control volume dual mesh |
991 | + |
992 | +the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution) |
993 | + |
994 | +it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct |
995 | + |
996 | +the darcy flow equation is |
997 | + |
998 | +sigma*darcy_vel = - grad P |
999 | + |
1000 | +sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options) |
1001 | + |
1002 | +python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity |
1003 | + |
1004 | +this case has a 1 region 1 material with pressure boundary condition inflow. |
1005 | + |
1006 | +to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input. |
1007 | + |
1008 | +the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) |
1009 | + |
1010 | +the geometry is 2d (although this is a 1d problem) and one time step is done (as nothing depends on time) |
1011 | + |
1012 | +the left and right boundary conditions of pressure are applied strongly. |
1013 | + |
1014 | +For 2d and 3d a second approach is also done via solving a pressure poisson equation then calculating the |
1015 | +velocity diagnostically from this. This uses generic scalar and vector fields and python diagnostics. |
1016 | +The latter didnt work in 1d. The same vel-press element pair is used as the mixed projection formulation.</comment> |
1017 | + </simulation_name> |
1018 | + <problem_type> |
1019 | + <string_value lines="1">fluids</string_value> |
1020 | + </problem_type> |
1021 | + <geometry> |
1022 | + <dimension> |
1023 | + <integer_value rank="0">2</integer_value> |
1024 | + </dimension> |
1025 | + <mesh name="CoordinateMesh"> |
1026 | + <from_file file_name="square"> |
1027 | + <format name="gmsh"/> |
1028 | + <stat> |
1029 | + <include_in_stat/> |
1030 | + </stat> |
1031 | + </from_file> |
1032 | + </mesh> |
1033 | + <mesh name="VelocityMesh"> |
1034 | + <from_mesh> |
1035 | + <mesh name="CoordinateMesh"/> |
1036 | + <mesh_shape> |
1037 | + <polynomial_degree> |
1038 | + <integer_value rank="0">0</integer_value> |
1039 | + </polynomial_degree> |
1040 | + </mesh_shape> |
1041 | + <mesh_continuity> |
1042 | + <string_value>discontinuous</string_value> |
1043 | + </mesh_continuity> |
1044 | + <stat> |
1045 | + <exclude_from_stat/> |
1046 | + </stat> |
1047 | + </from_mesh> |
1048 | + </mesh> |
1049 | + <mesh name="PressureMesh"> |
1050 | + <from_mesh> |
1051 | + <mesh name="CoordinateMesh"/> |
1052 | + <mesh_shape> |
1053 | + <polynomial_degree> |
1054 | + <integer_value rank="0">1</integer_value> |
1055 | + </polynomial_degree> |
1056 | + </mesh_shape> |
1057 | + <stat> |
1058 | + <exclude_from_stat/> |
1059 | + </stat> |
1060 | + </from_mesh> |
1061 | + </mesh> |
1062 | + <mesh name="DGP0"> |
1063 | + <from_mesh> |
1064 | + <mesh name="CoordinateMesh"/> |
1065 | + <mesh_shape> |
1066 | + <polynomial_degree> |
1067 | + <integer_value rank="0">0</integer_value> |
1068 | + </polynomial_degree> |
1069 | + </mesh_shape> |
1070 | + <mesh_continuity> |
1071 | + <string_value>discontinuous</string_value> |
1072 | + </mesh_continuity> |
1073 | + <stat> |
1074 | + <exclude_from_stat/> |
1075 | + </stat> |
1076 | + </from_mesh> |
1077 | + </mesh> |
1078 | + <mesh name="OutputMesh"> |
1079 | + <from_mesh> |
1080 | + <mesh name="CoordinateMesh"/> |
1081 | + <mesh_continuity> |
1082 | + <string_value>discontinuous</string_value> |
1083 | + </mesh_continuity> |
1084 | + <stat> |
1085 | + <exclude_from_stat/> |
1086 | + </stat> |
1087 | + </from_mesh> |
1088 | + </mesh> |
1089 | + <quadrature> |
1090 | + <degree> |
1091 | + <integer_value rank="0">3</integer_value> |
1092 | + </degree> |
1093 | + </quadrature> |
1094 | + </geometry> |
1095 | + <io> |
1096 | + <dump_format> |
1097 | + <string_value>vtk</string_value> |
1098 | + </dump_format> |
1099 | + <dump_period_in_timesteps> |
1100 | + <constant> |
1101 | + <integer_value rank="0">0</integer_value> |
1102 | + </constant> |
1103 | + </dump_period_in_timesteps> |
1104 | + <output_mesh name="OutputMesh"/> |
1105 | + <stat/> |
1106 | + </io> |
1107 | + <timestepping> |
1108 | + <current_time> |
1109 | + <real_value rank="0">0.0</real_value> |
1110 | + </current_time> |
1111 | + <timestep> |
1112 | + <real_value rank="0">1.0</real_value> |
1113 | + </timestep> |
1114 | + <finish_time> |
1115 | + <real_value rank="0">1.0</real_value> |
1116 | + </finish_time> |
1117 | + </timestepping> |
1118 | + <material_phase name="fluid"> |
1119 | + <scalar_field name="Pressure" rank="0"> |
1120 | + <prognostic> |
1121 | + <mesh name="PressureMesh"/> |
1122 | + <spatial_discretisation> |
1123 | + <continuous_galerkin> |
1124 | + <test_continuity_with_cv_dual/> |
1125 | + </continuous_galerkin> |
1126 | + </spatial_discretisation> |
1127 | + <scheme> |
1128 | + <poisson_pressure_solution> |
1129 | + <string_value lines="1">never</string_value> |
1130 | + </poisson_pressure_solution> |
1131 | + <use_projection_method/> |
1132 | + </scheme> |
1133 | + <solver> |
1134 | + <iterative_method name="cg"/> |
1135 | + <preconditioner name="sor"/> |
1136 | + <relative_error> |
1137 | + <real_value rank="0">1.0e-10</real_value> |
1138 | + </relative_error> |
1139 | + <max_iterations> |
1140 | + <integer_value rank="0">1000</integer_value> |
1141 | + </max_iterations> |
1142 | + <never_ignore_solver_failures/> |
1143 | + <diagnostics> |
1144 | + <monitors/> |
1145 | + </diagnostics> |
1146 | + </solver> |
1147 | + <boundary_conditions name="right_outflow"> |
1148 | + <surface_ids> |
1149 | + <integer_value shape="1" rank="1">8</integer_value> |
1150 | + </surface_ids> |
1151 | + <type name="dirichlet"> |
1152 | + <constant> |
1153 | + <real_value rank="0">0.0</real_value> |
1154 | + </constant> |
1155 | + </type> |
1156 | + </boundary_conditions> |
1157 | + <boundary_conditions name="left_inflow"> |
1158 | + <surface_ids> |
1159 | + <integer_value shape="1" rank="1">10</integer_value> |
1160 | + </surface_ids> |
1161 | + <type name="dirichlet"> |
1162 | + <apply_weakly/> |
1163 | + <constant> |
1164 | + <real_value rank="0">1.5e+09</real_value> |
1165 | + </constant> |
1166 | + </type> |
1167 | + </boundary_conditions> |
1168 | + <output/> |
1169 | + <stat/> |
1170 | + <convergence> |
1171 | + <include_in_convergence/> |
1172 | + </convergence> |
1173 | + <detectors> |
1174 | + <exclude_from_detectors/> |
1175 | + </detectors> |
1176 | + <steady_state> |
1177 | + <include_in_steady_state/> |
1178 | + </steady_state> |
1179 | + <no_interpolation/> |
1180 | + </prognostic> |
1181 | + </scalar_field> |
1182 | + <vector_field name="Velocity" rank="1"> |
1183 | + <prognostic> |
1184 | + <mesh name="VelocityMesh"/> |
1185 | + <equation name="Boussinesq"/> |
1186 | + <spatial_discretisation> |
1187 | + <discontinuous_galerkin> |
1188 | + <mass_terms> |
1189 | + <exclude_mass_terms/> |
1190 | + </mass_terms> |
1191 | + <viscosity_scheme> |
1192 | + <compact_discontinuous_galerkin/> |
1193 | + </viscosity_scheme> |
1194 | + <advection_scheme> |
1195 | + <none/> |
1196 | + <integrate_advection_by_parts> |
1197 | + <twice/> |
1198 | + </integrate_advection_by_parts> |
1199 | + </advection_scheme> |
1200 | + </discontinuous_galerkin> |
1201 | + <conservative_advection> |
1202 | + <real_value rank="0">0.0</real_value> |
1203 | + </conservative_advection> |
1204 | + </spatial_discretisation> |
1205 | + <temporal_discretisation> |
1206 | + <theta> |
1207 | + <real_value rank="0">1.0</real_value> |
1208 | + </theta> |
1209 | + <relaxation> |
1210 | + <real_value rank="0">1.0</real_value> |
1211 | + </relaxation> |
1212 | + </temporal_discretisation> |
1213 | + <solver> |
1214 | + <iterative_method name="gmres"> |
1215 | + <restart> |
1216 | + <integer_value rank="0">30</integer_value> |
1217 | + </restart> |
1218 | + </iterative_method> |
1219 | + <preconditioner name="sor"/> |
1220 | + <relative_error> |
1221 | + <real_value rank="0">1.0e-10</real_value> |
1222 | + </relative_error> |
1223 | + <max_iterations> |
1224 | + <integer_value rank="0">1000</integer_value> |
1225 | + </max_iterations> |
1226 | + <never_ignore_solver_failures/> |
1227 | + <diagnostics> |
1228 | + <monitors/> |
1229 | + </diagnostics> |
1230 | + </solver> |
1231 | + <initial_condition name="WholeMesh"> |
1232 | + <constant> |
1233 | + <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value> |
1234 | + </constant> |
1235 | + </initial_condition> |
1236 | + <boundary_conditions name="no_outflow_top_bottom"> |
1237 | + <surface_ids> |
1238 | + <integer_value shape="2" rank="1">7 9</integer_value> |
1239 | + </surface_ids> |
1240 | + <type name="dirichlet"> |
1241 | + <apply_weakly/> |
1242 | + <align_bc_with_cartesian> |
1243 | + <y_component> |
1244 | + <constant> |
1245 | + <real_value rank="0">0.0</real_value> |
1246 | + </constant> |
1247 | + </y_component> |
1248 | + </align_bc_with_cartesian> |
1249 | + </type> |
1250 | + </boundary_conditions> |
1251 | + <vector_field name="Absorption" rank="1"> |
1252 | + <diagnostic> |
1253 | + <mesh name="DGP0"/> |
1254 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
1255 | + <string_value lines="20" type="python"># get the prescribed fields |
1256 | +perm = states["fluid"].scalar_fields["Permeability"] |
1257 | +visc = states["fluid"].scalar_fields["Viscosity"] |
1258 | + |
1259 | +# loop the DMmesh_p0 nodes (which are element centred) |
1260 | +for n in range(field.node_count): |
1261 | + perm_n = perm.node_val(n) |
1262 | + visc_n = visc.node_val(n) |
1263 | + |
1264 | + # calc the absorption term |
1265 | + sigma_n = visc_n/perm_n |
1266 | + |
1267 | + # set the absorption term |
1268 | + field.set(n,sigma_n)</string_value> |
1269 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
1270 | + |
1271 | +also assume that the velocity is dg</comment> |
1272 | + </algorithm> |
1273 | + <output/> |
1274 | + <stat> |
1275 | + <include_in_stat/> |
1276 | + </stat> |
1277 | + <convergence> |
1278 | + <include_in_convergence/> |
1279 | + </convergence> |
1280 | + <detectors> |
1281 | + <include_in_detectors/> |
1282 | + </detectors> |
1283 | + <steady_state> |
1284 | + <include_in_steady_state/> |
1285 | + </steady_state> |
1286 | + </diagnostic> |
1287 | + <include_pressure_correction/> |
1288 | + </vector_field> |
1289 | + <output/> |
1290 | + <stat> |
1291 | + <include_in_stat/> |
1292 | + <previous_time_step> |
1293 | + <exclude_from_stat/> |
1294 | + </previous_time_step> |
1295 | + <nonlinear_field> |
1296 | + <exclude_from_stat/> |
1297 | + </nonlinear_field> |
1298 | + </stat> |
1299 | + <convergence> |
1300 | + <include_in_convergence/> |
1301 | + </convergence> |
1302 | + <detectors> |
1303 | + <include_in_detectors/> |
1304 | + </detectors> |
1305 | + <steady_state> |
1306 | + <include_in_steady_state/> |
1307 | + </steady_state> |
1308 | + <consistent_interpolation/> |
1309 | + </prognostic> |
1310 | + </vector_field> |
1311 | + <scalar_field name="Porosity" rank="0"> |
1312 | + <prescribed> |
1313 | + <mesh name="DGP0"/> |
1314 | + <value name="WholeMesh"> |
1315 | + <constant> |
1316 | + <real_value rank="0">0.5</real_value> |
1317 | + </constant> |
1318 | + </value> |
1319 | + <output/> |
1320 | + <stat/> |
1321 | + <detectors> |
1322 | + <exclude_from_detectors/> |
1323 | + </detectors> |
1324 | + </prescribed> |
1325 | + </scalar_field> |
1326 | + <scalar_field name="Permeability" rank="0"> |
1327 | + <prescribed> |
1328 | + <mesh name="DGP0"/> |
1329 | + <value name="WholeMesh"> |
1330 | + <constant> |
1331 | + <real_value rank="0">1.0e-10</real_value> |
1332 | + </constant> |
1333 | + </value> |
1334 | + <output/> |
1335 | + <stat/> |
1336 | + <detectors> |
1337 | + <exclude_from_detectors/> |
1338 | + </detectors> |
1339 | + </prescribed> |
1340 | + </scalar_field> |
1341 | + <scalar_field name="Viscosity" rank="0"> |
1342 | + <prescribed> |
1343 | + <mesh name="DGP0"/> |
1344 | + <value name="WholeMesh"> |
1345 | + <constant> |
1346 | + <real_value rank="0">1.0e-04</real_value> |
1347 | + </constant> |
1348 | + </value> |
1349 | + <output/> |
1350 | + <stat/> |
1351 | + <detectors> |
1352 | + <exclude_from_detectors/> |
1353 | + </detectors> |
1354 | + </prescribed> |
1355 | + </scalar_field> |
1356 | + <scalar_field name="ControlVolumeDivergence" rank="0"> |
1357 | + <diagnostic field_name="Velocity"> |
1358 | + <algorithm name="Internal" material_phase_support="multiple"/> |
1359 | + <mesh name="PressureMesh"/> |
1360 | + <output/> |
1361 | + <stat/> |
1362 | + <convergence> |
1363 | + <include_in_convergence/> |
1364 | + </convergence> |
1365 | + <detectors> |
1366 | + <include_in_detectors/> |
1367 | + </detectors> |
1368 | + <steady_state> |
1369 | + <include_in_steady_state/> |
1370 | + </steady_state> |
1371 | + </diagnostic> |
1372 | + </scalar_field> |
1373 | + <vector_field name="interstitial_velocity" rank="1"> |
1374 | + <diagnostic> |
1375 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
1376 | + <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"] |
1377 | +darcy_vel = states["fluid"].vector_fields["Velocity"] |
1378 | + |
1379 | +for ele in range(field.element_count): |
1380 | + phi_ele = phi.node_val(ele) |
1381 | + |
1382 | + factor_ele = 1.0/max(phi_ele,1.0e-08) |
1383 | + |
1384 | + vel_ele = [] |
1385 | + for i in range(len(darcy_vel.ele_val(ele))): |
1386 | + vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i]) |
1387 | + |
1388 | + field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value> |
1389 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
1390 | + |
1391 | +also assume that the velocity is dg</comment> |
1392 | + </algorithm> |
1393 | + <mesh name="VelocityMesh"/> |
1394 | + <output/> |
1395 | + <stat> |
1396 | + <include_in_stat/> |
1397 | + </stat> |
1398 | + <convergence> |
1399 | + <include_in_convergence/> |
1400 | + </convergence> |
1401 | + <detectors> |
1402 | + <include_in_detectors/> |
1403 | + </detectors> |
1404 | + <steady_state> |
1405 | + <include_in_steady_state/> |
1406 | + </steady_state> |
1407 | + </diagnostic> |
1408 | + </vector_field> |
1409 | + </material_phase> |
1410 | +</fluidity_options> |
1411 | |
1412 | === added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_3d.flml' |
1413 | --- tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_3d.flml 1970-01-01 00:00:00 +0000 |
1414 | +++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/darcy_p0p1_test_cty_cv_pressBCinlet_3d.flml 2011-12-23 12:05:25 +0000 |
1415 | @@ -0,0 +1,442 @@ |
1416 | +<?xml version='1.0' encoding='utf-8'?> |
1417 | +<fluidity_options> |
1418 | + <simulation_name> |
1419 | + <string_value lines="1">darcy_p0p1_test_cty_cv_pressBCinlet_3d</string_value> |
1420 | + <comment>a simple short test case for darcy flow using the element type p0p1 for velocity-pressure |
1421 | + |
1422 | +with the continuity equation tested with the control volume dual mesh |
1423 | + |
1424 | +the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution) |
1425 | + |
1426 | +it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct |
1427 | + |
1428 | +the darcy flow equation is |
1429 | + |
1430 | +sigma*darcy_vel = - grad P |
1431 | + |
1432 | +sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options) |
1433 | + |
1434 | +python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity |
1435 | + |
1436 | +this case has a 1 region 1 material with pressure boundary condition inflow. |
1437 | + |
1438 | +to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input. |
1439 | + |
1440 | +the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) |
1441 | + |
1442 | +the geometry is 3d (although this is a 1d problem) and one time step is done (as nothing depends on time) |
1443 | + |
1444 | +the left and right boundary conditions of pressure are applied strongly. |
1445 | + |
1446 | +For 2d and 3d a second approach is also done via solving a pressure poisson equation then calculating the |
1447 | +velocity diagnostically from this. This uses generic scalar and vector fields and python diagnostics. |
1448 | +The latter didnt work in 1d. The same vel-press element pair is used as the mixed projection formulation.</comment> |
1449 | + </simulation_name> |
1450 | + <problem_type> |
1451 | + <string_value lines="1">fluids</string_value> |
1452 | + </problem_type> |
1453 | + <geometry> |
1454 | + <dimension> |
1455 | + <integer_value rank="0">3</integer_value> |
1456 | + </dimension> |
1457 | + <mesh name="CoordinateMesh"> |
1458 | + <from_file file_name="cube"> |
1459 | + <format name="gmsh"/> |
1460 | + <stat> |
1461 | + <include_in_stat/> |
1462 | + </stat> |
1463 | + </from_file> |
1464 | + </mesh> |
1465 | + <mesh name="VelocityMesh"> |
1466 | + <from_mesh> |
1467 | + <mesh name="CoordinateMesh"/> |
1468 | + <mesh_shape> |
1469 | + <polynomial_degree> |
1470 | + <integer_value rank="0">0</integer_value> |
1471 | + </polynomial_degree> |
1472 | + </mesh_shape> |
1473 | + <mesh_continuity> |
1474 | + <string_value>discontinuous</string_value> |
1475 | + </mesh_continuity> |
1476 | + <stat> |
1477 | + <exclude_from_stat/> |
1478 | + </stat> |
1479 | + </from_mesh> |
1480 | + </mesh> |
1481 | + <mesh name="PressureMesh"> |
1482 | + <from_mesh> |
1483 | + <mesh name="CoordinateMesh"/> |
1484 | + <mesh_shape> |
1485 | + <polynomial_degree> |
1486 | + <integer_value rank="0">1</integer_value> |
1487 | + </polynomial_degree> |
1488 | + </mesh_shape> |
1489 | + <stat> |
1490 | + <exclude_from_stat/> |
1491 | + </stat> |
1492 | + </from_mesh> |
1493 | + </mesh> |
1494 | + <mesh name="DGP0"> |
1495 | + <from_mesh> |
1496 | + <mesh name="CoordinateMesh"/> |
1497 | + <mesh_shape> |
1498 | + <polynomial_degree> |
1499 | + <integer_value rank="0">0</integer_value> |
1500 | + </polynomial_degree> |
1501 | + </mesh_shape> |
1502 | + <mesh_continuity> |
1503 | + <string_value>discontinuous</string_value> |
1504 | + </mesh_continuity> |
1505 | + <stat> |
1506 | + <exclude_from_stat/> |
1507 | + </stat> |
1508 | + </from_mesh> |
1509 | + </mesh> |
1510 | + <mesh name="OutputMesh"> |
1511 | + <from_mesh> |
1512 | + <mesh name="CoordinateMesh"/> |
1513 | + <mesh_continuity> |
1514 | + <string_value>discontinuous</string_value> |
1515 | + </mesh_continuity> |
1516 | + <stat> |
1517 | + <exclude_from_stat/> |
1518 | + </stat> |
1519 | + </from_mesh> |
1520 | + </mesh> |
1521 | + <quadrature> |
1522 | + <degree> |
1523 | + <integer_value rank="0">3</integer_value> |
1524 | + </degree> |
1525 | + </quadrature> |
1526 | + </geometry> |
1527 | + <io> |
1528 | + <dump_format> |
1529 | + <string_value>vtk</string_value> |
1530 | + </dump_format> |
1531 | + <dump_period_in_timesteps> |
1532 | + <constant> |
1533 | + <integer_value rank="0">0</integer_value> |
1534 | + </constant> |
1535 | + </dump_period_in_timesteps> |
1536 | + <output_mesh name="OutputMesh"/> |
1537 | + <stat/> |
1538 | + </io> |
1539 | + <timestepping> |
1540 | + <current_time> |
1541 | + <real_value rank="0">0.0</real_value> |
1542 | + </current_time> |
1543 | + <timestep> |
1544 | + <real_value rank="0">1.0</real_value> |
1545 | + </timestep> |
1546 | + <finish_time> |
1547 | + <real_value rank="0">1.0</real_value> |
1548 | + </finish_time> |
1549 | + </timestepping> |
1550 | + <material_phase name="fluid"> |
1551 | + <scalar_field name="Pressure" rank="0"> |
1552 | + <prognostic> |
1553 | + <mesh name="PressureMesh"/> |
1554 | + <spatial_discretisation> |
1555 | + <continuous_galerkin> |
1556 | + <test_continuity_with_cv_dual/> |
1557 | + </continuous_galerkin> |
1558 | + </spatial_discretisation> |
1559 | + <scheme> |
1560 | + <poisson_pressure_solution> |
1561 | + <string_value lines="1">never</string_value> |
1562 | + </poisson_pressure_solution> |
1563 | + <use_projection_method/> |
1564 | + </scheme> |
1565 | + <solver> |
1566 | + <iterative_method name="cg"/> |
1567 | + <preconditioner name="sor"/> |
1568 | + <relative_error> |
1569 | + <real_value rank="0">1.0e-10</real_value> |
1570 | + </relative_error> |
1571 | + <max_iterations> |
1572 | + <integer_value rank="0">1000</integer_value> |
1573 | + </max_iterations> |
1574 | + <never_ignore_solver_failures/> |
1575 | + <diagnostics> |
1576 | + <monitors/> |
1577 | + </diagnostics> |
1578 | + </solver> |
1579 | + <boundary_conditions name="right_outflow"> |
1580 | + <surface_ids> |
1581 | + <integer_value shape="1" rank="1">4</integer_value> |
1582 | + </surface_ids> |
1583 | + <type name="dirichlet"> |
1584 | + <constant> |
1585 | + <real_value rank="0">0.0</real_value> |
1586 | + </constant> |
1587 | + </type> |
1588 | + </boundary_conditions> |
1589 | + <boundary_conditions name="left_inflow"> |
1590 | + <surface_ids> |
1591 | + <integer_value shape="1" rank="1">3</integer_value> |
1592 | + </surface_ids> |
1593 | + <type name="dirichlet"> |
1594 | + <apply_weakly/> |
1595 | + <constant> |
1596 | + <real_value rank="0">1.5e+09</real_value> |
1597 | + </constant> |
1598 | + </type> |
1599 | + </boundary_conditions> |
1600 | + <output/> |
1601 | + <stat/> |
1602 | + <convergence> |
1603 | + <include_in_convergence/> |
1604 | + </convergence> |
1605 | + <detectors> |
1606 | + <exclude_from_detectors/> |
1607 | + </detectors> |
1608 | + <steady_state> |
1609 | + <include_in_steady_state/> |
1610 | + </steady_state> |
1611 | + <no_interpolation/> |
1612 | + </prognostic> |
1613 | + </scalar_field> |
1614 | + <vector_field name="Velocity" rank="1"> |
1615 | + <prognostic> |
1616 | + <mesh name="VelocityMesh"/> |
1617 | + <equation name="Boussinesq"/> |
1618 | + <spatial_discretisation> |
1619 | + <discontinuous_galerkin> |
1620 | + <mass_terms> |
1621 | + <exclude_mass_terms/> |
1622 | + </mass_terms> |
1623 | + <viscosity_scheme> |
1624 | + <compact_discontinuous_galerkin/> |
1625 | + </viscosity_scheme> |
1626 | + <advection_scheme> |
1627 | + <none/> |
1628 | + <integrate_advection_by_parts> |
1629 | + <twice/> |
1630 | + </integrate_advection_by_parts> |
1631 | + </advection_scheme> |
1632 | + </discontinuous_galerkin> |
1633 | + <conservative_advection> |
1634 | + <real_value rank="0">0.0</real_value> |
1635 | + </conservative_advection> |
1636 | + </spatial_discretisation> |
1637 | + <temporal_discretisation> |
1638 | + <theta> |
1639 | + <real_value rank="0">1.0</real_value> |
1640 | + </theta> |
1641 | + <relaxation> |
1642 | + <real_value rank="0">1.0</real_value> |
1643 | + </relaxation> |
1644 | + </temporal_discretisation> |
1645 | + <solver> |
1646 | + <iterative_method name="gmres"> |
1647 | + <restart> |
1648 | + <integer_value rank="0">30</integer_value> |
1649 | + </restart> |
1650 | + </iterative_method> |
1651 | + <preconditioner name="sor"/> |
1652 | + <relative_error> |
1653 | + <real_value rank="0">1.0e-10</real_value> |
1654 | + </relative_error> |
1655 | + <max_iterations> |
1656 | + <integer_value rank="0">1000</integer_value> |
1657 | + </max_iterations> |
1658 | + <never_ignore_solver_failures/> |
1659 | + <diagnostics> |
1660 | + <monitors/> |
1661 | + </diagnostics> |
1662 | + </solver> |
1663 | + <initial_condition name="WholeMesh"> |
1664 | + <constant> |
1665 | + <real_value shape="3" dim1="dim" rank="1">0.0 0.0 0.0</real_value> |
1666 | + </constant> |
1667 | + </initial_condition> |
1668 | + <boundary_conditions name="no_outflow_top_bottom"> |
1669 | + <surface_ids> |
1670 | + <integer_value shape="2" rank="1">1 2</integer_value> |
1671 | + </surface_ids> |
1672 | + <type name="dirichlet"> |
1673 | + <apply_weakly/> |
1674 | + <align_bc_with_cartesian> |
1675 | + <z_component> |
1676 | + <constant> |
1677 | + <real_value rank="0">0.0</real_value> |
1678 | + </constant> |
1679 | + </z_component> |
1680 | + </align_bc_with_cartesian> |
1681 | + </type> |
1682 | + </boundary_conditions> |
1683 | + <boundary_conditions name="no_outflow_front_back"> |
1684 | + <surface_ids> |
1685 | + <integer_value shape="2" rank="1">5 6</integer_value> |
1686 | + </surface_ids> |
1687 | + <type name="dirichlet"> |
1688 | + <apply_weakly/> |
1689 | + <align_bc_with_cartesian> |
1690 | + <y_component> |
1691 | + <constant> |
1692 | + <real_value rank="0">0.0</real_value> |
1693 | + </constant> |
1694 | + </y_component> |
1695 | + </align_bc_with_cartesian> |
1696 | + </type> |
1697 | + </boundary_conditions> |
1698 | + <vector_field name="Absorption" rank="1"> |
1699 | + <diagnostic> |
1700 | + <mesh name="DGP0"/> |
1701 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
1702 | + <string_value lines="20" type="python"># get the prescribed fields |
1703 | +perm = states["fluid"].scalar_fields["Permeability"] |
1704 | +visc = states["fluid"].scalar_fields["Viscosity"] |
1705 | + |
1706 | +# loop the DMmesh_p0 nodes (which are element centred) |
1707 | +for n in range(field.node_count): |
1708 | + perm_n = perm.node_val(n) |
1709 | + visc_n = visc.node_val(n) |
1710 | + |
1711 | + # calc the absorption term |
1712 | + sigma_n = visc_n/perm_n |
1713 | + |
1714 | + # set the absorption term |
1715 | + field.set(n,sigma_n)</string_value> |
1716 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
1717 | + |
1718 | +also assume that the velocity is dg</comment> |
1719 | + </algorithm> |
1720 | + <output/> |
1721 | + <stat> |
1722 | + <include_in_stat/> |
1723 | + </stat> |
1724 | + <convergence> |
1725 | + <include_in_convergence/> |
1726 | + </convergence> |
1727 | + <detectors> |
1728 | + <include_in_detectors/> |
1729 | + </detectors> |
1730 | + <steady_state> |
1731 | + <include_in_steady_state/> |
1732 | + </steady_state> |
1733 | + </diagnostic> |
1734 | + <include_pressure_correction/> |
1735 | + </vector_field> |
1736 | + <output/> |
1737 | + <stat> |
1738 | + <include_in_stat/> |
1739 | + <previous_time_step> |
1740 | + <exclude_from_stat/> |
1741 | + </previous_time_step> |
1742 | + <nonlinear_field> |
1743 | + <exclude_from_stat/> |
1744 | + </nonlinear_field> |
1745 | + </stat> |
1746 | + <convergence> |
1747 | + <include_in_convergence/> |
1748 | + </convergence> |
1749 | + <detectors> |
1750 | + <include_in_detectors/> |
1751 | + </detectors> |
1752 | + <steady_state> |
1753 | + <include_in_steady_state/> |
1754 | + </steady_state> |
1755 | + <consistent_interpolation/> |
1756 | + </prognostic> |
1757 | + </vector_field> |
1758 | + <scalar_field name="Porosity" rank="0"> |
1759 | + <prescribed> |
1760 | + <mesh name="DGP0"/> |
1761 | + <value name="WholeMesh"> |
1762 | + <constant> |
1763 | + <real_value rank="0">0.5</real_value> |
1764 | + </constant> |
1765 | + </value> |
1766 | + <output/> |
1767 | + <stat/> |
1768 | + <detectors> |
1769 | + <exclude_from_detectors/> |
1770 | + </detectors> |
1771 | + </prescribed> |
1772 | + </scalar_field> |
1773 | + <scalar_field name="Permeability" rank="0"> |
1774 | + <prescribed> |
1775 | + <mesh name="DGP0"/> |
1776 | + <value name="WholeMesh"> |
1777 | + <constant> |
1778 | + <real_value rank="0">1.0e-10</real_value> |
1779 | + </constant> |
1780 | + </value> |
1781 | + <output/> |
1782 | + <stat/> |
1783 | + <detectors> |
1784 | + <exclude_from_detectors/> |
1785 | + </detectors> |
1786 | + </prescribed> |
1787 | + </scalar_field> |
1788 | + <scalar_field name="Viscosity" rank="0"> |
1789 | + <prescribed> |
1790 | + <mesh name="DGP0"/> |
1791 | + <value name="WholeMesh"> |
1792 | + <constant> |
1793 | + <real_value rank="0">1.0e-04</real_value> |
1794 | + </constant> |
1795 | + </value> |
1796 | + <output/> |
1797 | + <stat/> |
1798 | + <detectors> |
1799 | + <exclude_from_detectors/> |
1800 | + </detectors> |
1801 | + </prescribed> |
1802 | + </scalar_field> |
1803 | + <scalar_field name="ControlVolumeDivergence" rank="0"> |
1804 | + <diagnostic field_name="Velocity"> |
1805 | + <algorithm name="Internal" material_phase_support="multiple"/> |
1806 | + <mesh name="PressureMesh"/> |
1807 | + <output/> |
1808 | + <stat/> |
1809 | + <convergence> |
1810 | + <include_in_convergence/> |
1811 | + </convergence> |
1812 | + <detectors> |
1813 | + <include_in_detectors/> |
1814 | + </detectors> |
1815 | + <steady_state> |
1816 | + <include_in_steady_state/> |
1817 | + </steady_state> |
1818 | + </diagnostic> |
1819 | + </scalar_field> |
1820 | + <vector_field name="interstitial_velocity" rank="1"> |
1821 | + <diagnostic> |
1822 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
1823 | + <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"] |
1824 | +darcy_vel = states["fluid"].vector_fields["Velocity"] |
1825 | + |
1826 | +for ele in range(field.element_count): |
1827 | + phi_ele = phi.node_val(ele) |
1828 | + |
1829 | + factor_ele = 1.0/max(phi_ele,1.0e-08) |
1830 | + |
1831 | + vel_ele = [] |
1832 | + for i in range(len(darcy_vel.ele_val(ele))): |
1833 | + vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i]) |
1834 | + |
1835 | + field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value> |
1836 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
1837 | + |
1838 | +also assume that the velocity is dg</comment> |
1839 | + </algorithm> |
1840 | + <mesh name="VelocityMesh"/> |
1841 | + <output/> |
1842 | + <stat> |
1843 | + <include_in_stat/> |
1844 | + </stat> |
1845 | + <convergence> |
1846 | + <include_in_convergence/> |
1847 | + </convergence> |
1848 | + <detectors> |
1849 | + <include_in_detectors/> |
1850 | + </detectors> |
1851 | + <steady_state> |
1852 | + <include_in_steady_state/> |
1853 | + </steady_state> |
1854 | + </diagnostic> |
1855 | + </vector_field> |
1856 | + </material_phase> |
1857 | +</fluidity_options> |
1858 | |
1859 | === added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet/square.geo' |
1860 | --- tests/darcy_p0p1_test_cty_cv_pressBCinlet/square.geo 1970-01-01 00:00:00 +0000 |
1861 | +++ tests/darcy_p0p1_test_cty_cv_pressBCinlet/square.geo 2011-12-23 12:05:25 +0000 |
1862 | @@ -0,0 +1,22 @@ |
1863 | +// define a layer variable |
1864 | +lay = 5; |
1865 | + |
1866 | +// define a len variable |
1867 | +len = 300.0; |
1868 | + |
1869 | +Point(1) = {0.0, 0.0, 0.0, 1.0}; |
1870 | + |
1871 | +Extrude {len, 0.0, 0.0} { |
1872 | + Point{1}; Layers{lay}; |
1873 | +} |
1874 | + |
1875 | +Extrude {0.0, len, 0.0} { |
1876 | + Line{1}; Layers{lay}; |
1877 | +} |
1878 | + |
1879 | +Physical Line(7) = {1}; |
1880 | +Physical Line(8) = {4}; |
1881 | +Physical Line(9) = {2}; |
1882 | +Physical Line(10) = {3}; |
1883 | + |
1884 | +Physical Surface(11) = {5}; |
1885 | |
1886 | === added directory 'tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain' |
1887 | === added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/Makefile' |
1888 | --- tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/Makefile 1970-01-01 00:00:00 +0000 |
1889 | +++ tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/Makefile 2011-12-23 12:05:25 +0000 |
1890 | @@ -0,0 +1,23 @@ |
1891 | +SHELL = sh |
1892 | + |
1893 | +input: clean |
1894 | + gmsh -2 square_rotated.geo -o square_rotated.msh |
1895 | + |
1896 | +clean: |
1897 | + rm -f *.ele |
1898 | + rm -f *.node |
1899 | + rm -f *.bound |
1900 | + rm -f *.edge |
1901 | + rm -f *.face |
1902 | + rm -f *.vtu |
1903 | + rm -f *.pvtu |
1904 | + rm -f *.s |
1905 | + rm -f *.stat |
1906 | + rm -f *.log-0 |
1907 | + rm -f *.err-0 |
1908 | + rm -f *.msh |
1909 | + rm -f *.halo |
1910 | + rm -f fluidity.err* |
1911 | + rm -f fluidity.log* |
1912 | + rm -f matrixdump* |
1913 | + rm -f first_timestep_adapted_mesh* |
1914 | |
1915 | === added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain.xml' |
1916 | --- tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain.xml 1970-01-01 00:00:00 +0000 |
1917 | +++ tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain.xml 2011-12-23 12:05:25 +0000 |
1918 | @@ -0,0 +1,46 @@ |
1919 | +<?xml version="1.0" encoding="UTF-8" ?> |
1920 | +<!DOCTYPE testproblem SYSTEM "regressiontest.dtd"> |
1921 | + |
1922 | +<testproblem> |
1923 | + <name>darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain</name> |
1924 | + <owner userid="btollit"/> |
1925 | + <tags>flml darcy</tags> |
1926 | + <problem_definition length="short" nprocs="1"> |
1927 | + <command_line> |
1928 | +../../bin/fluidity darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d.flml |
1929 | + </command_line> |
1930 | + <!-- Two dimensional problem (actually 1d done in 2d) for darcy flow with one region and one material using p0p1_test_cty_cv element type testing the pressure gradient against analytic and the interstitial velocity. This tests using a diagnostic velocity absorption field associated with a mesh and momentum dg for darcy flow. This will use a rotated gmsh domain to test applying the weak pressure BC at a non cartesian aligned axis.--> |
1931 | + </problem_definition> |
1932 | + <variables> |
1933 | + <variable name="pressure_2d" language="python"> |
1934 | +import vtktools |
1935 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d_1.vtu") |
1936 | +pressure_2d = v.GetScalarRange("Pressure") |
1937 | + </variable> |
1938 | + <variable name="inter_vel_2d" language="python"> |
1939 | +import vtktools |
1940 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d_1.vtu") |
1941 | +inter_vel_2d = max(v.GetScalarRange("interstitial_velocity")) |
1942 | + </variable> |
1943 | + </variables> |
1944 | + <pass_tests> |
1945 | + <test name="change_P for 2d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python"> |
1946 | +change_P = abs(max(pressure_2d) - min(pressure_2d)) |
1947 | +visc = 1.0e-04 |
1948 | +darcy_vel_BC = 5.0 |
1949 | +perm = 1.0e-10 |
1950 | +domain_length = 300.0 |
1951 | +print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) |
1952 | +assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) < 1.0e-09 |
1953 | + </test> |
1954 | + <test name="interstitial velocity 2d should equal darcy_vel/porosity, check relative difference to be under tolerance 1.0e-09" language="python"> |
1955 | +import math |
1956 | +# the analytic velocity magnitude is 10.0, with an incline of 45 degrees (= pi/4.0 radians) |
1957 | +analytic_vel = 10.0 * math.sin(math.pi/4.0) |
1958 | +print 'Relative error of interstitial velocity: ',abs((inter_vel_2d - analytic_vel)/analytic_vel) |
1959 | +assert abs((inter_vel_2d - analytic_vel)/analytic_vel) < 1.0e-09 |
1960 | + </test> |
1961 | + </pass_tests> |
1962 | + <warn_tests> |
1963 | + </warn_tests> |
1964 | +</testproblem> |
1965 | |
1966 | === added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d.flml' |
1967 | --- tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d.flml 1970-01-01 00:00:00 +0000 |
1968 | +++ tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d.flml 2011-12-23 12:05:25 +0000 |
1969 | @@ -0,0 +1,372 @@ |
1970 | +<?xml version='1.0' encoding='utf-8'?> |
1971 | +<fluidity_options> |
1972 | + <simulation_name> |
1973 | + <string_value lines="1">darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain_2d</string_value> |
1974 | + </simulation_name> |
1975 | + <problem_type> |
1976 | + <string_value lines="1">fluids</string_value> |
1977 | + </problem_type> |
1978 | + <geometry> |
1979 | + <dimension> |
1980 | + <integer_value rank="0">2</integer_value> |
1981 | + </dimension> |
1982 | + <mesh name="CoordinateMesh"> |
1983 | + <from_file file_name="square_rotated"> |
1984 | + <format name="gmsh"/> |
1985 | + <stat> |
1986 | + <include_in_stat/> |
1987 | + </stat> |
1988 | + </from_file> |
1989 | + </mesh> |
1990 | + <mesh name="VelocityMesh"> |
1991 | + <from_mesh> |
1992 | + <mesh name="CoordinateMesh"/> |
1993 | + <mesh_shape> |
1994 | + <polynomial_degree> |
1995 | + <integer_value rank="0">0</integer_value> |
1996 | + </polynomial_degree> |
1997 | + </mesh_shape> |
1998 | + <mesh_continuity> |
1999 | + <string_value>discontinuous</string_value> |
2000 | + </mesh_continuity> |
2001 | + <stat> |
2002 | + <exclude_from_stat/> |
2003 | + </stat> |
2004 | + </from_mesh> |
2005 | + </mesh> |
2006 | + <mesh name="PressureMesh"> |
2007 | + <from_mesh> |
2008 | + <mesh name="CoordinateMesh"/> |
2009 | + <mesh_shape> |
2010 | + <polynomial_degree> |
2011 | + <integer_value rank="0">1</integer_value> |
2012 | + </polynomial_degree> |
2013 | + </mesh_shape> |
2014 | + <stat> |
2015 | + <exclude_from_stat/> |
2016 | + </stat> |
2017 | + </from_mesh> |
2018 | + </mesh> |
2019 | + <mesh name="DGP0"> |
2020 | + <from_mesh> |
2021 | + <mesh name="CoordinateMesh"/> |
2022 | + <mesh_shape> |
2023 | + <polynomial_degree> |
2024 | + <integer_value rank="0">0</integer_value> |
2025 | + </polynomial_degree> |
2026 | + </mesh_shape> |
2027 | + <mesh_continuity> |
2028 | + <string_value>discontinuous</string_value> |
2029 | + </mesh_continuity> |
2030 | + <stat> |
2031 | + <exclude_from_stat/> |
2032 | + </stat> |
2033 | + </from_mesh> |
2034 | + </mesh> |
2035 | + <mesh name="OutputMesh"> |
2036 | + <from_mesh> |
2037 | + <mesh name="CoordinateMesh"/> |
2038 | + <mesh_continuity> |
2039 | + <string_value>discontinuous</string_value> |
2040 | + </mesh_continuity> |
2041 | + <stat> |
2042 | + <exclude_from_stat/> |
2043 | + </stat> |
2044 | + </from_mesh> |
2045 | + </mesh> |
2046 | + <quadrature> |
2047 | + <degree> |
2048 | + <integer_value rank="0">5</integer_value> |
2049 | + </degree> |
2050 | + </quadrature> |
2051 | + </geometry> |
2052 | + <io> |
2053 | + <dump_format> |
2054 | + <string_value>vtk</string_value> |
2055 | + </dump_format> |
2056 | + <dump_period_in_timesteps> |
2057 | + <constant> |
2058 | + <integer_value rank="0">0</integer_value> |
2059 | + </constant> |
2060 | + </dump_period_in_timesteps> |
2061 | + <output_mesh name="OutputMesh"/> |
2062 | + <stat/> |
2063 | + </io> |
2064 | + <timestepping> |
2065 | + <current_time> |
2066 | + <real_value rank="0">0.0</real_value> |
2067 | + </current_time> |
2068 | + <timestep> |
2069 | + <real_value rank="0">1.0</real_value> |
2070 | + </timestep> |
2071 | + <finish_time> |
2072 | + <real_value rank="0">1.0</real_value> |
2073 | + </finish_time> |
2074 | + </timestepping> |
2075 | + <material_phase name="fluid"> |
2076 | + <scalar_field name="Pressure" rank="0"> |
2077 | + <prognostic> |
2078 | + <mesh name="PressureMesh"/> |
2079 | + <spatial_discretisation> |
2080 | + <continuous_galerkin> |
2081 | + <test_continuity_with_cv_dual/> |
2082 | + </continuous_galerkin> |
2083 | + </spatial_discretisation> |
2084 | + <scheme> |
2085 | + <poisson_pressure_solution> |
2086 | + <string_value lines="1">never</string_value> |
2087 | + </poisson_pressure_solution> |
2088 | + <use_projection_method/> |
2089 | + </scheme> |
2090 | + <solver> |
2091 | + <iterative_method name="cg"/> |
2092 | + <preconditioner name="sor"/> |
2093 | + <relative_error> |
2094 | + <real_value rank="0">1.0e-10</real_value> |
2095 | + </relative_error> |
2096 | + <max_iterations> |
2097 | + <integer_value rank="0">1000</integer_value> |
2098 | + </max_iterations> |
2099 | + <never_ignore_solver_failures/> |
2100 | + <diagnostics> |
2101 | + <monitors/> |
2102 | + </diagnostics> |
2103 | + </solver> |
2104 | + <boundary_conditions name="right_outflow"> |
2105 | + <surface_ids> |
2106 | + <integer_value shape="1" rank="1">8</integer_value> |
2107 | + </surface_ids> |
2108 | + <type name="dirichlet"> |
2109 | + <constant> |
2110 | + <real_value rank="0">0.0</real_value> |
2111 | + </constant> |
2112 | + </type> |
2113 | + </boundary_conditions> |
2114 | + <boundary_conditions name="left_inflow"> |
2115 | + <surface_ids> |
2116 | + <integer_value shape="1" rank="1">10</integer_value> |
2117 | + </surface_ids> |
2118 | + <type name="dirichlet"> |
2119 | + <apply_weakly/> |
2120 | + <constant> |
2121 | + <real_value rank="0">1.5e+09</real_value> |
2122 | + </constant> |
2123 | + </type> |
2124 | + </boundary_conditions> |
2125 | + <output/> |
2126 | + <stat/> |
2127 | + <convergence> |
2128 | + <include_in_convergence/> |
2129 | + </convergence> |
2130 | + <detectors> |
2131 | + <exclude_from_detectors/> |
2132 | + </detectors> |
2133 | + <steady_state> |
2134 | + <include_in_steady_state/> |
2135 | + </steady_state> |
2136 | + <no_interpolation/> |
2137 | + </prognostic> |
2138 | + </scalar_field> |
2139 | + <vector_field name="Velocity" rank="1"> |
2140 | + <prognostic> |
2141 | + <mesh name="VelocityMesh"/> |
2142 | + <equation name="Boussinesq"/> |
2143 | + <spatial_discretisation> |
2144 | + <discontinuous_galerkin> |
2145 | + <mass_terms> |
2146 | + <exclude_mass_terms/> |
2147 | + </mass_terms> |
2148 | + <viscosity_scheme> |
2149 | + <compact_discontinuous_galerkin/> |
2150 | + </viscosity_scheme> |
2151 | + <advection_scheme> |
2152 | + <none/> |
2153 | + <integrate_advection_by_parts> |
2154 | + <twice/> |
2155 | + </integrate_advection_by_parts> |
2156 | + </advection_scheme> |
2157 | + </discontinuous_galerkin> |
2158 | + <conservative_advection> |
2159 | + <real_value rank="0">0.0</real_value> |
2160 | + </conservative_advection> |
2161 | + </spatial_discretisation> |
2162 | + <temporal_discretisation> |
2163 | + <theta> |
2164 | + <real_value rank="0">1.0</real_value> |
2165 | + </theta> |
2166 | + <relaxation> |
2167 | + <real_value rank="0">1.0</real_value> |
2168 | + </relaxation> |
2169 | + </temporal_discretisation> |
2170 | + <solver> |
2171 | + <iterative_method name="gmres"> |
2172 | + <restart> |
2173 | + <integer_value rank="0">30</integer_value> |
2174 | + </restart> |
2175 | + </iterative_method> |
2176 | + <preconditioner name="sor"/> |
2177 | + <relative_error> |
2178 | + <real_value rank="0">1.0e-10</real_value> |
2179 | + </relative_error> |
2180 | + <max_iterations> |
2181 | + <integer_value rank="0">1000</integer_value> |
2182 | + </max_iterations> |
2183 | + <never_ignore_solver_failures/> |
2184 | + <diagnostics> |
2185 | + <monitors/> |
2186 | + </diagnostics> |
2187 | + </solver> |
2188 | + <initial_condition name="WholeMesh"> |
2189 | + <constant> |
2190 | + <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value> |
2191 | + </constant> |
2192 | + </initial_condition> |
2193 | + <boundary_conditions name="no_outflow_top_bottom"> |
2194 | + <surface_ids> |
2195 | + <integer_value shape="2" rank="1">7 9</integer_value> |
2196 | + </surface_ids> |
2197 | + <type name="no_normal_flow"/> |
2198 | + </boundary_conditions> |
2199 | + <vector_field name="Absorption" rank="1"> |
2200 | + <diagnostic> |
2201 | + <mesh name="DGP0"/> |
2202 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
2203 | + <string_value lines="20" type="python"># get the prescribed fields |
2204 | +perm = states["fluid"].scalar_fields["Permeability"] |
2205 | +visc = states["fluid"].scalar_fields["Viscosity"] |
2206 | + |
2207 | +# loop the DMmesh_p0 nodes (which are element centred) |
2208 | +for n in range(field.node_count): |
2209 | + perm_n = perm.node_val(n) |
2210 | + visc_n = visc.node_val(n) |
2211 | + |
2212 | + # calc the absorption term |
2213 | + sigma_n = visc_n/perm_n |
2214 | + |
2215 | + # set the absorption term |
2216 | + field.set(n,sigma_n)</string_value> |
2217 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
2218 | + |
2219 | +also assume that the velocity is dg</comment> |
2220 | + </algorithm> |
2221 | + <output/> |
2222 | + <stat> |
2223 | + <include_in_stat/> |
2224 | + </stat> |
2225 | + <convergence> |
2226 | + <include_in_convergence/> |
2227 | + </convergence> |
2228 | + <detectors> |
2229 | + <include_in_detectors/> |
2230 | + </detectors> |
2231 | + <steady_state> |
2232 | + <include_in_steady_state/> |
2233 | + </steady_state> |
2234 | + </diagnostic> |
2235 | + <include_pressure_correction/> |
2236 | + </vector_field> |
2237 | + <output/> |
2238 | + <stat> |
2239 | + <include_in_stat/> |
2240 | + <previous_time_step> |
2241 | + <exclude_from_stat/> |
2242 | + </previous_time_step> |
2243 | + <nonlinear_field> |
2244 | + <exclude_from_stat/> |
2245 | + </nonlinear_field> |
2246 | + </stat> |
2247 | + <convergence> |
2248 | + <include_in_convergence/> |
2249 | + </convergence> |
2250 | + <detectors> |
2251 | + <include_in_detectors/> |
2252 | + </detectors> |
2253 | + <steady_state> |
2254 | + <include_in_steady_state/> |
2255 | + </steady_state> |
2256 | + <consistent_interpolation/> |
2257 | + </prognostic> |
2258 | + </vector_field> |
2259 | + <scalar_field name="Porosity" rank="0"> |
2260 | + <prescribed> |
2261 | + <mesh name="DGP0"/> |
2262 | + <value name="WholeMesh"> |
2263 | + <constant> |
2264 | + <real_value rank="0">0.5</real_value> |
2265 | + </constant> |
2266 | + </value> |
2267 | + <output/> |
2268 | + <stat/> |
2269 | + <detectors> |
2270 | + <exclude_from_detectors/> |
2271 | + </detectors> |
2272 | + </prescribed> |
2273 | + </scalar_field> |
2274 | + <scalar_field name="Permeability" rank="0"> |
2275 | + <prescribed> |
2276 | + <mesh name="DGP0"/> |
2277 | + <value name="WholeMesh"> |
2278 | + <constant> |
2279 | + <real_value rank="0">1.0e-10</real_value> |
2280 | + </constant> |
2281 | + </value> |
2282 | + <output/> |
2283 | + <stat/> |
2284 | + <detectors> |
2285 | + <exclude_from_detectors/> |
2286 | + </detectors> |
2287 | + </prescribed> |
2288 | + </scalar_field> |
2289 | + <scalar_field name="Viscosity" rank="0"> |
2290 | + <prescribed> |
2291 | + <mesh name="DGP0"/> |
2292 | + <value name="WholeMesh"> |
2293 | + <constant> |
2294 | + <real_value rank="0">1.0e-04</real_value> |
2295 | + </constant> |
2296 | + </value> |
2297 | + <output/> |
2298 | + <stat/> |
2299 | + <detectors> |
2300 | + <exclude_from_detectors/> |
2301 | + </detectors> |
2302 | + </prescribed> |
2303 | + </scalar_field> |
2304 | + <vector_field name="interstitial_velocity" rank="1"> |
2305 | + <diagnostic> |
2306 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
2307 | + <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"] |
2308 | +darcy_vel = states["fluid"].vector_fields["Velocity"] |
2309 | + |
2310 | +for ele in range(field.element_count): |
2311 | + phi_ele = phi.node_val(ele) |
2312 | + |
2313 | + factor_ele = 1.0/max(phi_ele,1.0e-08) |
2314 | + |
2315 | + vel_ele = [] |
2316 | + for i in range(len(darcy_vel.ele_val(ele))): |
2317 | + vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i]) |
2318 | + |
2319 | + field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value> |
2320 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
2321 | + |
2322 | +also assume that the velocity is dg</comment> |
2323 | + </algorithm> |
2324 | + <mesh name="VelocityMesh"/> |
2325 | + <output/> |
2326 | + <stat> |
2327 | + <include_in_stat/> |
2328 | + </stat> |
2329 | + <convergence> |
2330 | + <include_in_convergence/> |
2331 | + </convergence> |
2332 | + <detectors> |
2333 | + <include_in_detectors/> |
2334 | + </detectors> |
2335 | + <steady_state> |
2336 | + <include_in_steady_state/> |
2337 | + </steady_state> |
2338 | + </diagnostic> |
2339 | + </vector_field> |
2340 | + </material_phase> |
2341 | +</fluidity_options> |
2342 | |
2343 | === added file 'tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/square_rotated.geo' |
2344 | --- tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/square_rotated.geo 1970-01-01 00:00:00 +0000 |
2345 | +++ tests/darcy_p0p1_test_cty_cv_pressBCinlet_rotated_domain/square_rotated.geo 2011-12-23 12:05:25 +0000 |
2346 | @@ -0,0 +1,38 @@ |
2347 | +// define a characteristic length |
2348 | +lc = 60.0; |
2349 | + |
2350 | +len = 300.0; |
2351 | + |
2352 | +// define points |
2353 | +Point(1) = {0.0,0.0,0.0,lc}; |
2354 | +Point(2) = {len,0.0,0.0,lc}; |
2355 | +Point(3) = {len,len,0.0,lc}; |
2356 | +Point(4) = {0.0,len,0.0,lc}; |
2357 | + |
2358 | +// define line's |
2359 | +Line(1) = {1,2}; |
2360 | +Line(2) = {2,3}; |
2361 | +Line(3) = {3,4}; |
2362 | +Line(4) = {4,1}; |
2363 | + |
2364 | +// define line loop |
2365 | +Line Loop(1) = {1,2,3,4}; |
2366 | + |
2367 | +// define surface |
2368 | +Plane Surface(1) = {1}; |
2369 | + |
2370 | +// Rotate the domain about the z axis by 45 degrees (Pi/4) |
2371 | +Rotate {{0, 0, 1}, {0, 0, 0}, Pi/4} { |
2372 | + Surface{1}; |
2373 | +} |
2374 | + |
2375 | +// Tag the lines and surface |
2376 | +Physical Line(7) = {1}; |
2377 | +Physical Line(8) = {2}; |
2378 | +Physical Line(9) = {3}; |
2379 | +Physical Line(10) = {4}; |
2380 | + |
2381 | +Physical Surface(11) = {1}; |
2382 | + |
2383 | +// Make the grid regular |
2384 | +Transfinite Surface"*"; |
2385 | |
2386 | === added directory 'tests/darcy_p0p1_test_cty_cv_velBCinlet' |
2387 | === added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/Makefile' |
2388 | --- tests/darcy_p0p1_test_cty_cv_velBCinlet/Makefile 1970-01-01 00:00:00 +0000 |
2389 | +++ tests/darcy_p0p1_test_cty_cv_velBCinlet/Makefile 2011-12-23 12:05:25 +0000 |
2390 | @@ -0,0 +1,32 @@ |
2391 | +SHELL = sh |
2392 | + |
2393 | +SIM = darcy_p0p1_test_cty_cv_velBCinlet |
2394 | + |
2395 | +# options for 1d mesh interval script |
2396 | +MESH_SIZE = 100.0 |
2397 | +LEFT_X = 0.0 |
2398 | +RIGHT_X = 300.0 |
2399 | + |
2400 | +input: clean |
2401 | + interval --dx=$(MESH_SIZE) $(LEFT_X) $(RIGHT_X) $(SIM)_1d |
2402 | + gmsh -2 square.geo -o square.msh |
2403 | + gmsh -3 cube.geo -o cube.msh |
2404 | + |
2405 | +clean: |
2406 | + rm -f *.ele |
2407 | + rm -f *.node |
2408 | + rm -f *.bound |
2409 | + rm -f *.edge |
2410 | + rm -f *.face |
2411 | + rm -f *.vtu |
2412 | + rm -f *.pvtu |
2413 | + rm -f *.s |
2414 | + rm -f *.stat |
2415 | + rm -f *.log-0 |
2416 | + rm -f *.err-0 |
2417 | + rm -f *.msh |
2418 | + rm -f *.halo |
2419 | + rm -f fluidity.err* |
2420 | + rm -f fluidity.log* |
2421 | + rm -f matrixdump* |
2422 | + rm -f first_timestep_adapted_mesh* |
2423 | |
2424 | === added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/cube.geo' |
2425 | --- tests/darcy_p0p1_test_cty_cv_velBCinlet/cube.geo 1970-01-01 00:00:00 +0000 |
2426 | +++ tests/darcy_p0p1_test_cty_cv_velBCinlet/cube.geo 2011-12-23 12:05:25 +0000 |
2427 | @@ -0,0 +1,39 @@ |
2428 | +// define a layer variable |
2429 | +lay = 5; |
2430 | + |
2431 | +// define a len variable |
2432 | +len = 300.0; |
2433 | + |
2434 | +Point(1) = {0.0, 0.0, 0.0, 1.0}; |
2435 | + |
2436 | +Extrude {len, 0.0, 0.0} { |
2437 | + Point{1}; Layers{lay}; |
2438 | +} |
2439 | + |
2440 | +Extrude {0.0, len, 0.0} { |
2441 | + Line{1}; Layers{lay}; |
2442 | +} |
2443 | + |
2444 | +Extrude {0.0, 0.0, len} { |
2445 | + Surface{5}; Layers{lay}; |
2446 | +} |
2447 | + |
2448 | +// bottom |
2449 | +Physical Surface(1) = {5}; |
2450 | + |
2451 | +// top |
2452 | +Physical Surface(2) = {27}; |
2453 | + |
2454 | +// left |
2455 | +Physical Surface(3) = {26}; |
2456 | + |
2457 | +// right |
2458 | +Physical Surface(4) = {18}; |
2459 | + |
2460 | +// front |
2461 | +Physical Surface(5) = {14}; |
2462 | + |
2463 | +// back |
2464 | +Physical Surface(6) = {22}; |
2465 | + |
2466 | +Physical Volume(1) = {1}; |
2467 | |
2468 | === added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet.xml' |
2469 | --- tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet.xml 1970-01-01 00:00:00 +0000 |
2470 | +++ tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet.xml 2011-12-23 12:05:25 +0000 |
2471 | @@ -0,0 +1,94 @@ |
2472 | +<?xml version="1.0" encoding="UTF-8" ?> |
2473 | +<!DOCTYPE testproblem SYSTEM "regressiontest.dtd"> |
2474 | + |
2475 | +<testproblem> |
2476 | + <name>darcy_p0p1_test_cty_cv_velBCinlet</name> |
2477 | + <owner userid="btollit"/> |
2478 | + <tags>flml darcy</tags> |
2479 | + <problem_definition length="short" nprocs="1"> |
2480 | + <command_line> |
2481 | +../../bin/fluidity darcy_p0p1_test_cty_cv_velBCinlet_1d.flml |
2482 | +../../bin/fluidity darcy_p0p1_test_cty_cv_velBCinlet_2d.flml |
2483 | +../../bin/fluidity darcy_p0p1_test_cty_cv_velBCinlet_3d.flml |
2484 | + </command_line> |
2485 | + <!-- One/two/three dimensional problem for darcy flow with one region and one material using p0p1_test_cty_cv element type testing the pressure gradient against analytic and the interstitial velocity. This tests using a diagnostic velocity absorption field associated with a mesh and momentum dg for darcy flow--> |
2486 | + </problem_definition> |
2487 | + <variables> |
2488 | + <variable name="pressure_1d" language="python"> |
2489 | +import vtktools |
2490 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_velBCinlet_1d_1.vtu") |
2491 | +pressure_1d = v.GetScalarRange("Pressure") |
2492 | + </variable> |
2493 | + <variable name="inter_vel_1d" language="python"> |
2494 | +import vtktools |
2495 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_velBCinlet_1d_1.vtu") |
2496 | +inter_vel_1d = max(v.GetScalarRange("interstitial_velocity")) |
2497 | + </variable> |
2498 | + <variable name="pressure_2d" language="python"> |
2499 | +import vtktools |
2500 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_velBCinlet_2d_1.vtu") |
2501 | +pressure_2d = v.GetScalarRange("Pressure") |
2502 | + </variable> |
2503 | + <variable name="inter_vel_2d" language="python"> |
2504 | +import vtktools |
2505 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_velBCinlet_2d_1.vtu") |
2506 | +inter_vel_2d = max(v.GetScalarRange("interstitial_velocity")) |
2507 | + </variable> |
2508 | + <variable name="pressure_3d" language="python"> |
2509 | +import vtktools |
2510 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_velBCinlet_3d_1.vtu") |
2511 | +pressure_3d = v.GetScalarRange("Pressure") |
2512 | + </variable> |
2513 | + <variable name="inter_vel_3d" language="python"> |
2514 | +import vtktools |
2515 | +v = vtktools.vtu("darcy_p0p1_test_cty_cv_velBCinlet_3d_1.vtu") |
2516 | +inter_vel_3d = max(v.GetScalarRange("interstitial_velocity")) |
2517 | + </variable> |
2518 | + </variables> |
2519 | + <pass_tests> |
2520 | + <test name="change_P for 1d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python"> |
2521 | +change_P = abs(max(pressure_1d) - min(pressure_1d)) |
2522 | +visc = 1.0e-04 |
2523 | +darcy_vel_BC = 5.0 |
2524 | +perm = 1.0e-10 |
2525 | +domain_length = 300.0 |
2526 | +print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) |
2527 | +assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) < 1.0e-09 |
2528 | + </test> |
2529 | + <test name="interstitial velocity 1d should equal darcy_vel/porosity), check relative difference to be under tolerance 1.0e-09" language="python"> |
2530 | +analytic_vel = 10.0 |
2531 | +print 'Relative error of interstitial velocity: ',abs((inter_vel_1d - analytic_vel)/analytic_vel) |
2532 | +assert abs((inter_vel_1d - analytic_vel)/analytic_vel) < 1.0e-09 |
2533 | + </test> |
2534 | + <test name="change_P for 2d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python"> |
2535 | +change_P = abs(max(pressure_2d) - min(pressure_2d)) |
2536 | +visc = 1.0e-04 |
2537 | +darcy_vel_BC = 5.0 |
2538 | +perm = 1.0e-10 |
2539 | +domain_length = 300.0 |
2540 | +print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) |
2541 | +assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) < 1.0e-09 |
2542 | + </test> |
2543 | + <test name="interstitial velocity 2d should equal darcy_vel/porosity), check relative difference to be under tolerance 1.0e-09" language="python"> |
2544 | +analytic_vel = 10.0 |
2545 | +print 'Relative error of interstitial velocity: ',abs((inter_vel_2d - analytic_vel)/analytic_vel) |
2546 | +assert abs((inter_vel_2d - analytic_vel)/analytic_vel) < 1.0e-09 |
2547 | + </test> |
2548 | + <test name="change_P for 3d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python"> |
2549 | +change_P = abs(max(pressure_3d) - min(pressure_3d)) |
2550 | +visc = 1.0e-04 |
2551 | +darcy_vel_BC = 5.0 |
2552 | +perm = 1.0e-10 |
2553 | +domain_length = 300.0 |
2554 | +print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) |
2555 | +assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) < 1.0e-09 |
2556 | + </test> |
2557 | + <test name="interstitial velocity 3d should equal darcy_vel/porosity), check relative difference to be under tolerance 1.0e-09" language="python"> |
2558 | +analytic_vel = 10.0 |
2559 | +print 'Relative error of interstitial velocity: ',abs((inter_vel_3d - analytic_vel)/analytic_vel) |
2560 | +assert abs((inter_vel_3d - analytic_vel)/analytic_vel) < 1.0e-09 |
2561 | + </test> |
2562 | + </pass_tests> |
2563 | + <warn_tests> |
2564 | + </warn_tests> |
2565 | +</testproblem> |
2566 | |
2567 | === added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_1d.flml' |
2568 | --- tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_1d.flml 1970-01-01 00:00:00 +0000 |
2569 | +++ tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_1d.flml 2011-12-23 12:05:25 +0000 |
2570 | @@ -0,0 +1,391 @@ |
2571 | +<?xml version='1.0' encoding='utf-8'?> |
2572 | +<fluidity_options> |
2573 | + <simulation_name> |
2574 | + <string_value lines="1">darcy_p0p1_test_cty_cv_velBCinlet_1d</string_value> |
2575 | + <comment>a simple short test case for darcy flow using the element type p0p1_test_cty_cv for velocity-pressure |
2576 | + |
2577 | +the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution) |
2578 | + |
2579 | +it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct |
2580 | + |
2581 | +the darcy flow equation is |
2582 | + |
2583 | +sigma*darcy_vel = - grad P |
2584 | + |
2585 | +sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options) |
2586 | + |
2587 | +python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity |
2588 | + |
2589 | +this case has a 1 region 1 material with velocity boundary condition inflow. |
2590 | + |
2591 | +to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input. |
2592 | + |
2593 | +the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) |
2594 | + |
2595 | +the geometry is 1d and one time step is done (as nothing depends on time)</comment> |
2596 | + </simulation_name> |
2597 | + <problem_type> |
2598 | + <string_value lines="1">fluids</string_value> |
2599 | + </problem_type> |
2600 | + <geometry> |
2601 | + <dimension> |
2602 | + <integer_value rank="0">1</integer_value> |
2603 | + </dimension> |
2604 | + <mesh name="CoordinateMesh"> |
2605 | + <from_file file_name="darcy_p0p1_test_cty_cv_velBCinlet_1d"> |
2606 | + <format name="triangle"/> |
2607 | + <stat> |
2608 | + <include_in_stat/> |
2609 | + </stat> |
2610 | + </from_file> |
2611 | + </mesh> |
2612 | + <mesh name="VelocityMesh"> |
2613 | + <from_mesh> |
2614 | + <mesh name="CoordinateMesh"/> |
2615 | + <mesh_shape> |
2616 | + <polynomial_degree> |
2617 | + <integer_value rank="0">0</integer_value> |
2618 | + </polynomial_degree> |
2619 | + </mesh_shape> |
2620 | + <mesh_continuity> |
2621 | + <string_value>discontinuous</string_value> |
2622 | + </mesh_continuity> |
2623 | + <stat> |
2624 | + <exclude_from_stat/> |
2625 | + </stat> |
2626 | + </from_mesh> |
2627 | + </mesh> |
2628 | + <mesh name="PressureMesh"> |
2629 | + <from_mesh> |
2630 | + <mesh name="CoordinateMesh"/> |
2631 | + <mesh_shape> |
2632 | + <polynomial_degree> |
2633 | + <integer_value rank="0">1</integer_value> |
2634 | + </polynomial_degree> |
2635 | + </mesh_shape> |
2636 | + <stat> |
2637 | + <exclude_from_stat/> |
2638 | + </stat> |
2639 | + </from_mesh> |
2640 | + </mesh> |
2641 | + <mesh name="DGP0"> |
2642 | + <from_mesh> |
2643 | + <mesh name="CoordinateMesh"/> |
2644 | + <mesh_shape> |
2645 | + <polynomial_degree> |
2646 | + <integer_value rank="0">0</integer_value> |
2647 | + </polynomial_degree> |
2648 | + </mesh_shape> |
2649 | + <mesh_continuity> |
2650 | + <string_value>discontinuous</string_value> |
2651 | + </mesh_continuity> |
2652 | + <stat> |
2653 | + <exclude_from_stat/> |
2654 | + </stat> |
2655 | + </from_mesh> |
2656 | + </mesh> |
2657 | + <mesh name="OutputMesh"> |
2658 | + <from_mesh> |
2659 | + <mesh name="CoordinateMesh"/> |
2660 | + <mesh_continuity> |
2661 | + <string_value>discontinuous</string_value> |
2662 | + </mesh_continuity> |
2663 | + <stat> |
2664 | + <exclude_from_stat/> |
2665 | + </stat> |
2666 | + </from_mesh> |
2667 | + </mesh> |
2668 | + <quadrature> |
2669 | + <degree> |
2670 | + <integer_value rank="0">2</integer_value> |
2671 | + </degree> |
2672 | + </quadrature> |
2673 | + </geometry> |
2674 | + <io> |
2675 | + <dump_format> |
2676 | + <string_value>vtk</string_value> |
2677 | + </dump_format> |
2678 | + <dump_period_in_timesteps> |
2679 | + <constant> |
2680 | + <integer_value rank="0">0</integer_value> |
2681 | + </constant> |
2682 | + </dump_period_in_timesteps> |
2683 | + <output_mesh name="OutputMesh"/> |
2684 | + <stat/> |
2685 | + </io> |
2686 | + <timestepping> |
2687 | + <current_time> |
2688 | + <real_value rank="0">0.0</real_value> |
2689 | + </current_time> |
2690 | + <timestep> |
2691 | + <real_value rank="0">1.0</real_value> |
2692 | + </timestep> |
2693 | + <finish_time> |
2694 | + <real_value rank="0">1.0</real_value> |
2695 | + </finish_time> |
2696 | + </timestepping> |
2697 | + <material_phase name="fluid"> |
2698 | + <scalar_field name="Pressure" rank="0"> |
2699 | + <prognostic> |
2700 | + <mesh name="PressureMesh"/> |
2701 | + <spatial_discretisation> |
2702 | + <continuous_galerkin> |
2703 | + <test_continuity_with_cv_dual/> |
2704 | + </continuous_galerkin> |
2705 | + </spatial_discretisation> |
2706 | + <scheme> |
2707 | + <poisson_pressure_solution> |
2708 | + <string_value lines="1">never</string_value> |
2709 | + </poisson_pressure_solution> |
2710 | + <use_projection_method/> |
2711 | + </scheme> |
2712 | + <solver> |
2713 | + <iterative_method name="cg"/> |
2714 | + <preconditioner name="sor"/> |
2715 | + <relative_error> |
2716 | + <real_value rank="0">1.0e-10</real_value> |
2717 | + </relative_error> |
2718 | + <max_iterations> |
2719 | + <integer_value rank="0">1000</integer_value> |
2720 | + </max_iterations> |
2721 | + <never_ignore_solver_failures/> |
2722 | + <diagnostics> |
2723 | + <monitors/> |
2724 | + </diagnostics> |
2725 | + </solver> |
2726 | + <boundary_conditions name="right_outflow"> |
2727 | + <surface_ids> |
2728 | + <integer_value shape="1" rank="1">2</integer_value> |
2729 | + </surface_ids> |
2730 | + <type name="dirichlet"> |
2731 | + <constant> |
2732 | + <real_value rank="0">0.0</real_value> |
2733 | + </constant> |
2734 | + </type> |
2735 | + </boundary_conditions> |
2736 | + <output/> |
2737 | + <stat/> |
2738 | + <convergence> |
2739 | + <include_in_convergence/> |
2740 | + </convergence> |
2741 | + <detectors> |
2742 | + <exclude_from_detectors/> |
2743 | + </detectors> |
2744 | + <steady_state> |
2745 | + <include_in_steady_state/> |
2746 | + </steady_state> |
2747 | + <no_interpolation/> |
2748 | + </prognostic> |
2749 | + </scalar_field> |
2750 | + <vector_field name="Velocity" rank="1"> |
2751 | + <prognostic> |
2752 | + <mesh name="VelocityMesh"/> |
2753 | + <equation name="Boussinesq"/> |
2754 | + <spatial_discretisation> |
2755 | + <discontinuous_galerkin> |
2756 | + <mass_terms> |
2757 | + <exclude_mass_terms/> |
2758 | + </mass_terms> |
2759 | + <viscosity_scheme> |
2760 | + <compact_discontinuous_galerkin/> |
2761 | + </viscosity_scheme> |
2762 | + <advection_scheme> |
2763 | + <none/> |
2764 | + <integrate_advection_by_parts> |
2765 | + <twice/> |
2766 | + </integrate_advection_by_parts> |
2767 | + </advection_scheme> |
2768 | + </discontinuous_galerkin> |
2769 | + <conservative_advection> |
2770 | + <real_value rank="0">0.0</real_value> |
2771 | + </conservative_advection> |
2772 | + </spatial_discretisation> |
2773 | + <temporal_discretisation> |
2774 | + <theta> |
2775 | + <real_value rank="0">1.0</real_value> |
2776 | + </theta> |
2777 | + <relaxation> |
2778 | + <real_value rank="0">1.0</real_value> |
2779 | + </relaxation> |
2780 | + </temporal_discretisation> |
2781 | + <solver> |
2782 | + <iterative_method name="gmres"> |
2783 | + <restart> |
2784 | + <integer_value rank="0">30</integer_value> |
2785 | + </restart> |
2786 | + </iterative_method> |
2787 | + <preconditioner name="sor"/> |
2788 | + <relative_error> |
2789 | + <real_value rank="0">1.0e-10</real_value> |
2790 | + </relative_error> |
2791 | + <max_iterations> |
2792 | + <integer_value rank="0">1000</integer_value> |
2793 | + </max_iterations> |
2794 | + <never_ignore_solver_failures/> |
2795 | + <diagnostics> |
2796 | + <monitors/> |
2797 | + </diagnostics> |
2798 | + </solver> |
2799 | + <initial_condition name="WholeMesh"> |
2800 | + <constant> |
2801 | + <real_value shape="1" dim1="dim" rank="1">0.0</real_value> |
2802 | + </constant> |
2803 | + </initial_condition> |
2804 | + <boundary_conditions name="left_inflow"> |
2805 | + <surface_ids> |
2806 | + <integer_value shape="1" rank="1">1</integer_value> |
2807 | + </surface_ids> |
2808 | + <type name="dirichlet"> |
2809 | + <apply_weakly/> |
2810 | + <align_bc_with_cartesian> |
2811 | + <x_component> |
2812 | + <constant> |
2813 | + <real_value rank="0">5.0</real_value> |
2814 | + </constant> |
2815 | + </x_component> |
2816 | + </align_bc_with_cartesian> |
2817 | + </type> |
2818 | + </boundary_conditions> |
2819 | + <vector_field name="Absorption" rank="1"> |
2820 | + <diagnostic> |
2821 | + <mesh name="DGP0"/> |
2822 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
2823 | + <string_value lines="20" type="python"># get the prescribed fields |
2824 | +perm = states["fluid"].scalar_fields["Permeability"] |
2825 | +visc = states["fluid"].scalar_fields["Viscosity"] |
2826 | + |
2827 | +# loop the DMmesh_p0 nodes (which are element centred) |
2828 | +for n in range(field.node_count): |
2829 | + perm_n = perm.node_val(n) |
2830 | + visc_n = visc.node_val(n) |
2831 | + |
2832 | + # calc the absorption term |
2833 | + sigma_n = visc_n/perm_n |
2834 | + |
2835 | + # set the absorption term |
2836 | + field.set(n,sigma_n)</string_value> |
2837 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
2838 | + |
2839 | +also assume that the velocity is dg</comment> |
2840 | + </algorithm> |
2841 | + <output/> |
2842 | + <stat> |
2843 | + <include_in_stat/> |
2844 | + </stat> |
2845 | + <convergence> |
2846 | + <include_in_convergence/> |
2847 | + </convergence> |
2848 | + <detectors> |
2849 | + <include_in_detectors/> |
2850 | + </detectors> |
2851 | + <steady_state> |
2852 | + <include_in_steady_state/> |
2853 | + </steady_state> |
2854 | + </diagnostic> |
2855 | + <include_pressure_correction/> |
2856 | + </vector_field> |
2857 | + <output/> |
2858 | + <stat> |
2859 | + <include_in_stat/> |
2860 | + <previous_time_step> |
2861 | + <exclude_from_stat/> |
2862 | + </previous_time_step> |
2863 | + <nonlinear_field> |
2864 | + <exclude_from_stat/> |
2865 | + </nonlinear_field> |
2866 | + </stat> |
2867 | + <convergence> |
2868 | + <include_in_convergence/> |
2869 | + </convergence> |
2870 | + <detectors> |
2871 | + <include_in_detectors/> |
2872 | + </detectors> |
2873 | + <steady_state> |
2874 | + <include_in_steady_state/> |
2875 | + </steady_state> |
2876 | + <consistent_interpolation/> |
2877 | + </prognostic> |
2878 | + </vector_field> |
2879 | + <scalar_field name="Porosity" rank="0"> |
2880 | + <prescribed> |
2881 | + <mesh name="DGP0"/> |
2882 | + <value name="WholeMesh"> |
2883 | + <constant> |
2884 | + <real_value rank="0">0.5</real_value> |
2885 | + </constant> |
2886 | + </value> |
2887 | + <output/> |
2888 | + <stat/> |
2889 | + <detectors> |
2890 | + <exclude_from_detectors/> |
2891 | + </detectors> |
2892 | + </prescribed> |
2893 | + </scalar_field> |
2894 | + <scalar_field name="Permeability" rank="0"> |
2895 | + <prescribed> |
2896 | + <mesh name="DGP0"/> |
2897 | + <value name="WholeMesh"> |
2898 | + <constant> |
2899 | + <real_value rank="0">1.0e-10</real_value> |
2900 | + </constant> |
2901 | + </value> |
2902 | + <output/> |
2903 | + <stat/> |
2904 | + <detectors> |
2905 | + <exclude_from_detectors/> |
2906 | + </detectors> |
2907 | + </prescribed> |
2908 | + </scalar_field> |
2909 | + <scalar_field name="Viscosity" rank="0"> |
2910 | + <prescribed> |
2911 | + <mesh name="DGP0"/> |
2912 | + <value name="WholeMesh"> |
2913 | + <constant> |
2914 | + <real_value rank="0">1.0e-04</real_value> |
2915 | + </constant> |
2916 | + </value> |
2917 | + <output/> |
2918 | + <stat/> |
2919 | + <detectors> |
2920 | + <exclude_from_detectors/> |
2921 | + </detectors> |
2922 | + </prescribed> |
2923 | + </scalar_field> |
2924 | + <vector_field name="interstitial_velocity" rank="1"> |
2925 | + <diagnostic> |
2926 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
2927 | + <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"] |
2928 | +darcy_vel = states["fluid"].vector_fields["Velocity"] |
2929 | + |
2930 | +for ele in range(field.element_count): |
2931 | + phi_ele = phi.node_val(ele) |
2932 | + |
2933 | + factor_ele = 1.0/max(phi_ele,1.0e-08) |
2934 | + |
2935 | + vel_ele = [] |
2936 | + for i in range(len(darcy_vel.ele_val(ele))): |
2937 | + vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i]) |
2938 | + |
2939 | + field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value> |
2940 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
2941 | + |
2942 | +also assume that the velocity is dg</comment> |
2943 | + </algorithm> |
2944 | + <mesh name="VelocityMesh"/> |
2945 | + <output/> |
2946 | + <stat> |
2947 | + <include_in_stat/> |
2948 | + </stat> |
2949 | + <convergence> |
2950 | + <include_in_convergence/> |
2951 | + </convergence> |
2952 | + <detectors> |
2953 | + <include_in_detectors/> |
2954 | + </detectors> |
2955 | + <steady_state> |
2956 | + <include_in_steady_state/> |
2957 | + </steady_state> |
2958 | + </diagnostic> |
2959 | + </vector_field> |
2960 | + </material_phase> |
2961 | +</fluidity_options> |
2962 | |
2963 | === added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_2d.flml' |
2964 | --- tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_2d.flml 1970-01-01 00:00:00 +0000 |
2965 | +++ tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_2d.flml 2011-12-23 12:05:25 +0000 |
2966 | @@ -0,0 +1,407 @@ |
2967 | +<?xml version='1.0' encoding='utf-8'?> |
2968 | +<fluidity_options> |
2969 | + <simulation_name> |
2970 | + <string_value lines="1">darcy_p0p1_test_cty_cv_velBCinlet_2d</string_value> |
2971 | + <comment>a simple short test case for darcy flow using the element type p0p1_test_cty_cv for velocity-pressure |
2972 | + |
2973 | +the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution) |
2974 | + |
2975 | +it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct |
2976 | + |
2977 | +the darcy flow equation is |
2978 | + |
2979 | +sigma*darcy_vel = - grad P |
2980 | + |
2981 | +sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options) |
2982 | + |
2983 | +python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity |
2984 | + |
2985 | +this case has a 1 region 1 material with velocity boundary condition inflow. |
2986 | + |
2987 | +to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input. |
2988 | + |
2989 | +the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) |
2990 | + |
2991 | +the geometry is 2d (although this is a 1d problem) and one time step is done (as nothing depends on time)</comment> |
2992 | + </simulation_name> |
2993 | + <problem_type> |
2994 | + <string_value lines="1">fluids</string_value> |
2995 | + </problem_type> |
2996 | + <geometry> |
2997 | + <dimension> |
2998 | + <integer_value rank="0">2</integer_value> |
2999 | + </dimension> |
3000 | + <mesh name="CoordinateMesh"> |
3001 | + <from_file file_name="square"> |
3002 | + <format name="gmsh"/> |
3003 | + <stat> |
3004 | + <include_in_stat/> |
3005 | + </stat> |
3006 | + </from_file> |
3007 | + </mesh> |
3008 | + <mesh name="VelocityMesh"> |
3009 | + <from_mesh> |
3010 | + <mesh name="CoordinateMesh"/> |
3011 | + <mesh_shape> |
3012 | + <polynomial_degree> |
3013 | + <integer_value rank="0">0</integer_value> |
3014 | + </polynomial_degree> |
3015 | + </mesh_shape> |
3016 | + <mesh_continuity> |
3017 | + <string_value>discontinuous</string_value> |
3018 | + </mesh_continuity> |
3019 | + <stat> |
3020 | + <exclude_from_stat/> |
3021 | + </stat> |
3022 | + </from_mesh> |
3023 | + </mesh> |
3024 | + <mesh name="PressureMesh"> |
3025 | + <from_mesh> |
3026 | + <mesh name="CoordinateMesh"/> |
3027 | + <mesh_shape> |
3028 | + <polynomial_degree> |
3029 | + <integer_value rank="0">1</integer_value> |
3030 | + </polynomial_degree> |
3031 | + </mesh_shape> |
3032 | + <stat> |
3033 | + <exclude_from_stat/> |
3034 | + </stat> |
3035 | + </from_mesh> |
3036 | + </mesh> |
3037 | + <mesh name="DGP0"> |
3038 | + <from_mesh> |
3039 | + <mesh name="CoordinateMesh"/> |
3040 | + <mesh_shape> |
3041 | + <polynomial_degree> |
3042 | + <integer_value rank="0">0</integer_value> |
3043 | + </polynomial_degree> |
3044 | + </mesh_shape> |
3045 | + <mesh_continuity> |
3046 | + <string_value>discontinuous</string_value> |
3047 | + </mesh_continuity> |
3048 | + <stat> |
3049 | + <exclude_from_stat/> |
3050 | + </stat> |
3051 | + </from_mesh> |
3052 | + </mesh> |
3053 | + <mesh name="OutputMesh"> |
3054 | + <from_mesh> |
3055 | + <mesh name="CoordinateMesh"/> |
3056 | + <mesh_continuity> |
3057 | + <string_value>discontinuous</string_value> |
3058 | + </mesh_continuity> |
3059 | + <stat> |
3060 | + <exclude_from_stat/> |
3061 | + </stat> |
3062 | + </from_mesh> |
3063 | + </mesh> |
3064 | + <quadrature> |
3065 | + <degree> |
3066 | + <integer_value rank="0">2</integer_value> |
3067 | + </degree> |
3068 | + </quadrature> |
3069 | + </geometry> |
3070 | + <io> |
3071 | + <dump_format> |
3072 | + <string_value>vtk</string_value> |
3073 | + </dump_format> |
3074 | + <dump_period_in_timesteps> |
3075 | + <constant> |
3076 | + <integer_value rank="0">0</integer_value> |
3077 | + </constant> |
3078 | + </dump_period_in_timesteps> |
3079 | + <output_mesh name="OutputMesh"/> |
3080 | + <stat/> |
3081 | + </io> |
3082 | + <timestepping> |
3083 | + <current_time> |
3084 | + <real_value rank="0">0.0</real_value> |
3085 | + </current_time> |
3086 | + <timestep> |
3087 | + <real_value rank="0">1.0</real_value> |
3088 | + </timestep> |
3089 | + <finish_time> |
3090 | + <real_value rank="0">1.0</real_value> |
3091 | + </finish_time> |
3092 | + </timestepping> |
3093 | + <material_phase name="fluid"> |
3094 | + <scalar_field name="Pressure" rank="0"> |
3095 | + <prognostic> |
3096 | + <mesh name="PressureMesh"/> |
3097 | + <spatial_discretisation> |
3098 | + <continuous_galerkin> |
3099 | + <test_continuity_with_cv_dual/> |
3100 | + </continuous_galerkin> |
3101 | + </spatial_discretisation> |
3102 | + <scheme> |
3103 | + <poisson_pressure_solution> |
3104 | + <string_value lines="1">never</string_value> |
3105 | + </poisson_pressure_solution> |
3106 | + <use_projection_method/> |
3107 | + </scheme> |
3108 | + <solver> |
3109 | + <iterative_method name="cg"/> |
3110 | + <preconditioner name="sor"/> |
3111 | + <relative_error> |
3112 | + <real_value rank="0">1.0e-10</real_value> |
3113 | + </relative_error> |
3114 | + <max_iterations> |
3115 | + <integer_value rank="0">1000</integer_value> |
3116 | + </max_iterations> |
3117 | + <never_ignore_solver_failures/> |
3118 | + <diagnostics> |
3119 | + <monitors/> |
3120 | + </diagnostics> |
3121 | + </solver> |
3122 | + <boundary_conditions name="right_outflow"> |
3123 | + <surface_ids> |
3124 | + <integer_value shape="1" rank="1">8</integer_value> |
3125 | + </surface_ids> |
3126 | + <type name="dirichlet"> |
3127 | + <apply_weakly/> |
3128 | + <constant> |
3129 | + <real_value rank="0">0.0</real_value> |
3130 | + </constant> |
3131 | + </type> |
3132 | + </boundary_conditions> |
3133 | + <output/> |
3134 | + <stat/> |
3135 | + <convergence> |
3136 | + <include_in_convergence/> |
3137 | + </convergence> |
3138 | + <detectors> |
3139 | + <exclude_from_detectors/> |
3140 | + </detectors> |
3141 | + <steady_state> |
3142 | + <include_in_steady_state/> |
3143 | + </steady_state> |
3144 | + <no_interpolation/> |
3145 | + </prognostic> |
3146 | + </scalar_field> |
3147 | + <vector_field name="Velocity" rank="1"> |
3148 | + <prognostic> |
3149 | + <mesh name="VelocityMesh"/> |
3150 | + <equation name="Boussinesq"/> |
3151 | + <spatial_discretisation> |
3152 | + <discontinuous_galerkin> |
3153 | + <mass_terms> |
3154 | + <exclude_mass_terms/> |
3155 | + </mass_terms> |
3156 | + <viscosity_scheme> |
3157 | + <compact_discontinuous_galerkin/> |
3158 | + </viscosity_scheme> |
3159 | + <advection_scheme> |
3160 | + <none/> |
3161 | + <integrate_advection_by_parts> |
3162 | + <twice/> |
3163 | + </integrate_advection_by_parts> |
3164 | + </advection_scheme> |
3165 | + </discontinuous_galerkin> |
3166 | + <conservative_advection> |
3167 | + <real_value rank="0">0.0</real_value> |
3168 | + </conservative_advection> |
3169 | + </spatial_discretisation> |
3170 | + <temporal_discretisation> |
3171 | + <theta> |
3172 | + <real_value rank="0">1.0</real_value> |
3173 | + </theta> |
3174 | + <relaxation> |
3175 | + <real_value rank="0">1.0</real_value> |
3176 | + </relaxation> |
3177 | + </temporal_discretisation> |
3178 | + <solver> |
3179 | + <iterative_method name="gmres"> |
3180 | + <restart> |
3181 | + <integer_value rank="0">30</integer_value> |
3182 | + </restart> |
3183 | + </iterative_method> |
3184 | + <preconditioner name="sor"/> |
3185 | + <relative_error> |
3186 | + <real_value rank="0">1.0e-10</real_value> |
3187 | + </relative_error> |
3188 | + <max_iterations> |
3189 | + <integer_value rank="0">1000</integer_value> |
3190 | + </max_iterations> |
3191 | + <never_ignore_solver_failures/> |
3192 | + <diagnostics> |
3193 | + <monitors/> |
3194 | + </diagnostics> |
3195 | + </solver> |
3196 | + <initial_condition name="WholeMesh"> |
3197 | + <constant> |
3198 | + <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value> |
3199 | + </constant> |
3200 | + </initial_condition> |
3201 | + <boundary_conditions name="left_inflow"> |
3202 | + <surface_ids> |
3203 | + <integer_value shape="1" rank="1">10</integer_value> |
3204 | + </surface_ids> |
3205 | + <type name="dirichlet"> |
3206 | + <apply_weakly/> |
3207 | + <align_bc_with_cartesian> |
3208 | + <x_component> |
3209 | + <constant> |
3210 | + <real_value rank="0">5.0</real_value> |
3211 | + </constant> |
3212 | + </x_component> |
3213 | + </align_bc_with_cartesian> |
3214 | + </type> |
3215 | + </boundary_conditions> |
3216 | + <boundary_conditions name="no_outflow_top_bottom"> |
3217 | + <surface_ids> |
3218 | + <integer_value shape="2" rank="1">7 9</integer_value> |
3219 | + </surface_ids> |
3220 | + <type name="dirichlet"> |
3221 | + <apply_weakly/> |
3222 | + <align_bc_with_cartesian> |
3223 | + <y_component> |
3224 | + <constant> |
3225 | + <real_value rank="0">0.0</real_value> |
3226 | + </constant> |
3227 | + </y_component> |
3228 | + </align_bc_with_cartesian> |
3229 | + </type> |
3230 | + </boundary_conditions> |
3231 | + <vector_field name="Absorption" rank="1"> |
3232 | + <diagnostic> |
3233 | + <mesh name="DGP0"/> |
3234 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
3235 | + <string_value lines="20" type="python"># get the prescribed fields |
3236 | +perm = states["fluid"].scalar_fields["Permeability"] |
3237 | +visc = states["fluid"].scalar_fields["Viscosity"] |
3238 | + |
3239 | +# loop the DMmesh_p0 nodes (which are element centred) |
3240 | +for n in range(field.node_count): |
3241 | + perm_n = perm.node_val(n) |
3242 | + visc_n = visc.node_val(n) |
3243 | + |
3244 | + # calc the absorption term |
3245 | + sigma_n = visc_n/perm_n |
3246 | + |
3247 | + # set the absorption term |
3248 | + field.set(n,sigma_n)</string_value> |
3249 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
3250 | + |
3251 | +also assume that the velocity is dg</comment> |
3252 | + </algorithm> |
3253 | + <output/> |
3254 | + <stat> |
3255 | + <include_in_stat/> |
3256 | + </stat> |
3257 | + <convergence> |
3258 | + <include_in_convergence/> |
3259 | + </convergence> |
3260 | + <detectors> |
3261 | + <include_in_detectors/> |
3262 | + </detectors> |
3263 | + <steady_state> |
3264 | + <include_in_steady_state/> |
3265 | + </steady_state> |
3266 | + </diagnostic> |
3267 | + <include_pressure_correction/> |
3268 | + </vector_field> |
3269 | + <output/> |
3270 | + <stat> |
3271 | + <include_in_stat/> |
3272 | + <previous_time_step> |
3273 | + <exclude_from_stat/> |
3274 | + </previous_time_step> |
3275 | + <nonlinear_field> |
3276 | + <exclude_from_stat/> |
3277 | + </nonlinear_field> |
3278 | + </stat> |
3279 | + <convergence> |
3280 | + <include_in_convergence/> |
3281 | + </convergence> |
3282 | + <detectors> |
3283 | + <include_in_detectors/> |
3284 | + </detectors> |
3285 | + <steady_state> |
3286 | + <include_in_steady_state/> |
3287 | + </steady_state> |
3288 | + <consistent_interpolation/> |
3289 | + </prognostic> |
3290 | + </vector_field> |
3291 | + <scalar_field name="Porosity" rank="0"> |
3292 | + <prescribed> |
3293 | + <mesh name="DGP0"/> |
3294 | + <value name="WholeMesh"> |
3295 | + <constant> |
3296 | + <real_value rank="0">0.5</real_value> |
3297 | + </constant> |
3298 | + </value> |
3299 | + <output/> |
3300 | + <stat/> |
3301 | + <detectors> |
3302 | + <exclude_from_detectors/> |
3303 | + </detectors> |
3304 | + </prescribed> |
3305 | + </scalar_field> |
3306 | + <scalar_field name="Permeability" rank="0"> |
3307 | + <prescribed> |
3308 | + <mesh name="DGP0"/> |
3309 | + <value name="WholeMesh"> |
3310 | + <constant> |
3311 | + <real_value rank="0">1.0e-10</real_value> |
3312 | + </constant> |
3313 | + </value> |
3314 | + <output/> |
3315 | + <stat/> |
3316 | + <detectors> |
3317 | + <exclude_from_detectors/> |
3318 | + </detectors> |
3319 | + </prescribed> |
3320 | + </scalar_field> |
3321 | + <scalar_field name="Viscosity" rank="0"> |
3322 | + <prescribed> |
3323 | + <mesh name="DGP0"/> |
3324 | + <value name="WholeMesh"> |
3325 | + <constant> |
3326 | + <real_value rank="0">1.0e-04</real_value> |
3327 | + </constant> |
3328 | + </value> |
3329 | + <output/> |
3330 | + <stat/> |
3331 | + <detectors> |
3332 | + <exclude_from_detectors/> |
3333 | + </detectors> |
3334 | + </prescribed> |
3335 | + </scalar_field> |
3336 | + <vector_field name="interstitial_velocity" rank="1"> |
3337 | + <diagnostic> |
3338 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
3339 | + <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"] |
3340 | +darcy_vel = states["fluid"].vector_fields["Velocity"] |
3341 | + |
3342 | +for ele in range(field.element_count): |
3343 | + phi_ele = phi.node_val(ele) |
3344 | + |
3345 | + factor_ele = 1.0/max(phi_ele,1.0e-08) |
3346 | + |
3347 | + vel_ele = [] |
3348 | + for i in range(len(darcy_vel.ele_val(ele))): |
3349 | + vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i]) |
3350 | + |
3351 | + field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value> |
3352 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
3353 | + |
3354 | +also assume that the velocity is dg</comment> |
3355 | + </algorithm> |
3356 | + <mesh name="VelocityMesh"/> |
3357 | + <output/> |
3358 | + <stat> |
3359 | + <include_in_stat/> |
3360 | + </stat> |
3361 | + <convergence> |
3362 | + <include_in_convergence/> |
3363 | + </convergence> |
3364 | + <detectors> |
3365 | + <include_in_detectors/> |
3366 | + </detectors> |
3367 | + <steady_state> |
3368 | + <include_in_steady_state/> |
3369 | + </steady_state> |
3370 | + </diagnostic> |
3371 | + </vector_field> |
3372 | + </material_phase> |
3373 | +</fluidity_options> |
3374 | |
3375 | === added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_3d.flml' |
3376 | --- tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_3d.flml 1970-01-01 00:00:00 +0000 |
3377 | +++ tests/darcy_p0p1_test_cty_cv_velBCinlet/darcy_p0p1_test_cty_cv_velBCinlet_3d.flml 2011-12-23 12:05:25 +0000 |
3378 | @@ -0,0 +1,404 @@ |
3379 | +<?xml version='1.0' encoding='utf-8'?> |
3380 | +<fluidity_options> |
3381 | + <simulation_name> |
3382 | + <string_value lines="1">darcy_p0p1_test_cty_cv_velBCinlet_3d</string_value> |
3383 | + <comment>a simple short test case for darcy flow using the element type p0p1_test_cty_cv for velocity-pressure |
3384 | + |
3385 | +the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution) |
3386 | + |
3387 | +it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct |
3388 | + |
3389 | +the darcy flow equation is |
3390 | + |
3391 | +sigma*darcy_vel = - grad P |
3392 | + |
3393 | +sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options) |
3394 | + |
3395 | +python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity |
3396 | + |
3397 | +this case has a 1 region 1 material with velocity boundary condition inflow. |
3398 | + |
3399 | +to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input. |
3400 | + |
3401 | +the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) |
3402 | + |
3403 | +the geometry is 3d (although this is a 1d problem) and one time step is done (as nothing depends on time)</comment> |
3404 | + </simulation_name> |
3405 | + <problem_type> |
3406 | + <string_value lines="1">fluids</string_value> |
3407 | + </problem_type> |
3408 | + <geometry> |
3409 | + <dimension> |
3410 | + <integer_value rank="0">3</integer_value> |
3411 | + </dimension> |
3412 | + <mesh name="CoordinateMesh"> |
3413 | + <from_file file_name="cube"> |
3414 | + <format name="gmsh"/> |
3415 | + <stat> |
3416 | + <include_in_stat/> |
3417 | + </stat> |
3418 | + </from_file> |
3419 | + </mesh> |
3420 | + <mesh name="VelocityMesh"> |
3421 | + <from_mesh> |
3422 | + <mesh name="CoordinateMesh"/> |
3423 | + <mesh_shape> |
3424 | + <polynomial_degree> |
3425 | + <integer_value rank="0">0</integer_value> |
3426 | + </polynomial_degree> |
3427 | + </mesh_shape> |
3428 | + <mesh_continuity> |
3429 | + <string_value>discontinuous</string_value> |
3430 | + </mesh_continuity> |
3431 | + <stat> |
3432 | + <exclude_from_stat/> |
3433 | + </stat> |
3434 | + </from_mesh> |
3435 | + </mesh> |
3436 | + <mesh name="PressureMesh"> |
3437 | + <from_mesh> |
3438 | + <mesh name="CoordinateMesh"/> |
3439 | + <mesh_shape> |
3440 | + <polynomial_degree> |
3441 | + <integer_value rank="0">1</integer_value> |
3442 | + </polynomial_degree> |
3443 | + </mesh_shape> |
3444 | + <stat> |
3445 | + <exclude_from_stat/> |
3446 | + </stat> |
3447 | + </from_mesh> |
3448 | + </mesh> |
3449 | + <mesh name="DGP0"> |
3450 | + <from_mesh> |
3451 | + <mesh name="CoordinateMesh"/> |
3452 | + <mesh_shape> |
3453 | + <polynomial_degree> |
3454 | + <integer_value rank="0">0</integer_value> |
3455 | + </polynomial_degree> |
3456 | + </mesh_shape> |
3457 | + <mesh_continuity> |
3458 | + <string_value>discontinuous</string_value> |
3459 | + </mesh_continuity> |
3460 | + <stat> |
3461 | + <exclude_from_stat/> |
3462 | + </stat> |
3463 | + </from_mesh> |
3464 | + </mesh> |
3465 | + <mesh name="OutputMesh"> |
3466 | + <from_mesh> |
3467 | + <mesh name="CoordinateMesh"/> |
3468 | + <mesh_continuity> |
3469 | + <string_value>discontinuous</string_value> |
3470 | + </mesh_continuity> |
3471 | + <stat> |
3472 | + <exclude_from_stat/> |
3473 | + </stat> |
3474 | + </from_mesh> |
3475 | + </mesh> |
3476 | + <quadrature> |
3477 | + <degree> |
3478 | + <integer_value rank="0">2</integer_value> |
3479 | + </degree> |
3480 | + </quadrature> |
3481 | + </geometry> |
3482 | + <io> |
3483 | + <dump_format> |
3484 | + <string_value>vtk</string_value> |
3485 | + </dump_format> |
3486 | + <dump_period_in_timesteps> |
3487 | + <constant> |
3488 | + <integer_value rank="0">0</integer_value> |
3489 | + </constant> |
3490 | + </dump_period_in_timesteps> |
3491 | + <output_mesh name="OutputMesh"/> |
3492 | + <stat/> |
3493 | + </io> |
3494 | + <timestepping> |
3495 | + <current_time> |
3496 | + <real_value rank="0">0.0</real_value> |
3497 | + </current_time> |
3498 | + <timestep> |
3499 | + <real_value rank="0">1.0</real_value> |
3500 | + </timestep> |
3501 | + <finish_time> |
3502 | + <real_value rank="0">1.0</real_value> |
3503 | + </finish_time> |
3504 | + </timestepping> |
3505 | + <material_phase name="fluid"> |
3506 | + <scalar_field name="Pressure" rank="0"> |
3507 | + <prognostic> |
3508 | + <mesh name="PressureMesh"/> |
3509 | + <spatial_discretisation> |
3510 | + <continuous_galerkin> |
3511 | + <test_continuity_with_cv_dual/> |
3512 | + </continuous_galerkin> |
3513 | + </spatial_discretisation> |
3514 | + <scheme> |
3515 | + <poisson_pressure_solution> |
3516 | + <string_value lines="1">never</string_value> |
3517 | + </poisson_pressure_solution> |
3518 | + <use_projection_method/> |
3519 | + </scheme> |
3520 | + <solver> |
3521 | + <iterative_method name="cg"/> |
3522 | + <preconditioner name="sor"/> |
3523 | + <relative_error> |
3524 | + <real_value rank="0">1.0e-10</real_value> |
3525 | + </relative_error> |
3526 | + <max_iterations> |
3527 | + <integer_value rank="0">1000</integer_value> |
3528 | + </max_iterations> |
3529 | + <never_ignore_solver_failures/> |
3530 | + <diagnostics> |
3531 | + <monitors/> |
3532 | + </diagnostics> |
3533 | + </solver> |
3534 | + <boundary_conditions name="right_outflow"> |
3535 | + <surface_ids> |
3536 | + <integer_value shape="1" rank="1">4</integer_value> |
3537 | + </surface_ids> |
3538 | + <type name="dirichlet"> |
3539 | + <apply_weakly/> |
3540 | + <constant> |
3541 | + <real_value rank="0">0.0</real_value> |
3542 | + </constant> |
3543 | + </type> |
3544 | + </boundary_conditions> |
3545 | + <output/> |
3546 | + <stat/> |
3547 | + <convergence> |
3548 | + <include_in_convergence/> |
3549 | + </convergence> |
3550 | + <detectors> |
3551 | + <exclude_from_detectors/> |
3552 | + </detectors> |
3553 | + <steady_state> |
3554 | + <include_in_steady_state/> |
3555 | + </steady_state> |
3556 | + <no_interpolation/> |
3557 | + </prognostic> |
3558 | + </scalar_field> |
3559 | + <vector_field name="Velocity" rank="1"> |
3560 | + <prognostic> |
3561 | + <mesh name="VelocityMesh"/> |
3562 | + <equation name="Boussinesq"/> |
3563 | + <spatial_discretisation> |
3564 | + <discontinuous_galerkin> |
3565 | + <mass_terms> |
3566 | + <exclude_mass_terms/> |
3567 | + </mass_terms> |
3568 | + <viscosity_scheme> |
3569 | + <compact_discontinuous_galerkin/> |
3570 | + </viscosity_scheme> |
3571 | + <advection_scheme> |
3572 | + <none/> |
3573 | + <integrate_advection_by_parts> |
3574 | + <twice/> |
3575 | + </integrate_advection_by_parts> |
3576 | + </advection_scheme> |
3577 | + </discontinuous_galerkin> |
3578 | + <conservative_advection> |
3579 | + <real_value rank="0">0.0</real_value> |
3580 | + </conservative_advection> |
3581 | + </spatial_discretisation> |
3582 | + <temporal_discretisation> |
3583 | + <theta> |
3584 | + <real_value rank="0">1.0</real_value> |
3585 | + </theta> |
3586 | + <relaxation> |
3587 | + <real_value rank="0">1.0</real_value> |
3588 | + </relaxation> |
3589 | + </temporal_discretisation> |
3590 | + <solver> |
3591 | + <iterative_method name="gmres"> |
3592 | + <restart> |
3593 | + <integer_value rank="0">30</integer_value> |
3594 | + </restart> |
3595 | + </iterative_method> |
3596 | + <preconditioner name="sor"/> |
3597 | + <relative_error> |
3598 | + <real_value rank="0">1.0e-10</real_value> |
3599 | + </relative_error> |
3600 | + <max_iterations> |
3601 | + <integer_value rank="0">1000</integer_value> |
3602 | + </max_iterations> |
3603 | + <never_ignore_solver_failures/> |
3604 | + <diagnostics> |
3605 | + <monitors/> |
3606 | + </diagnostics> |
3607 | + </solver> |
3608 | + <initial_condition name="WholeMesh"> |
3609 | + <constant> |
3610 | + <real_value shape="3" dim1="dim" rank="1">0.0 0.0 0.0</real_value> |
3611 | + </constant> |
3612 | + </initial_condition> |
3613 | + <boundary_conditions name="left_inflow"> |
3614 | + <surface_ids> |
3615 | + <integer_value shape="1" rank="1">3</integer_value> |
3616 | + </surface_ids> |
3617 | + <type name="dirichlet"> |
3618 | + <apply_weakly/> |
3619 | + <align_bc_with_cartesian> |
3620 | + <x_component> |
3621 | + <constant> |
3622 | + <real_value rank="0">5.0</real_value> |
3623 | + </constant> |
3624 | + </x_component> |
3625 | + </align_bc_with_cartesian> |
3626 | + </type> |
3627 | + </boundary_conditions> |
3628 | + <boundary_conditions name="no_outflow_top_bottom"> |
3629 | + <surface_ids> |
3630 | + <integer_value shape="2" rank="1">1 2</integer_value> |
3631 | + </surface_ids> |
3632 | + <type name="no_normal_flow"/> |
3633 | + </boundary_conditions> |
3634 | + <boundary_conditions name="no_outflow_front_back"> |
3635 | + <surface_ids> |
3636 | + <integer_value shape="2" rank="1">5 6</integer_value> |
3637 | + </surface_ids> |
3638 | + <type name="no_normal_flow"/> |
3639 | + </boundary_conditions> |
3640 | + <vector_field name="Absorption" rank="1"> |
3641 | + <diagnostic> |
3642 | + <mesh name="DGP0"/> |
3643 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
3644 | + <string_value lines="20" type="python"># get the prescribed fields |
3645 | +perm = states["fluid"].scalar_fields["Permeability"] |
3646 | +visc = states["fluid"].scalar_fields["Viscosity"] |
3647 | + |
3648 | +# loop the DMmesh_p0 nodes (which are element centred) |
3649 | +for n in range(field.node_count): |
3650 | + perm_n = perm.node_val(n) |
3651 | + visc_n = visc.node_val(n) |
3652 | + |
3653 | + # calc the absorption term |
3654 | + sigma_n = visc_n/perm_n |
3655 | + |
3656 | + # set the absorption term |
3657 | + field.set(n,sigma_n)</string_value> |
3658 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
3659 | + |
3660 | +also assume that the velocity is dg</comment> |
3661 | + </algorithm> |
3662 | + <output/> |
3663 | + <stat> |
3664 | + <include_in_stat/> |
3665 | + </stat> |
3666 | + <convergence> |
3667 | + <include_in_convergence/> |
3668 | + </convergence> |
3669 | + <detectors> |
3670 | + <include_in_detectors/> |
3671 | + </detectors> |
3672 | + <steady_state> |
3673 | + <include_in_steady_state/> |
3674 | + </steady_state> |
3675 | + </diagnostic> |
3676 | + <include_pressure_correction/> |
3677 | + </vector_field> |
3678 | + <output/> |
3679 | + <stat> |
3680 | + <include_in_stat/> |
3681 | + <previous_time_step> |
3682 | + <exclude_from_stat/> |
3683 | + </previous_time_step> |
3684 | + <nonlinear_field> |
3685 | + <exclude_from_stat/> |
3686 | + </nonlinear_field> |
3687 | + </stat> |
3688 | + <convergence> |
3689 | + <include_in_convergence/> |
3690 | + </convergence> |
3691 | + <detectors> |
3692 | + <include_in_detectors/> |
3693 | + </detectors> |
3694 | + <steady_state> |
3695 | + <include_in_steady_state/> |
3696 | + </steady_state> |
3697 | + <consistent_interpolation/> |
3698 | + </prognostic> |
3699 | + </vector_field> |
3700 | + <scalar_field name="Porosity" rank="0"> |
3701 | + <prescribed> |
3702 | + <mesh name="DGP0"/> |
3703 | + <value name="WholeMesh"> |
3704 | + <constant> |
3705 | + <real_value rank="0">0.5</real_value> |
3706 | + </constant> |
3707 | + </value> |
3708 | + <output/> |
3709 | + <stat/> |
3710 | + <detectors> |
3711 | + <exclude_from_detectors/> |
3712 | + </detectors> |
3713 | + </prescribed> |
3714 | + </scalar_field> |
3715 | + <scalar_field name="Permeability" rank="0"> |
3716 | + <prescribed> |
3717 | + <mesh name="DGP0"/> |
3718 | + <value name="WholeMesh"> |
3719 | + <constant> |
3720 | + <real_value rank="0">1.0e-10</real_value> |
3721 | + </constant> |
3722 | + </value> |
3723 | + <output/> |
3724 | + <stat/> |
3725 | + <detectors> |
3726 | + <exclude_from_detectors/> |
3727 | + </detectors> |
3728 | + </prescribed> |
3729 | + </scalar_field> |
3730 | + <scalar_field name="Viscosity" rank="0"> |
3731 | + <prescribed> |
3732 | + <mesh name="DGP0"/> |
3733 | + <value name="WholeMesh"> |
3734 | + <constant> |
3735 | + <real_value rank="0">1.0e-04</real_value> |
3736 | + </constant> |
3737 | + </value> |
3738 | + <output/> |
3739 | + <stat/> |
3740 | + <detectors> |
3741 | + <exclude_from_detectors/> |
3742 | + </detectors> |
3743 | + </prescribed> |
3744 | + </scalar_field> |
3745 | + <vector_field name="interstitial_velocity" rank="1"> |
3746 | + <diagnostic> |
3747 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
3748 | + <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"] |
3749 | +darcy_vel = states["fluid"].vector_fields["Velocity"] |
3750 | + |
3751 | +for ele in range(field.element_count): |
3752 | + phi_ele = phi.node_val(ele) |
3753 | + |
3754 | + factor_ele = 1.0/max(phi_ele,1.0e-08) |
3755 | + |
3756 | + vel_ele = [] |
3757 | + for i in range(len(darcy_vel.ele_val(ele))): |
3758 | + vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i]) |
3759 | + |
3760 | + field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value> |
3761 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
3762 | + |
3763 | +also assume that the velocity is dg</comment> |
3764 | + </algorithm> |
3765 | + <mesh name="VelocityMesh"/> |
3766 | + <output/> |
3767 | + <stat> |
3768 | + <include_in_stat/> |
3769 | + </stat> |
3770 | + <convergence> |
3771 | + <include_in_convergence/> |
3772 | + </convergence> |
3773 | + <detectors> |
3774 | + <include_in_detectors/> |
3775 | + </detectors> |
3776 | + <steady_state> |
3777 | + <include_in_steady_state/> |
3778 | + </steady_state> |
3779 | + </diagnostic> |
3780 | + </vector_field> |
3781 | + </material_phase> |
3782 | +</fluidity_options> |
3783 | |
3784 | === added file 'tests/darcy_p0p1_test_cty_cv_velBCinlet/square.geo' |
3785 | --- tests/darcy_p0p1_test_cty_cv_velBCinlet/square.geo 1970-01-01 00:00:00 +0000 |
3786 | +++ tests/darcy_p0p1_test_cty_cv_velBCinlet/square.geo 2011-12-23 12:05:25 +0000 |
3787 | @@ -0,0 +1,22 @@ |
3788 | +// define a layer variable |
3789 | +lay = 5; |
3790 | + |
3791 | +// define a len variable |
3792 | +len = 300.0; |
3793 | + |
3794 | +Point(1) = {0.0, 0.0, 0.0, 1.0}; |
3795 | + |
3796 | +Extrude {len, 0.0, 0.0} { |
3797 | + Point{1}; Layers{lay}; |
3798 | +} |
3799 | + |
3800 | +Extrude {0.0, len, 0.0} { |
3801 | + Line{1}; Layers{lay}; |
3802 | +} |
3803 | + |
3804 | +Physical Line(7) = {1}; |
3805 | +Physical Line(8) = {4}; |
3806 | +Physical Line(9) = {2}; |
3807 | +Physical Line(10) = {3}; |
3808 | + |
3809 | +Physical Surface(11) = {5}; |
3810 | |
3811 | === added directory 'tests/darcy_p0p1cv_pressBCinlet' |
3812 | === added file 'tests/darcy_p0p1cv_pressBCinlet/Makefile' |
3813 | --- tests/darcy_p0p1cv_pressBCinlet/Makefile 1970-01-01 00:00:00 +0000 |
3814 | +++ tests/darcy_p0p1cv_pressBCinlet/Makefile 2011-12-23 12:05:25 +0000 |
3815 | @@ -0,0 +1,32 @@ |
3816 | +SHELL = sh |
3817 | + |
3818 | +SIM = darcy_p0p1cv_pressBCinlet |
3819 | + |
3820 | +# options for 1d mesh interval script |
3821 | +MESH_SIZE = 100.0 |
3822 | +LEFT_X = 0.0 |
3823 | +RIGHT_X = 300.0 |
3824 | + |
3825 | +input: clean |
3826 | + interval --dx=$(MESH_SIZE) $(LEFT_X) $(RIGHT_X) $(SIM)_1d |
3827 | + gmsh -2 square.geo -o square.msh |
3828 | + gmsh -3 cube.geo -o cube.msh |
3829 | + |
3830 | +clean: |
3831 | + rm -f *.ele |
3832 | + rm -f *.node |
3833 | + rm -f *.bound |
3834 | + rm -f *.edge |
3835 | + rm -f *.face |
3836 | + rm -f *.vtu |
3837 | + rm -f *.pvtu |
3838 | + rm -f *.s |
3839 | + rm -f *.stat |
3840 | + rm -f *.log-0 |
3841 | + rm -f *.err-0 |
3842 | + rm -f *.msh |
3843 | + rm -f *.halo |
3844 | + rm -f fluidity.err* |
3845 | + rm -f fluidity.log* |
3846 | + rm -f matrixdump* |
3847 | + rm -f first_timestep_adapted_mesh* |
3848 | |
3849 | === added file 'tests/darcy_p0p1cv_pressBCinlet/cube.geo' |
3850 | --- tests/darcy_p0p1cv_pressBCinlet/cube.geo 1970-01-01 00:00:00 +0000 |
3851 | +++ tests/darcy_p0p1cv_pressBCinlet/cube.geo 2011-12-23 12:05:25 +0000 |
3852 | @@ -0,0 +1,39 @@ |
3853 | +// define a layer variable |
3854 | +lay = 5; |
3855 | + |
3856 | +// define a len variable |
3857 | +len = 300.0; |
3858 | + |
3859 | +Point(1) = {0.0, 0.0, 0.0, 1.0}; |
3860 | + |
3861 | +Extrude {len, 0.0, 0.0} { |
3862 | + Point{1}; Layers{lay}; |
3863 | +} |
3864 | + |
3865 | +Extrude {0.0, len, 0.0} { |
3866 | + Line{1}; Layers{lay}; |
3867 | +} |
3868 | + |
3869 | +Extrude {0.0, 0.0, len} { |
3870 | + Surface{5}; Layers{lay}; |
3871 | +} |
3872 | + |
3873 | +// bottom |
3874 | +Physical Surface(1) = {5}; |
3875 | + |
3876 | +// top |
3877 | +Physical Surface(2) = {27}; |
3878 | + |
3879 | +// left |
3880 | +Physical Surface(3) = {26}; |
3881 | + |
3882 | +// right |
3883 | +Physical Surface(4) = {18}; |
3884 | + |
3885 | +// front |
3886 | +Physical Surface(5) = {14}; |
3887 | + |
3888 | +// back |
3889 | +Physical Surface(6) = {22}; |
3890 | + |
3891 | +Physical Volume(1) = {1}; |
3892 | |
3893 | === added file 'tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet.xml' |
3894 | --- tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet.xml 1970-01-01 00:00:00 +0000 |
3895 | +++ tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet.xml 2011-12-23 12:05:25 +0000 |
3896 | @@ -0,0 +1,94 @@ |
3897 | +<?xml version="1.0" encoding="UTF-8" ?> |
3898 | +<!DOCTYPE testproblem SYSTEM "regressiontest.dtd"> |
3899 | + |
3900 | +<testproblem> |
3901 | + <name>darcy_p0p1cv_pressBCinlet</name> |
3902 | + <owner userid="btollit"/> |
3903 | + <tags>flml darcy</tags> |
3904 | + <problem_definition length="short" nprocs="1"> |
3905 | + <command_line> |
3906 | +../../bin/fluidity darcy_p0p1cv_pressBCinlet_1d.flml |
3907 | +../../bin/fluidity darcy_p0p1cv_pressBCinlet_2d.flml |
3908 | +../../bin/fluidity darcy_p0p1cv_pressBCinlet_3d.flml |
3909 | + </command_line> |
3910 | + <!-- One/two/three dimensional problem for darcy flow with one region and one material using p0p1cv element type testing the pressure gradient against analytic and the interstitial velocity. This tests using a diagnostic velocity absorption field associated with a mesh and momentum dg for darcy flow. This also tests the use of pressure BC for CV pressure.--> |
3911 | + </problem_definition> |
3912 | + <variables> |
3913 | + <variable name="pressure_1d" language="python"> |
3914 | +import vtktools |
3915 | +v = vtktools.vtu("darcy_p0p1cv_pressBCinlet_1d_1.vtu") |
3916 | +pressure_1d = v.GetScalarRange("Pressure") |
3917 | + </variable> |
3918 | + <variable name="inter_vel_1d" language="python"> |
3919 | +import vtktools |
3920 | +v = vtktools.vtu("darcy_p0p1cv_pressBCinlet_1d_1.vtu") |
3921 | +inter_vel_1d = max(v.GetScalarRange("interstitial_velocity")) |
3922 | + </variable> |
3923 | + <variable name="pressure_2d" language="python"> |
3924 | +import vtktools |
3925 | +v = vtktools.vtu("darcy_p0p1cv_pressBCinlet_2d_1.vtu") |
3926 | +pressure_2d = v.GetScalarRange("Pressure") |
3927 | + </variable> |
3928 | + <variable name="inter_vel_2d" language="python"> |
3929 | +import vtktools |
3930 | +v = vtktools.vtu("darcy_p0p1cv_pressBCinlet_2d_1.vtu") |
3931 | +inter_vel_2d = max(v.GetScalarRange("interstitial_velocity")) |
3932 | + </variable> |
3933 | + <variable name="pressure_3d" language="python"> |
3934 | +import vtktools |
3935 | +v = vtktools.vtu("darcy_p0p1cv_pressBCinlet_3d_1.vtu") |
3936 | +pressure_3d = v.GetScalarRange("Pressure") |
3937 | + </variable> |
3938 | + <variable name="inter_vel_3d" language="python"> |
3939 | +import vtktools |
3940 | +v = vtktools.vtu("darcy_p0p1cv_pressBCinlet_3d_1.vtu") |
3941 | +inter_vel_3d = max(v.GetScalarRange("interstitial_velocity")) |
3942 | + </variable> |
3943 | + </variables> |
3944 | + <pass_tests> |
3945 | + <test name="change_P for 1d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python"> |
3946 | +change_P = abs(max(pressure_1d) - min(pressure_1d)) |
3947 | +visc = 1.0e-04 |
3948 | +darcy_vel_BC = 5.0 |
3949 | +perm = 1.0e-10 |
3950 | +domain_length = 300.0 |
3951 | +print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) |
3952 | +assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) < 1.0e-09 |
3953 | + </test> |
3954 | + <test name="interstitial velocity 1d should equal darcy_vel/porosity, check relative difference to be under tolerance 1.0e-09" language="python"> |
3955 | +analytic_vel = 10.0 |
3956 | +print 'Relative error of interstitial velocity: ',abs((inter_vel_1d - analytic_vel)/analytic_vel) |
3957 | +assert abs((inter_vel_1d - analytic_vel)/analytic_vel) < 1.0e-09 |
3958 | + </test> |
3959 | + <test name="change_P for 2d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python"> |
3960 | +change_P = abs(max(pressure_2d) - min(pressure_2d)) |
3961 | +visc = 1.0e-04 |
3962 | +darcy_vel_BC = 5.0 |
3963 | +perm = 1.0e-10 |
3964 | +domain_length = 300.0 |
3965 | +print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) |
3966 | +assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) < 1.0e-09 |
3967 | + </test> |
3968 | + <test name="interstitial velocity 2d should equal darcy_vel/porosity, check relative difference to be under tolerance 1.0e-09" language="python"> |
3969 | +analytic_vel = 10.0 |
3970 | +print 'Relative error of interstitial velocity: ',abs((inter_vel_2d - analytic_vel)/analytic_vel) |
3971 | +assert abs((inter_vel_2d - analytic_vel)/analytic_vel) < 1.0e-09 |
3972 | + </test> |
3973 | + <test name="change_P for 3d should equal domain_length*visc*darcy_vel_BC/perm, check relative difference to be under tolerance 1.0e-09" language="python"> |
3974 | +change_P = abs(max(pressure_3d) - min(pressure_3d)) |
3975 | +visc = 1.0e-04 |
3976 | +darcy_vel_BC = 5.0 |
3977 | +perm = 1.0e-10 |
3978 | +domain_length = 300.0 |
3979 | +print 'Relative error of pressure difference: ',abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) |
3980 | +assert abs((change_P - domain_length*visc*darcy_vel_BC/perm)/(domain_length*visc*darcy_vel_BC/perm)) < 1.0e-09 |
3981 | + </test> |
3982 | + <test name="interstitial velocity 3d should equal darcy_vel/porosity, check relative difference to be under tolerance 1.0e-09" language="python"> |
3983 | +analytic_vel = 10.0 |
3984 | +print 'Relative error of interstitial velocity: ',abs((inter_vel_3d - analytic_vel)/analytic_vel) |
3985 | +assert abs((inter_vel_3d - analytic_vel)/analytic_vel) < 1.0e-09 |
3986 | + </test> |
3987 | + </pass_tests> |
3988 | + <warn_tests> |
3989 | + </warn_tests> |
3990 | +</testproblem> |
3991 | |
3992 | === added file 'tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_1d.flml' |
3993 | --- tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_1d.flml 1970-01-01 00:00:00 +0000 |
3994 | +++ tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_1d.flml 2011-12-23 12:05:25 +0000 |
3995 | @@ -0,0 +1,384 @@ |
3996 | +<?xml version='1.0' encoding='utf-8'?> |
3997 | +<fluidity_options> |
3998 | + <simulation_name> |
3999 | + <string_value lines="1">darcy_p0p1cv_pressBCinlet_1d</string_value> |
4000 | + <comment>a simple short test case for darcy flow using the element type p0p1cv for velocity-pressure |
4001 | + |
4002 | +the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution) |
4003 | + |
4004 | +it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct |
4005 | + |
4006 | +the darcy flow equation is |
4007 | + |
4008 | +sigma*darcy_vel = - grad P |
4009 | + |
4010 | +sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options) |
4011 | + |
4012 | +python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity |
4013 | + |
4014 | +this case has a 1 region 1 material with pressure boundary condition inflow. |
4015 | + |
4016 | +to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input. |
4017 | + |
4018 | +the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) |
4019 | + |
4020 | +the geometry is 1d and one time step is done (as nothing depends on time).</comment> |
4021 | + </simulation_name> |
4022 | + <problem_type> |
4023 | + <string_value lines="1">fluids</string_value> |
4024 | + </problem_type> |
4025 | + <geometry> |
4026 | + <dimension> |
4027 | + <integer_value rank="0">1</integer_value> |
4028 | + </dimension> |
4029 | + <mesh name="CoordinateMesh"> |
4030 | + <from_file file_name="darcy_p0p1cv_pressBCinlet_1d"> |
4031 | + <format name="triangle"/> |
4032 | + <stat> |
4033 | + <include_in_stat/> |
4034 | + </stat> |
4035 | + </from_file> |
4036 | + </mesh> |
4037 | + <mesh name="VelocityMesh"> |
4038 | + <from_mesh> |
4039 | + <mesh name="CoordinateMesh"/> |
4040 | + <mesh_shape> |
4041 | + <polynomial_degree> |
4042 | + <integer_value rank="0">0</integer_value> |
4043 | + </polynomial_degree> |
4044 | + </mesh_shape> |
4045 | + <mesh_continuity> |
4046 | + <string_value>discontinuous</string_value> |
4047 | + </mesh_continuity> |
4048 | + <stat> |
4049 | + <exclude_from_stat/> |
4050 | + </stat> |
4051 | + </from_mesh> |
4052 | + </mesh> |
4053 | + <mesh name="PressureMesh"> |
4054 | + <from_mesh> |
4055 | + <mesh name="CoordinateMesh"/> |
4056 | + <mesh_shape> |
4057 | + <polynomial_degree> |
4058 | + <integer_value rank="0">1</integer_value> |
4059 | + </polynomial_degree> |
4060 | + </mesh_shape> |
4061 | + <stat> |
4062 | + <exclude_from_stat/> |
4063 | + </stat> |
4064 | + </from_mesh> |
4065 | + </mesh> |
4066 | + <mesh name="DGP0"> |
4067 | + <from_mesh> |
4068 | + <mesh name="CoordinateMesh"/> |
4069 | + <mesh_shape> |
4070 | + <polynomial_degree> |
4071 | + <integer_value rank="0">0</integer_value> |
4072 | + </polynomial_degree> |
4073 | + </mesh_shape> |
4074 | + <mesh_continuity> |
4075 | + <string_value>discontinuous</string_value> |
4076 | + </mesh_continuity> |
4077 | + <stat> |
4078 | + <exclude_from_stat/> |
4079 | + </stat> |
4080 | + </from_mesh> |
4081 | + </mesh> |
4082 | + <mesh name="OutputMesh"> |
4083 | + <from_mesh> |
4084 | + <mesh name="CoordinateMesh"/> |
4085 | + <mesh_continuity> |
4086 | + <string_value>discontinuous</string_value> |
4087 | + </mesh_continuity> |
4088 | + <stat> |
4089 | + <exclude_from_stat/> |
4090 | + </stat> |
4091 | + </from_mesh> |
4092 | + </mesh> |
4093 | + <quadrature> |
4094 | + <degree> |
4095 | + <integer_value rank="0">3</integer_value> |
4096 | + </degree> |
4097 | + </quadrature> |
4098 | + </geometry> |
4099 | + <io> |
4100 | + <dump_format> |
4101 | + <string_value>vtk</string_value> |
4102 | + </dump_format> |
4103 | + <dump_period_in_timesteps> |
4104 | + <constant> |
4105 | + <integer_value rank="0">0</integer_value> |
4106 | + </constant> |
4107 | + </dump_period_in_timesteps> |
4108 | + <output_mesh name="OutputMesh"/> |
4109 | + <stat/> |
4110 | + </io> |
4111 | + <timestepping> |
4112 | + <current_time> |
4113 | + <real_value rank="0">0.0</real_value> |
4114 | + </current_time> |
4115 | + <timestep> |
4116 | + <real_value rank="0">1.0</real_value> |
4117 | + </timestep> |
4118 | + <finish_time> |
4119 | + <real_value rank="0">1.0</real_value> |
4120 | + </finish_time> |
4121 | + </timestepping> |
4122 | + <material_phase name="fluid"> |
4123 | + <scalar_field name="Pressure" rank="0"> |
4124 | + <prognostic> |
4125 | + <mesh name="PressureMesh"/> |
4126 | + <spatial_discretisation> |
4127 | + <control_volumes/> |
4128 | + </spatial_discretisation> |
4129 | + <scheme> |
4130 | + <poisson_pressure_solution> |
4131 | + <string_value lines="1">never</string_value> |
4132 | + </poisson_pressure_solution> |
4133 | + <use_projection_method/> |
4134 | + </scheme> |
4135 | + <solver> |
4136 | + <iterative_method name="cg"/> |
4137 | + <preconditioner name="sor"/> |
4138 | + <relative_error> |
4139 | + <real_value rank="0">1.0e-10</real_value> |
4140 | + </relative_error> |
4141 | + <max_iterations> |
4142 | + <integer_value rank="0">1000</integer_value> |
4143 | + </max_iterations> |
4144 | + <never_ignore_solver_failures/> |
4145 | + <diagnostics> |
4146 | + <monitors/> |
4147 | + </diagnostics> |
4148 | + </solver> |
4149 | + <boundary_conditions name="right_outflow"> |
4150 | + <surface_ids> |
4151 | + <integer_value shape="1" rank="1">2</integer_value> |
4152 | + </surface_ids> |
4153 | + <type name="dirichlet"> |
4154 | + <constant> |
4155 | + <real_value rank="0">0.0</real_value> |
4156 | + </constant> |
4157 | + </type> |
4158 | + </boundary_conditions> |
4159 | + <boundary_conditions name="left_inflow"> |
4160 | + <surface_ids> |
4161 | + <integer_value shape="1" rank="1">1</integer_value> |
4162 | + </surface_ids> |
4163 | + <type name="dirichlet"> |
4164 | + <constant> |
4165 | + <real_value rank="0">1.5e+09</real_value> |
4166 | + </constant> |
4167 | + </type> |
4168 | + </boundary_conditions> |
4169 | + <output/> |
4170 | + <stat/> |
4171 | + <convergence> |
4172 | + <include_in_convergence/> |
4173 | + </convergence> |
4174 | + <detectors> |
4175 | + <exclude_from_detectors/> |
4176 | + </detectors> |
4177 | + <steady_state> |
4178 | + <include_in_steady_state/> |
4179 | + </steady_state> |
4180 | + <no_interpolation/> |
4181 | + </prognostic> |
4182 | + </scalar_field> |
4183 | + <vector_field name="Velocity" rank="1"> |
4184 | + <prognostic> |
4185 | + <mesh name="VelocityMesh"/> |
4186 | + <equation name="Boussinesq"/> |
4187 | + <spatial_discretisation> |
4188 | + <discontinuous_galerkin> |
4189 | + <mass_terms> |
4190 | + <exclude_mass_terms/> |
4191 | + </mass_terms> |
4192 | + <viscosity_scheme> |
4193 | + <compact_discontinuous_galerkin/> |
4194 | + </viscosity_scheme> |
4195 | + <advection_scheme> |
4196 | + <none/> |
4197 | + <integrate_advection_by_parts> |
4198 | + <twice/> |
4199 | + </integrate_advection_by_parts> |
4200 | + </advection_scheme> |
4201 | + </discontinuous_galerkin> |
4202 | + <conservative_advection> |
4203 | + <real_value rank="0">0.0</real_value> |
4204 | + </conservative_advection> |
4205 | + </spatial_discretisation> |
4206 | + <temporal_discretisation> |
4207 | + <theta> |
4208 | + <real_value rank="0">1.0</real_value> |
4209 | + </theta> |
4210 | + <relaxation> |
4211 | + <real_value rank="0">1.0</real_value> |
4212 | + </relaxation> |
4213 | + </temporal_discretisation> |
4214 | + <solver> |
4215 | + <iterative_method name="gmres"> |
4216 | + <restart> |
4217 | + <integer_value rank="0">30</integer_value> |
4218 | + </restart> |
4219 | + </iterative_method> |
4220 | + <preconditioner name="sor"/> |
4221 | + <relative_error> |
4222 | + <real_value rank="0">1.0e-10</real_value> |
4223 | + </relative_error> |
4224 | + <max_iterations> |
4225 | + <integer_value rank="0">1000</integer_value> |
4226 | + </max_iterations> |
4227 | + <never_ignore_solver_failures/> |
4228 | + <diagnostics> |
4229 | + <monitors/> |
4230 | + </diagnostics> |
4231 | + </solver> |
4232 | + <initial_condition name="WholeMesh"> |
4233 | + <constant> |
4234 | + <real_value shape="1" dim1="dim" rank="1">0.0</real_value> |
4235 | + </constant> |
4236 | + </initial_condition> |
4237 | + <vector_field name="Absorption" rank="1"> |
4238 | + <diagnostic> |
4239 | + <mesh name="DGP0"/> |
4240 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
4241 | + <string_value lines="20" type="python"># get the prescribed fields |
4242 | +perm = states["fluid"].scalar_fields["Permeability"] |
4243 | +visc = states["fluid"].scalar_fields["Viscosity"] |
4244 | + |
4245 | +# loop the DMmesh_p0 nodes (which are element centred) |
4246 | +for n in range(field.node_count): |
4247 | + perm_n = perm.node_val(n) |
4248 | + visc_n = visc.node_val(n) |
4249 | + |
4250 | + # calc the absorption term |
4251 | + sigma_n = visc_n/perm_n |
4252 | + |
4253 | + # set the absorption term |
4254 | + field.set(n,sigma_n)</string_value> |
4255 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
4256 | + |
4257 | +also assume that the velocity is dg</comment> |
4258 | + </algorithm> |
4259 | + <output/> |
4260 | + <stat> |
4261 | + <include_in_stat/> |
4262 | + </stat> |
4263 | + <convergence> |
4264 | + <include_in_convergence/> |
4265 | + </convergence> |
4266 | + <detectors> |
4267 | + <include_in_detectors/> |
4268 | + </detectors> |
4269 | + <steady_state> |
4270 | + <include_in_steady_state/> |
4271 | + </steady_state> |
4272 | + </diagnostic> |
4273 | + <include_pressure_correction/> |
4274 | + </vector_field> |
4275 | + <output/> |
4276 | + <stat> |
4277 | + <include_in_stat/> |
4278 | + <previous_time_step> |
4279 | + <exclude_from_stat/> |
4280 | + </previous_time_step> |
4281 | + <nonlinear_field> |
4282 | + <exclude_from_stat/> |
4283 | + </nonlinear_field> |
4284 | + </stat> |
4285 | + <convergence> |
4286 | + <include_in_convergence/> |
4287 | + </convergence> |
4288 | + <detectors> |
4289 | + <include_in_detectors/> |
4290 | + </detectors> |
4291 | + <steady_state> |
4292 | + <include_in_steady_state/> |
4293 | + </steady_state> |
4294 | + <consistent_interpolation/> |
4295 | + </prognostic> |
4296 | + </vector_field> |
4297 | + <scalar_field name="Porosity" rank="0"> |
4298 | + <prescribed> |
4299 | + <mesh name="DGP0"/> |
4300 | + <value name="WholeMesh"> |
4301 | + <constant> |
4302 | + <real_value rank="0">0.5</real_value> |
4303 | + </constant> |
4304 | + </value> |
4305 | + <output/> |
4306 | + <stat/> |
4307 | + <detectors> |
4308 | + <exclude_from_detectors/> |
4309 | + </detectors> |
4310 | + </prescribed> |
4311 | + </scalar_field> |
4312 | + <scalar_field name="Permeability" rank="0"> |
4313 | + <prescribed> |
4314 | + <mesh name="DGP0"/> |
4315 | + <value name="WholeMesh"> |
4316 | + <constant> |
4317 | + <real_value rank="0">1.0e-10</real_value> |
4318 | + </constant> |
4319 | + </value> |
4320 | + <output/> |
4321 | + <stat/> |
4322 | + <detectors> |
4323 | + <exclude_from_detectors/> |
4324 | + </detectors> |
4325 | + </prescribed> |
4326 | + </scalar_field> |
4327 | + <scalar_field name="Viscosity" rank="0"> |
4328 | + <prescribed> |
4329 | + <mesh name="DGP0"/> |
4330 | + <value name="WholeMesh"> |
4331 | + <constant> |
4332 | + <real_value rank="0">1.0e-04</real_value> |
4333 | + </constant> |
4334 | + </value> |
4335 | + <output/> |
4336 | + <stat/> |
4337 | + <detectors> |
4338 | + <exclude_from_detectors/> |
4339 | + </detectors> |
4340 | + </prescribed> |
4341 | + </scalar_field> |
4342 | + <vector_field name="interstitial_velocity" rank="1"> |
4343 | + <diagnostic> |
4344 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
4345 | + <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"] |
4346 | +darcy_vel = states["fluid"].vector_fields["Velocity"] |
4347 | + |
4348 | +for ele in range(field.element_count): |
4349 | + phi_ele = phi.node_val(ele) |
4350 | + |
4351 | + factor_ele = 1.0/max(phi_ele,1.0e-08) |
4352 | + |
4353 | + vel_ele = [] |
4354 | + for i in range(len(darcy_vel.ele_val(ele))): |
4355 | + vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i]) |
4356 | + |
4357 | + field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value> |
4358 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
4359 | + |
4360 | +also assume that the velocity is dg</comment> |
4361 | + </algorithm> |
4362 | + <mesh name="VelocityMesh"/> |
4363 | + <output/> |
4364 | + <stat> |
4365 | + <include_in_stat/> |
4366 | + </stat> |
4367 | + <convergence> |
4368 | + <include_in_convergence/> |
4369 | + </convergence> |
4370 | + <detectors> |
4371 | + <include_in_detectors/> |
4372 | + </detectors> |
4373 | + <steady_state> |
4374 | + <include_in_steady_state/> |
4375 | + </steady_state> |
4376 | + </diagnostic> |
4377 | + </vector_field> |
4378 | + </material_phase> |
4379 | +</fluidity_options> |
4380 | |
4381 | === added file 'tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_2d.flml' |
4382 | --- tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_2d.flml 1970-01-01 00:00:00 +0000 |
4383 | +++ tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_2d.flml 2011-12-23 12:05:25 +0000 |
4384 | @@ -0,0 +1,392 @@ |
4385 | +<?xml version='1.0' encoding='utf-8'?> |
4386 | +<fluidity_options> |
4387 | + <simulation_name> |
4388 | + <string_value lines="1">darcy_p0p1cv_pressBCinlet_2d</string_value> |
4389 | + <comment>a simple short test case for darcy flow using the element type p0p1cv for velocity-pressure |
4390 | + |
4391 | +the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution) |
4392 | + |
4393 | +it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct |
4394 | + |
4395 | +the darcy flow equation is |
4396 | + |
4397 | +sigma*darcy_vel = - grad P |
4398 | + |
4399 | +sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options) |
4400 | + |
4401 | +python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity |
4402 | + |
4403 | +this case has a 1 region 1 material with pressure boundary condition inflow. |
4404 | + |
4405 | +to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input. |
4406 | + |
4407 | +the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) |
4408 | + |
4409 | +the geometry is 2d (although this is a 1d problem) and one time step is done (as nothing depends on time).</comment> |
4410 | + </simulation_name> |
4411 | + <problem_type> |
4412 | + <string_value lines="1">fluids</string_value> |
4413 | + </problem_type> |
4414 | + <geometry> |
4415 | + <dimension> |
4416 | + <integer_value rank="0">2</integer_value> |
4417 | + </dimension> |
4418 | + <mesh name="CoordinateMesh"> |
4419 | + <from_file file_name="square"> |
4420 | + <format name="gmsh"/> |
4421 | + <stat> |
4422 | + <include_in_stat/> |
4423 | + </stat> |
4424 | + </from_file> |
4425 | + </mesh> |
4426 | + <mesh name="VelocityMesh"> |
4427 | + <from_mesh> |
4428 | + <mesh name="CoordinateMesh"/> |
4429 | + <mesh_shape> |
4430 | + <polynomial_degree> |
4431 | + <integer_value rank="0">0</integer_value> |
4432 | + </polynomial_degree> |
4433 | + </mesh_shape> |
4434 | + <mesh_continuity> |
4435 | + <string_value>discontinuous</string_value> |
4436 | + </mesh_continuity> |
4437 | + <stat> |
4438 | + <exclude_from_stat/> |
4439 | + </stat> |
4440 | + </from_mesh> |
4441 | + </mesh> |
4442 | + <mesh name="PressureMesh"> |
4443 | + <from_mesh> |
4444 | + <mesh name="CoordinateMesh"/> |
4445 | + <mesh_shape> |
4446 | + <polynomial_degree> |
4447 | + <integer_value rank="0">1</integer_value> |
4448 | + </polynomial_degree> |
4449 | + </mesh_shape> |
4450 | + <stat> |
4451 | + <exclude_from_stat/> |
4452 | + </stat> |
4453 | + </from_mesh> |
4454 | + </mesh> |
4455 | + <mesh name="DGP0"> |
4456 | + <from_mesh> |
4457 | + <mesh name="CoordinateMesh"/> |
4458 | + <mesh_shape> |
4459 | + <polynomial_degree> |
4460 | + <integer_value rank="0">0</integer_value> |
4461 | + </polynomial_degree> |
4462 | + </mesh_shape> |
4463 | + <mesh_continuity> |
4464 | + <string_value>discontinuous</string_value> |
4465 | + </mesh_continuity> |
4466 | + <stat> |
4467 | + <exclude_from_stat/> |
4468 | + </stat> |
4469 | + </from_mesh> |
4470 | + </mesh> |
4471 | + <mesh name="OutputMesh"> |
4472 | + <from_mesh> |
4473 | + <mesh name="CoordinateMesh"/> |
4474 | + <mesh_continuity> |
4475 | + <string_value>discontinuous</string_value> |
4476 | + </mesh_continuity> |
4477 | + <stat> |
4478 | + <exclude_from_stat/> |
4479 | + </stat> |
4480 | + </from_mesh> |
4481 | + </mesh> |
4482 | + <quadrature> |
4483 | + <degree> |
4484 | + <integer_value rank="0">3</integer_value> |
4485 | + </degree> |
4486 | + </quadrature> |
4487 | + </geometry> |
4488 | + <io> |
4489 | + <dump_format> |
4490 | + <string_value>vtk</string_value> |
4491 | + </dump_format> |
4492 | + <dump_period_in_timesteps> |
4493 | + <constant> |
4494 | + <integer_value rank="0">0</integer_value> |
4495 | + </constant> |
4496 | + </dump_period_in_timesteps> |
4497 | + <output_mesh name="OutputMesh"/> |
4498 | + <stat/> |
4499 | + </io> |
4500 | + <timestepping> |
4501 | + <current_time> |
4502 | + <real_value rank="0">0.0</real_value> |
4503 | + </current_time> |
4504 | + <timestep> |
4505 | + <real_value rank="0">1.0</real_value> |
4506 | + </timestep> |
4507 | + <finish_time> |
4508 | + <real_value rank="0">1.0</real_value> |
4509 | + </finish_time> |
4510 | + </timestepping> |
4511 | + <material_phase name="fluid"> |
4512 | + <scalar_field name="Pressure" rank="0"> |
4513 | + <prognostic> |
4514 | + <mesh name="PressureMesh"/> |
4515 | + <spatial_discretisation> |
4516 | + <control_volumes/> |
4517 | + </spatial_discretisation> |
4518 | + <scheme> |
4519 | + <poisson_pressure_solution> |
4520 | + <string_value lines="1">never</string_value> |
4521 | + </poisson_pressure_solution> |
4522 | + <use_projection_method/> |
4523 | + </scheme> |
4524 | + <solver> |
4525 | + <iterative_method name="cg"/> |
4526 | + <preconditioner name="sor"/> |
4527 | + <relative_error> |
4528 | + <real_value rank="0">1.0e-10</real_value> |
4529 | + </relative_error> |
4530 | + <max_iterations> |
4531 | + <integer_value rank="0">1000</integer_value> |
4532 | + </max_iterations> |
4533 | + <never_ignore_solver_failures/> |
4534 | + <diagnostics> |
4535 | + <monitors/> |
4536 | + </diagnostics> |
4537 | + </solver> |
4538 | + <boundary_conditions name="right_outflow"> |
4539 | + <surface_ids> |
4540 | + <integer_value shape="1" rank="1">8</integer_value> |
4541 | + </surface_ids> |
4542 | + <type name="dirichlet"> |
4543 | + <apply_weakly/> |
4544 | + <constant> |
4545 | + <real_value rank="0">0.0</real_value> |
4546 | + </constant> |
4547 | + </type> |
4548 | + </boundary_conditions> |
4549 | + <boundary_conditions name="left_inflow"> |
4550 | + <surface_ids> |
4551 | + <integer_value shape="1" rank="1">10</integer_value> |
4552 | + </surface_ids> |
4553 | + <type name="dirichlet"> |
4554 | + <apply_weakly/> |
4555 | + <constant> |
4556 | + <real_value rank="0">1.5e+09</real_value> |
4557 | + </constant> |
4558 | + </type> |
4559 | + </boundary_conditions> |
4560 | + <output/> |
4561 | + <stat/> |
4562 | + <convergence> |
4563 | + <include_in_convergence/> |
4564 | + </convergence> |
4565 | + <detectors> |
4566 | + <exclude_from_detectors/> |
4567 | + </detectors> |
4568 | + <steady_state> |
4569 | + <include_in_steady_state/> |
4570 | + </steady_state> |
4571 | + <no_interpolation/> |
4572 | + </prognostic> |
4573 | + </scalar_field> |
4574 | + <vector_field name="Velocity" rank="1"> |
4575 | + <prognostic> |
4576 | + <mesh name="VelocityMesh"/> |
4577 | + <equation name="Boussinesq"/> |
4578 | + <spatial_discretisation> |
4579 | + <discontinuous_galerkin> |
4580 | + <mass_terms> |
4581 | + <exclude_mass_terms/> |
4582 | + </mass_terms> |
4583 | + <viscosity_scheme> |
4584 | + <compact_discontinuous_galerkin/> |
4585 | + </viscosity_scheme> |
4586 | + <advection_scheme> |
4587 | + <none/> |
4588 | + <integrate_advection_by_parts> |
4589 | + <twice/> |
4590 | + </integrate_advection_by_parts> |
4591 | + </advection_scheme> |
4592 | + </discontinuous_galerkin> |
4593 | + <conservative_advection> |
4594 | + <real_value rank="0">0.0</real_value> |
4595 | + </conservative_advection> |
4596 | + </spatial_discretisation> |
4597 | + <temporal_discretisation> |
4598 | + <theta> |
4599 | + <real_value rank="0">1.0</real_value> |
4600 | + </theta> |
4601 | + <relaxation> |
4602 | + <real_value rank="0">1.0</real_value> |
4603 | + </relaxation> |
4604 | + </temporal_discretisation> |
4605 | + <solver> |
4606 | + <iterative_method name="gmres"> |
4607 | + <restart> |
4608 | + <integer_value rank="0">30</integer_value> |
4609 | + </restart> |
4610 | + </iterative_method> |
4611 | + <preconditioner name="sor"/> |
4612 | + <relative_error> |
4613 | + <real_value rank="0">1.0e-10</real_value> |
4614 | + </relative_error> |
4615 | + <max_iterations> |
4616 | + <integer_value rank="0">1000</integer_value> |
4617 | + </max_iterations> |
4618 | + <never_ignore_solver_failures/> |
4619 | + <diagnostics> |
4620 | + <monitors/> |
4621 | + </diagnostics> |
4622 | + </solver> |
4623 | + <initial_condition name="WholeMesh"> |
4624 | + <constant> |
4625 | + <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value> |
4626 | + </constant> |
4627 | + </initial_condition> |
4628 | + <boundary_conditions name="no_outflow_top_bottom"> |
4629 | + <surface_ids> |
4630 | + <integer_value shape="2" rank="1">7 9</integer_value> |
4631 | + </surface_ids> |
4632 | + <type name="no_normal_flow"/> |
4633 | + </boundary_conditions> |
4634 | + <vector_field name="Absorption" rank="1"> |
4635 | + <diagnostic> |
4636 | + <mesh name="DGP0"/> |
4637 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
4638 | + <string_value lines="20" type="python"># get the prescribed fields |
4639 | +perm = states["fluid"].scalar_fields["Permeability"] |
4640 | +visc = states["fluid"].scalar_fields["Viscosity"] |
4641 | + |
4642 | +# loop the DMmesh_p0 nodes (which are element centred) |
4643 | +for n in range(field.node_count): |
4644 | + perm_n = perm.node_val(n) |
4645 | + visc_n = visc.node_val(n) |
4646 | + |
4647 | + # calc the absorption term |
4648 | + sigma_n = visc_n/perm_n |
4649 | + |
4650 | + # set the absorption term |
4651 | + field.set(n,sigma_n)</string_value> |
4652 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
4653 | + |
4654 | +also assume that the velocity is dg</comment> |
4655 | + </algorithm> |
4656 | + <output/> |
4657 | + <stat> |
4658 | + <include_in_stat/> |
4659 | + </stat> |
4660 | + <convergence> |
4661 | + <include_in_convergence/> |
4662 | + </convergence> |
4663 | + <detectors> |
4664 | + <include_in_detectors/> |
4665 | + </detectors> |
4666 | + <steady_state> |
4667 | + <include_in_steady_state/> |
4668 | + </steady_state> |
4669 | + </diagnostic> |
4670 | + <include_pressure_correction/> |
4671 | + </vector_field> |
4672 | + <output/> |
4673 | + <stat> |
4674 | + <include_in_stat/> |
4675 | + <previous_time_step> |
4676 | + <exclude_from_stat/> |
4677 | + </previous_time_step> |
4678 | + <nonlinear_field> |
4679 | + <exclude_from_stat/> |
4680 | + </nonlinear_field> |
4681 | + </stat> |
4682 | + <convergence> |
4683 | + <include_in_convergence/> |
4684 | + </convergence> |
4685 | + <detectors> |
4686 | + <include_in_detectors/> |
4687 | + </detectors> |
4688 | + <steady_state> |
4689 | + <include_in_steady_state/> |
4690 | + </steady_state> |
4691 | + <consistent_interpolation/> |
4692 | + </prognostic> |
4693 | + </vector_field> |
4694 | + <scalar_field name="Porosity" rank="0"> |
4695 | + <prescribed> |
4696 | + <mesh name="DGP0"/> |
4697 | + <value name="WholeMesh"> |
4698 | + <constant> |
4699 | + <real_value rank="0">0.5</real_value> |
4700 | + </constant> |
4701 | + </value> |
4702 | + <output/> |
4703 | + <stat/> |
4704 | + <detectors> |
4705 | + <exclude_from_detectors/> |
4706 | + </detectors> |
4707 | + </prescribed> |
4708 | + </scalar_field> |
4709 | + <scalar_field name="Permeability" rank="0"> |
4710 | + <prescribed> |
4711 | + <mesh name="DGP0"/> |
4712 | + <value name="WholeMesh"> |
4713 | + <constant> |
4714 | + <real_value rank="0">1.0e-10</real_value> |
4715 | + </constant> |
4716 | + </value> |
4717 | + <output/> |
4718 | + <stat/> |
4719 | + <detectors> |
4720 | + <exclude_from_detectors/> |
4721 | + </detectors> |
4722 | + </prescribed> |
4723 | + </scalar_field> |
4724 | + <scalar_field name="Viscosity" rank="0"> |
4725 | + <prescribed> |
4726 | + <mesh name="DGP0"/> |
4727 | + <value name="WholeMesh"> |
4728 | + <constant> |
4729 | + <real_value rank="0">1.0e-04</real_value> |
4730 | + </constant> |
4731 | + </value> |
4732 | + <output/> |
4733 | + <stat/> |
4734 | + <detectors> |
4735 | + <exclude_from_detectors/> |
4736 | + </detectors> |
4737 | + </prescribed> |
4738 | + </scalar_field> |
4739 | + <vector_field name="interstitial_velocity" rank="1"> |
4740 | + <diagnostic> |
4741 | + <algorithm name="vector_python_diagnostic" material_phase_support="multiple"> |
4742 | + <string_value lines="20" type="python">phi = states["fluid"].scalar_fields["Porosity"] |
4743 | +darcy_vel = states["fluid"].vector_fields["Velocity"] |
4744 | + |
4745 | +for ele in range(field.element_count): |
4746 | + phi_ele = phi.node_val(ele) |
4747 | + |
4748 | + factor_ele = 1.0/max(phi_ele,1.0e-08) |
4749 | + |
4750 | + vel_ele = [] |
4751 | + for i in range(len(darcy_vel.ele_val(ele))): |
4752 | + vel_ele.append(factor_ele*darcy_vel.ele_val(ele)[i]) |
4753 | + |
4754 | + field.set(darcy_vel.ele_nodes(ele),vel_ele)</string_value> |
4755 | + <comment>this assumes that the porosity, absorption, permeability, and viscosity are defined on a p0 dg mesh (ie. element wise) such that the number of nodes is also the number of elements |
4756 | + |
4757 | +also assume that the velocity is dg</comment> |
4758 | + </algorithm> |
4759 | + <mesh name="VelocityMesh"/> |
4760 | + <output/> |
4761 | + <stat> |
4762 | + <include_in_stat/> |
4763 | + </stat> |
4764 | + <convergence> |
4765 | + <include_in_convergence/> |
4766 | + </convergence> |
4767 | + <detectors> |
4768 | + <include_in_detectors/> |
4769 | + </detectors> |
4770 | + <steady_state> |
4771 | + <include_in_steady_state/> |
4772 | + </steady_state> |
4773 | + </diagnostic> |
4774 | + </vector_field> |
4775 | + </material_phase> |
4776 | +</fluidity_options> |
4777 | |
4778 | === added file 'tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_3d.flml' |
4779 | --- tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_3d.flml 1970-01-01 00:00:00 +0000 |
4780 | +++ tests/darcy_p0p1cv_pressBCinlet/darcy_p0p1cv_pressBCinlet_3d.flml 2011-12-23 12:05:25 +0000 |
4781 | @@ -0,0 +1,398 @@ |
4782 | +<?xml version='1.0' encoding='utf-8'?> |
4783 | +<fluidity_options> |
4784 | + <simulation_name> |
4785 | + <string_value lines="1">darcy_p0p1cv_pressBCinlet_3d</string_value> |
4786 | + <comment>a simple short test case for darcy flow using the element type p0p1cv for velocity-pressure |
4787 | + |
4788 | +the darcy velocity is solved for and the interstitial velocity is a diagnostic from this (assuming a dg solution) |
4789 | + |
4790 | +it compares the pressure gradient against the analytic, as well as checking that the interstitial velocity is correct |
4791 | + |
4792 | +the darcy flow equation is |
4793 | + |
4794 | +sigma*darcy_vel = - grad P |
4795 | + |
4796 | +sigma is a function of viscosity and permeability which are two input fields that are named in the flml - not on fluidity special names (ie. this doesnt use the porous media object options) |
4797 | + |
4798 | +python diagnostics are used to form the sigma term that is assoicated with a p0 dg mesh (ie. element wise) and the interstitial velocity |
4799 | + |
4800 | +this case has a 1 region 1 material with pressure boundary condition inflow. |
4801 | + |
4802 | +to get a darcy equation the time, stress and advection terms are not included in a bousinessq formulation (to avoid needing to specify a density). Ideally it would be better to actually have a darcy momentum option to simplify the input. |
4803 | + |
4804 | +the absorption term is included in the pressure correction (being a necessity as it is the only term in the matrix equation) |
4805 | + |
4806 | +the geometry is 3d (although this is a 1d problem) and one time step is done (as nothing depends on time).</comment> |
4807 | + </simulation_name> |
4808 | + <problem_type> |
4809 | + <string_value lines="1">fluids</string_value> |
4810 | + </problem_type> |
4811 | + <geometry> |
4812 | + <dimension> |
4813 | + <integer_value rank="0">3</integer_value> |
4814 | + </dimension> |
4815 | + <mesh name="CoordinateMesh"> |
4816 | + <from_file file_name="cube"> |
4817 | + <format name="gmsh"/> |
4818 | + <stat> |
4819 | + <include_in_stat/> |
4820 | + </stat> |
4821 | + </from_file> |
4822 | + </mesh> |
4823 | + <mesh name="VelocityMesh"> |
4824 | + <from_mesh> |
4825 | + <mesh name="CoordinateMesh"/> |
4826 | + <mesh_shape> |
4827 | + <polynomial_degree> |
4828 | + <integer_value rank="0">0</integer_value> |
4829 | + </polynomial_degree> |
4830 | + </mesh_shape> |
4831 | + <mesh_continuity> |
4832 | + <string_value>discontinuous</string_value> |
4833 | + </mesh_continuity> |
4834 | + <stat> |
4835 | + <exclude_from_stat/> |
4836 | + </stat> |
4837 | + </from_mesh> |
4838 | + </mesh> |
4839 | + <mesh name="PressureMesh"> |
4840 | + <from_mesh> |
4841 | + <mesh name="CoordinateMesh"/> |
4842 | + <mesh_shape> |
4843 | + <polynomial_degree> |
4844 | + <integer_value rank="0">1</integer_value> |
4845 | + </polynomial_degree> |
4846 | + </mesh_shape> |
4847 | + <stat> |
4848 | + <exclude_from_stat/> |
4849 | + </stat> |
4850 | + </from_mesh> |
4851 | + </mesh> |
4852 | + <mesh name="DGP0"> |
4853 | + <from_mesh> |
4854 | + <mesh name="CoordinateMesh"/> |
4855 | + <mesh_shape> |
4856 | + <polynomial_degree> |
4857 | + <integer_value rank="0">0</integer_value> |
4858 | + </polynomial_degree> |
4859 | + </mesh_shape> |
4860 | + <mesh_continuity> |
4861 | + <string_value>discontinuous</string_value> |
4862 | + </mesh_continuity> |
4863 | + <stat> |
4864 | + <exclude_from_stat/> |
4865 | + </stat> |
4866 | + </from_mesh> |
4867 | + </mesh> |
4868 | + <mesh name="OutputMesh"> |
4869 | + <from_mesh> |
4870 | + <mesh name="CoordinateMesh"/> |
4871 | + <mesh_continuity> |
4872 | + <string_value>discontinuous</string_value> |
4873 | + </mesh_continuity> |
4874 | + <stat> |
4875 | + <exclude_from_stat/> |
4876 | + </stat> |
4877 | + </from_mesh> |
4878 | + </mesh> |
4879 | + <quadrature> |
4880 | + <degree> |
4881 | + <integer_value rank="0">3</integer_value> |
4882 | + </degree> |
4883 | + </quadrature> |
4884 | + </geometry> |
4885 | + <io> |
4886 | + <dump_format> |
4887 | + <string_value>vtk</string_value> |
4888 | + </dump_format> |
4889 | + <dump_period_in_timesteps> |
4890 | + <constant> |
4891 | + <integer_value rank="0">0</integer_value> |
4892 | + </constant> |
4893 | + </dump_period_in_timesteps> |
4894 | + <output_mesh name="OutputMesh"/> |
4895 | + <stat/> |
4896 | + </io> |
4897 | + <timestepping> |
4898 | + <current_time> |
4899 | + <real_value rank="0">0.0</real_value> |
4900 | + </current_time> |
4901 | + <timestep> |
4902 | + <real_value rank="0">1.0</real_value> |
4903 | + </timestep> |
4904 | + <finish_time> |
4905 | + <real_value rank="0">1.0</real_value> |
4906 | + </finish_time> |
4907 | + </timestepping> |
4908 | + <material_phase name="fluid"> |
4909 | + <scalar_field name="Pressure" rank="0"> |
4910 | + <prognostic> |
4911 | + <mesh name="PressureMesh"/> |
4912 | + <spatial_discretisation> |
4913 | + <control_volumes/> |
4914 | + </spatial_discretisation> |
4915 | + <scheme> |
4916 | + <poisson_pressure_solution> |
4917 | + <string_value lines="1">never</string_value> |
4918 | + </poisson_pressure_solution> |
4919 | + <use_projection_method/> |
4920 | + </scheme> |
4921 | + <solver> |
4922 | + <iterative_method name="cg"/> |
4923 | + <preconditioner name="sor"/> |
4924 | + <relative_error> |
4925 | + <real_value rank="0">1.0e-10</real_value> |
4926 | + </relative_error> |
4927 | + <max_iterations> |
4928 | + <integer_value rank="0">1000</integer_value> |
4929 | + </max_iterations> |
4930 | + <never_ignore_solver_failures/> |
4931 | + <diagnostics> |
4932 | + <monitors/> |
4933 | + </diagnostics> |
4934 | + </solver> |
4935 | + <boundary_conditions name="right_outflow"> |
4936 | + <surface_ids> |
4937 | + <integer_value shape="1" rank="1">4</integer_value> |
4938 | + </surface_ids> |
4939 | + <type name="dirichlet"> |
4940 | + <apply_weakly/> |
4941 | + <constant> |
4942 | + <real_value rank="0">0.0</real_value> |
4943 | + </constant> |
4944 | + </type> |
4945 | + </boundary_conditions> |
4946 | + <boundary_conditions name="left_inflow"> |
4947 | + <surface_ids> |
4948 | + <integer_value shape="1" rank="1">3</integer_value> |
4949 | + </surface_ids> |
4950 | + <type name="dirichlet"> |
4951 | + <apply_weakly/> |
4952 | + <constant> |
4953 | + <real_value rank="0">1.5e+09</real_value> |
4954 | + </constant> |
4955 | + </type> |
4956 | + </boundary_conditions> |
4957 | + <output/> |
4958 | + <stat/> |
4959 | + <convergence> |
4960 | + <include_in_convergence/> |
4961 | + </convergence> |
4962 | + <detectors> |
4963 | + <exclude_from_detectors/> |
4964 | + </detectors> |
4965 | + <steady_state> |
4966 | + <include_in_steady_state/> |
4967 | + </steady_state> |
4968 | + <no_interpolation/> |
4969 | + </prognostic> |
4970 | + </scalar_field> |
4971 | + <vector_field name="Velocity" rank="1"> |
4972 | + <prognostic> |
4973 | + <mesh name="VelocityMesh"/> |
4974 | + <equation name="Boussinesq"/> |
4975 | + <spatial_discretisation> |
4976 | + <discontinuous_galerkin> |
4977 | + <mass_terms> |
4978 | + <exclude_mass_terms/> |
4979 | + </mass_terms> |
4980 | + <viscosity_scheme> |
4981 | + <compact_discontinuous_galerkin/> |
4982 | + </viscosity_scheme> |
4983 | + <advection_scheme> |
4984 | + <none/> |
4985 | + <integrate_advection_by_parts> |
4986 | + <twice/> |
4987 | + </integrate_advection_by_parts> |
4988 | + </advection_scheme> |
4989 | + </discontinuous_galerkin> |
4990 | + <conservative_advection> |
4991 | + <real_value rank="0">0.0</real_value> |
4992 | + </conservative_advection> |
4993 | + </spatial_discretisation> |
4994 | + <temporal_discretisation> |
4995 | + <theta> |
4996 | + <real_value rank="0">1.0</real_value> |
4997 | + </theta> |
4998 | + <relaxation> |
4999 | + <real_value rank="0">1.0</real_value> |
5000 | + </relaxation> |
Since original merge request, added a very small change to a comment.
Merged latest trunk into branch and re-ran branch buildbot which is green.
Needs a volunteer to review please!