Merge lp:~ctjacobs-multiphase/fluidity/compressible-multiphase into lp:fluidity
- compressible-multiphase
- Merge into dev-trunk
Status: | Merged |
---|---|
Approved by: | Christian Jacobs |
Approved revision: | 4067 |
Merged at revision: | 4024 |
Proposed branch: | lp:~ctjacobs-multiphase/fluidity/compressible-multiphase |
Merge into: | lp:fluidity |
Diff against target: |
46667 lines (+44935/-226) 90 files modified
assemble/Advection_Diffusion_CG.F90 (+159/-38) assemble/Compressible_Projection.F90 (+81/-28) assemble/Diagnostic_Fields_Matrices.F90 (+161/-13) assemble/Diagnostic_fields_wrapper.F90 (+15/-0) assemble/Divergence_Matrix_CG.F90 (+146/-33) assemble/Makefile.dependencies (+10/-10) assemble/Momentum_Equation.F90 (+37/-27) assemble/Multiphase.F90 (+391/-40) main/Fluids.F90 (+1/-1) main/Shallow_Water.F90 (+17/-16) manual/bibliography.bib (+29/-0) manual/configuring_fluidity.tex (+52/-13) manual/model_equations.tex (+21/-5) manual/visualisation_and_diagnostics.tex (+1/-0) schemas/fluidity_options.rnc (+35/-2) schemas/fluidity_options.rng (+40/-0) schemas/multiphase_interaction.rnc (+58/-0) schemas/multiphase_interaction.rng (+65/-0) tests/flux_bc/Makefile (+22/-0) tests/flux_bc/flux_bc.xml (+44/-0) tests/flux_bc/flux_bc_2d.flml (+396/-0) tests/flux_bc/flux_bc_3d.flml (+405/-0) tests/flux_bc/src/flux_bc_2d.geo (+19/-0) tests/flux_bc/src/flux_bc_3d.geo (+69/-0) tests/mphase_dusty_gas_shock_tube/Makefile (+18/-0) tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml (+733/-0) tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.xml (+156/-0) tests/mphase_dusty_gas_shock_tube/plot.py (+93/-0) tests/mphase_dusty_gas_shock_tube/shocktube.py (+157/-0) tests/mphase_dusty_gas_shock_tube/single_phase_frozen_flow_test.flml (+379/-0) tests/mphase_dusty_gas_shock_tube/single_phase_frozen_flow_test.py (+75/-0) tests/mphase_inlet_velocity_bc_incompressible/Makefile (+17/-0) tests/mphase_inlet_velocity_bc_incompressible/mphase_inlet_velocity_bc_incompressible.flml (+527/-0) tests/mphase_inlet_velocity_bc_incompressible/mphase_inlet_velocity_bc_incompressible.xml (+63/-0) tests/mphase_mms_p2p1_compressible/MMS_A.flml (+933/-0) tests/mphase_mms_p2p1_compressible/MMS_B.flml (+933/-0) tests/mphase_mms_p2p1_compressible/MMS_C.flml (+933/-0) tests/mphase_mms_p2p1_compressible/MMS_D.flml (+933/-0) tests/mphase_mms_p2p1_compressible/Makefile (+22/-0) tests/mphase_mms_p2p1_compressible/get_residual_compressible.sage (+74/-0) tests/mphase_mms_p2p1_compressible/mphase_mms_p2p1_compressible.xml (+349/-0) tests/mphase_mms_p2p1_compressible/src/MMS_A.geo (+15/-0) tests/mphase_mms_p2p1_compressible/src/MMS_A.msh (+190/-0) tests/mphase_mms_p2p1_compressible/src/MMS_B.geo (+15/-0) tests/mphase_mms_p2p1_compressible/src/MMS_B.msh (+739/-0) tests/mphase_mms_p2p1_compressible/src/MMS_C.geo (+15/-0) tests/mphase_mms_p2p1_compressible/src/MMS_C.msh (+3013/-0) tests/mphase_mms_p2p1_compressible/src/MMS_D.geo (+15/-0) tests/mphase_mms_p2p1_compressible/src/MMS_D.msh (+12055/-0) tests/mphase_mms_p2p1_compressible_ie/MMS_A.flml (+1233/-0) tests/mphase_mms_p2p1_compressible_ie/MMS_B.flml (+1233/-0) tests/mphase_mms_p2p1_compressible_ie/MMS_C.flml (+1233/-0) tests/mphase_mms_p2p1_compressible_ie/Makefile (+20/-0) tests/mphase_mms_p2p1_compressible_ie/get_residual_compressible.sage (+75/-0) tests/mphase_mms_p2p1_compressible_ie/mphase_mms_p2p1_compressible_ie.xml (+311/-0) tests/mphase_mms_p2p1_compressible_ie/src/MMS_A.geo (+15/-0) tests/mphase_mms_p2p1_compressible_ie/src/MMS_A.msh (+238/-0) tests/mphase_mms_p2p1_compressible_ie/src/MMS_B.geo (+15/-0) tests/mphase_mms_p2p1_compressible_ie/src/MMS_B.msh (+832/-0) tests/mphase_mms_p2p1_compressible_ie/src/MMS_C.geo (+15/-0) tests/mphase_mms_p2p1_compressible_ie/src/MMS_C.msh (+3238/-0) tests/mphase_mms_p2p1_compressible_ie_heat_transfer/MMS_A.flml (+1316/-0) tests/mphase_mms_p2p1_compressible_ie_heat_transfer/MMS_B.flml (+1316/-0) tests/mphase_mms_p2p1_compressible_ie_heat_transfer/MMS_C.flml (+1316/-0) tests/mphase_mms_p2p1_compressible_ie_heat_transfer/Makefile (+20/-0) tests/mphase_mms_p2p1_compressible_ie_heat_transfer/get_residual_compressible.sage (+75/-0) tests/mphase_mms_p2p1_compressible_ie_heat_transfer/mphase_mms_p2p1_compressible_ie_heat_transfer.xml (+311/-0) tests/mphase_mms_p2p1_compressible_ie_heat_transfer/src/MMS_A.geo (+15/-0) tests/mphase_mms_p2p1_compressible_ie_heat_transfer/src/MMS_A.msh (+238/-0) tests/mphase_mms_p2p1_compressible_ie_heat_transfer/src/MMS_B.geo (+15/-0) tests/mphase_mms_p2p1_compressible_ie_heat_transfer/src/MMS_B.msh (+832/-0) tests/mphase_mms_p2p1_compressible_ie_heat_transfer/src/MMS_C.geo (+15/-0) tests/mphase_mms_p2p1_compressible_ie_heat_transfer/src/MMS_C.msh (+3238/-0) tests/mphase_rogue_shock_tube_dense_bed_glass/Makefile (+18/-0) tests/mphase_rogue_shock_tube_dense_bed_glass/mphase_rogue_shock_tube_dense_bed_glass.flml (+824/-0) tests/mphase_rogue_shock_tube_dense_bed_glass/mphase_rogue_shock_tube_dense_bed_glass.xml (+108/-0) tests/mphase_rogue_shock_tube_dense_bed_glass/plot.py (+70/-0) tests/mphase_rogue_shock_tube_dense_bed_glass/src/mphase_rogue_shock_tube_dense_bed_glass.geo (+19/-0) tests/mphase_rogue_shock_tube_dense_bed_nylon/Makefile (+17/-0) tests/mphase_rogue_shock_tube_dense_bed_nylon/mphase_rogue_shock_tube_dense_bed_nylon.flml (+824/-0) tests/mphase_rogue_shock_tube_dense_bed_nylon/mphase_rogue_shock_tube_dense_bed_nylon.xml (+108/-0) tests/mphase_rogue_shock_tube_dense_bed_nylon/plot.py (+96/-0) tests/mphase_rogue_shock_tube_dense_bed_nylon/src/mphase_rogue_shock_tube_dense_bed_nylon.geo (+19/-0) tests/mphase_stokes_law/mphase_stokes_law.flml (+5/-0) tests/mphase_tephra_settling/mphase_tephra_settling.flml (+5/-0) tests/mphase_wen_yu_drag_correlation/Makefile (+17/-0) tests/mphase_wen_yu_drag_correlation/c_against_re.py (+26/-0) tests/mphase_wen_yu_drag_correlation/mphase_wen_yu_drag_correlation.flml (+793/-0) tests/mphase_wen_yu_drag_correlation/mphase_wen_yu_drag_correlation.xml (+76/-0) tests/mphase_wen_yu_drag_correlation/src/mphase_wen_yu_drag_correlation.geo (+19/-0) |
To merge this branch: | bzr merge lp:~ctjacobs-multiphase/fluidity/compressible-multiphase |
Related bugs: |
Reviewer | Review Type | Date Requested | Status |
---|---|---|---|
Stephan Kramer | Approve | ||
Cian Wilson | Pending | ||
Review via email: mp+114294@code.launchpad.net |
Commit message
Description of the change
Proposing merge of compressible-
Summary of changes:
===================
- Added multiphase support for the compressible projection method (currently only for CG Pressure).
- Small changes to the logic in Momentum_
- Added multiphase support for the InternalEnergy equation. Note that the whole state array is now passed to solve_field_
- Implemented the heat transfer term by Gunn (1978).
- Added three P2-P1 multiphase MMS tests for model verification: one without InternalEnergy fields, one with InternalEnergy fields, and one with InternalEnergy fields and heat transfer. The tests will become longtests after merging into the trunk.
- The manual now has documentation on the compressible multiphase model, including the multiphase version of the InternalEnergy equation.
- Added a diagnostic field called CompressibleCon
- Added two gas-solid shock tube tests to help to validate the model: mphase_
- Added another gas-solid shock tube test (mphase_
- Implemented two new correlations for fluid-particle drag. The one by Wen & Yu (1966) can support higher Reynolds numbers (compared to the existing Stokes' law correlation). The one by Ergun (1952) can support dense flows where the volume fraction of the particle phase is greater than 0.2, like in the dense bed shock tube tests (and in the pyroclastic flow tests which will hopefully come later).
- Changed the structure of the schema. Users now explicitly select which interaction terms they wish to include, as well as the drag correlation, under "/multiphase_
- Added a test for the Wen & Yu (1966) drag correlation. This is similar to mphase_stokes_law.
- Added a test to check if inlet velocity boundary conditions work with the incompressible multiphase flow model.
- Added a test for the flux boundary condition.
Stephan Kramer (s-kramer) wrote : | # |
- 4051. By Christian Jacobs
-
Improvements to Advection_
Diffusion_ CG.F90 based on Stephan's feedback. - 4052. By Christian Jacobs
-
Moving the call to add_heat_transfer up one level.
- 4053. By Christian Jacobs
-
Merging trunk revisions (from r4011 to r4019 inclusive) into branch.
Christian Jacobs (ctjacobs) wrote : | # |
Thanks for the feedback Stephan - improvements have now been made to Advection_
Stephan Kramer (s-kramer) wrote : | # |
Some more comments trickling through:
* assemble/
* assemble/
multiphase variable there, it's not used afaics and a bit confusing...
* assemble/
but it doesn't seem to be used? Also in the same routine you've added vfrac
and temp_vfrac that aren't used.
* calculate_
better maybe be compressible_
what "compressible_
- 4054. By Christian Jacobs
-
Some more improvements based on Stephan's feedback.
- 4055. By Christian Jacobs
-
Added a test for the CompressibleCon
tinuityResidual diagnostic field. - 4056. By Christian Jacobs
-
More tests for the CompressibleCon
tinuityResidual diagnostic field.
Christian Jacobs (ctjacobs) wrote : | # |
Thanks again - improvements made and the (renamed) CompressibleCon
- 4057. By Christian Jacobs
-
Merging trunk revisions (from r4020 to r4022 inclusive) into branch.
Stephan Kramer (s-kramer) wrote : | # |
Some final, small remarks:
* assemble/
* assemble/
* assemble/
a Density, but do have a prognostic Velocity and Pressure?
* assemble/
sensible, but I'm not sure and I don't like second-guessing. A solution would be to move the select case up one level, and set an integer with some "enum"-like parameters DRAG_CORRELATIO
* assemble/
* assemble/
Great stuff, well documented, lots of tests. Next time, could we maybe merge this in smaller bits?
- 4058. By Christian Jacobs
-
Using state%option_path instead.
- 4059. By Christian Jacobs
-
- Don't try to extract the Density field here, as it's not necessary.
- Another change to use state%option_path instead. - 4060. By Christian Jacobs
-
Call FLExit if use_theta_pg or use_theta_
divergence are true when running a compressible multiphase flow simulation. - 4061. By Christian Jacobs
-
Introduced integer parameters to represent the different drag correlations.
- 4062. By Christian Jacobs
-
Merging trunk revision r4023 into branch.
- 4063. By Christian Jacobs
-
Corrected a typo introduced when changing to state%option_path.
- 4064. By Christian Jacobs
-
- Corrected a few typos in the XML files.
- Set the mphase_rogue_shock_ tube_dense_ bed_glass test's length to "medium". - 4065. By Christian Jacobs
-
Line continuation fix. The code now compiles on the Intel compiler.
- 4066. By Christian Jacobs
-
Fixed an assert statement. Only noticed this when debugging was enabled.
- 4067. By Christian Jacobs
-
Updates from "make makefiles".
Christian Jacobs (ctjacobs) wrote : | # |
The build queue now has a green light on buildbot. Thanks for the review.
Preview Diff
1 | === modified file 'assemble/Advection_Diffusion_CG.F90' |
2 | --- assemble/Advection_Diffusion_CG.F90 2012-04-16 16:41:59 +0000 |
3 | +++ assemble/Advection_Diffusion_CG.F90 2012-08-09 00:46:20 +0000 |
4 | @@ -47,6 +47,7 @@ |
5 | use upwind_stabilisation |
6 | use sparsity_patterns_meshes |
7 | use porous_media |
8 | + use multiphase_module |
9 | use sparse_tools_petsc |
10 | use colouring |
11 | #ifdef _OPENMP |
12 | @@ -117,16 +118,19 @@ |
13 | logical :: move_mesh |
14 | ! Include porosity? |
15 | logical :: include_porosity |
16 | + ! Are we running a multiphase flow simulation? |
17 | + logical :: multiphase |
18 | |
19 | contains |
20 | |
21 | - subroutine solve_field_equation_cg(field_name, state, dt, velocity_name, iterations_taken) |
22 | + subroutine solve_field_equation_cg(field_name, state, istate, dt, velocity_name, iterations_taken) |
23 | !!< Construct and solve the advection-diffusion equation for the given |
24 | !!< field using a continuous Galerkin discretisation. Based on |
25 | !!< Advection_Diffusion_DG and Momentum_CG. |
26 | |
27 | character(len = *), intent(in) :: field_name |
28 | - type(state_type), intent(inout) :: state |
29 | + type(state_type), dimension(:), intent(inout) :: state |
30 | + integer, intent(in) :: istate |
31 | real, intent(in) :: dt |
32 | character(len = *), optional, intent(in) :: velocity_name |
33 | integer, intent(out), optional :: iterations_taken |
34 | @@ -138,16 +142,22 @@ |
35 | ewrite(1, *) "In solve_field_equation_cg" |
36 | |
37 | ewrite(2, *) "Solving advection-diffusion equation for field " // & |
38 | - & trim(field_name) // " in state " // trim(state%name) |
39 | + & trim(field_name) // " in state " // trim(state(istate)%name) |
40 | |
41 | - call initialise_advection_diffusion_cg(field_name, t, delta_t, matrix, rhs, state) |
42 | + call initialise_advection_diffusion_cg(field_name, t, delta_t, matrix, rhs, state(istate)) |
43 | |
44 | call profiler_tic(t, "assembly") |
45 | - call assemble_advection_diffusion_cg(t, matrix, rhs, state, dt, velocity_name = velocity_name) |
46 | + call assemble_advection_diffusion_cg(t, matrix, rhs, state(istate), dt, velocity_name = velocity_name) |
47 | + |
48 | + ! Note: the assembly of the heat transfer term is done here to avoid |
49 | + ! passing in the whole state array to assemble_advection_diffusion_cg. |
50 | + if(have_option("/multiphase_interaction/heat_transfer")) then |
51 | + call add_heat_transfer(state, istate, t, matrix, rhs) |
52 | + end if |
53 | call profiler_toc(t, "assembly") |
54 | - |
55 | + |
56 | call profiler_tic(t, "solve_total") |
57 | - call solve_advection_diffusion_cg(t, delta_t, matrix, rhs, state, & |
58 | + call solve_advection_diffusion_cg(t, delta_t, matrix, rhs, state(istate), & |
59 | iterations_taken = iterations_taken) |
60 | call profiler_toc(t, "solve_total") |
61 | |
62 | @@ -234,6 +244,11 @@ |
63 | type(scalar_field), pointer :: density, olddensity |
64 | character(len = FIELD_NAME_LEN) :: density_name |
65 | type(scalar_field), pointer :: pressure |
66 | + |
67 | + ! Volume fraction fields for multiphase flow simulation |
68 | + type(scalar_field), pointer :: vfrac |
69 | + type(scalar_field) :: nvfrac ! Non-linear version |
70 | + |
71 | ! Porosity field |
72 | type(scalar_field) :: porosity_theta |
73 | |
74 | @@ -373,6 +388,7 @@ |
75 | call allocate(porosity_theta, t%mesh, field_type=FIELD_TYPE_CONSTANT) |
76 | call set(porosity_theta, 1.0) |
77 | end if |
78 | + |
79 | |
80 | ! Step 2: Pull options out of the options tree |
81 | |
82 | @@ -462,6 +478,21 @@ |
83 | stabilisation_scheme = STABILISATION_NONE |
84 | end if |
85 | |
86 | + ! PhaseVolumeFraction for multiphase flow simulations |
87 | + if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then |
88 | + multiphase = .true. |
89 | + vfrac => extract_scalar_field(state, "PhaseVolumeFraction") |
90 | + call allocate(nvfrac, vfrac%mesh, "NonlinearPhaseVolumeFraction") |
91 | + call zero(nvfrac) |
92 | + call get_nonlinear_volume_fraction(state, nvfrac) |
93 | + |
94 | + ewrite_minmax(nvfrac) |
95 | + else |
96 | + multiphase = .false. |
97 | + call allocate(nvfrac, t%mesh, "DummyNonlinearPhaseVolumeFraction", field_type=FIELD_TYPE_CONSTANT) |
98 | + call set(nvfrac, 1.0) |
99 | + end if |
100 | + |
101 | call allocate(dummydensity, t%mesh, "DummyDensity", field_type=FIELD_TYPE_CONSTANT) |
102 | call set(dummydensity, 1.0) |
103 | ! find out equation type and hence if density is needed or not |
104 | @@ -479,19 +510,28 @@ |
105 | if(move_mesh) then |
106 | FLExit("Haven't implemented a moving mesh energy equation yet.") |
107 | end if |
108 | + |
109 | + ! Get old and current densities |
110 | call get_option(trim(t%option_path)//'/prognostic/equation[0]/density[0]/name', & |
111 | density_name) |
112 | density=>extract_scalar_field(state, trim(density_name)) |
113 | ewrite_minmax(density) |
114 | - |
115 | olddensity=>extract_scalar_field(state, "Old"//trim(density_name)) |
116 | ewrite_minmax(olddensity) |
117 | |
118 | - call get_option(trim(density%option_path)//"/prognostic/temporal_discretisation/theta", & |
119 | - density_theta) |
120 | + if(.not.multiphase .or. have_option('/material_phase::'//trim(state%name)//& |
121 | + &'/equation_of_state/compressible')) then |
122 | + call get_option(trim(density%option_path)//"/prognostic/temporal_discretisation/theta", & |
123 | + density_theta) |
124 | + else |
125 | + ! Since the particle phase is always incompressible then its Density |
126 | + ! will not be prognostic. Just use a fixed theta value of 1.0. |
127 | + density_theta = 1.0 |
128 | + end if |
129 | |
130 | pressure=>extract_scalar_field(state, "Pressure") |
131 | ewrite_minmax(pressure) |
132 | + |
133 | case default |
134 | FLExit("Unknown field equation type for cg advection diffusion.") |
135 | end select |
136 | @@ -531,7 +571,7 @@ |
137 | positions, old_positions, new_positions, & |
138 | velocity, grid_velocity, & |
139 | source, absorption, diffusivity, & |
140 | - density, olddensity, pressure, porosity_theta, & |
141 | + density, olddensity, pressure, porosity_theta, nvfrac, & |
142 | supg_element(thread_num+1)) |
143 | end do element_loop |
144 | !$OMP END DO |
145 | @@ -545,6 +585,7 @@ |
146 | ! which must be included before dirichlet BC's. |
147 | if (add_src_directly_to_rhs) call addto(rhs, source) |
148 | |
149 | + |
150 | ! Step 4: Boundary conditions |
151 | |
152 | if( & |
153 | @@ -571,7 +612,7 @@ |
154 | call assemble_advection_diffusion_face_cg(i, t_bc_types(i), t, t_bc, t_bc_2, & |
155 | matrix, rhs, & |
156 | positions, velocity, grid_velocity, & |
157 | - density, olddensity) |
158 | + density, olddensity, nvfrac) |
159 | end do |
160 | |
161 | call deallocate(t_bc) |
162 | @@ -586,6 +627,7 @@ |
163 | ewrite_minmax(rhs) |
164 | |
165 | call deallocate(velocity) |
166 | + call deallocate(nvfrac) |
167 | call deallocate(dummydensity) |
168 | if (stabilisation_scheme == STABILISATION_SUPG) then |
169 | do i = 1, num_threads |
170 | @@ -645,7 +687,7 @@ |
171 | positions, old_positions, new_positions, & |
172 | velocity, grid_velocity, & |
173 | source, absorption, diffusivity, & |
174 | - density, olddensity, pressure, porosity_theta, supg_shape) |
175 | + density, olddensity, pressure, porosity_theta, nvfrac, supg_shape) |
176 | integer, intent(in) :: ele |
177 | type(scalar_field), intent(in) :: t |
178 | type(csr_matrix), intent(inout) :: matrix |
179 | @@ -661,6 +703,7 @@ |
180 | type(scalar_field), intent(in) :: olddensity |
181 | type(scalar_field), intent(in) :: pressure |
182 | type(scalar_field), intent(in) :: porosity_theta |
183 | + type(scalar_field), intent(in) :: nvfrac |
184 | type(element_type), intent(inout) :: supg_shape |
185 | |
186 | integer, dimension(:), pointer :: element_nodes |
187 | @@ -669,6 +712,9 @@ |
188 | real, dimension(ele_loc(density, ele), ele_ngi(density, ele), mesh_dim(density)) :: drho_t |
189 | real, dimension(ele_loc(velocity, ele), ele_ngi(velocity, ele), mesh_dim(t)) :: du_t |
190 | real, dimension(ele_loc(positions, ele), ele_ngi(velocity, ele), mesh_dim(t)) :: dug_t |
191 | + ! Derivative of shape function for nvfrac field |
192 | + real, dimension(ele_loc(nvfrac, ele), ele_ngi(nvfrac, ele), mesh_dim(nvfrac)) :: dnvfrac_t |
193 | + |
194 | real, dimension(mesh_dim(t), mesh_dim(t), ele_ngi(t, ele)) :: j_mat |
195 | type(element_type) :: test_function |
196 | type(element_type), pointer :: t_shape |
197 | @@ -741,6 +787,16 @@ |
198 | end if |
199 | end if |
200 | |
201 | + if(have_advection .and. multiphase .and. (equation_type==FIELD_EQUATION_INTERNALENERGY)) then |
202 | + ! If the field and nvfrac meshes are different, then we need to |
203 | + ! compute the derivatives of the nvfrac shape functions. |
204 | + if(ele_shape(nvfrac, ele) == t_shape) then |
205 | + dnvfrac_t = dt_t |
206 | + else |
207 | + call transform_to_physical(positions, ele, ele_shape(nvfrac, ele), dshape=dnvfrac_t) |
208 | + end if |
209 | + end if |
210 | + |
211 | ! Step 2: Set up test function |
212 | |
213 | select case(stabilisation_scheme) |
214 | @@ -763,13 +819,13 @@ |
215 | ! Step 3: Assemble contributions |
216 | |
217 | ! Mass |
218 | - if(have_mass) call add_mass_element_cg(ele, test_function, t, density, olddensity, porosity_theta, detwei, detwei_old, detwei_new, matrix_addto, rhs_addto) |
219 | + if(have_mass) call add_mass_element_cg(ele, test_function, t, density, olddensity, porosity_theta, nvfrac, detwei, detwei_old, detwei_new, matrix_addto, rhs_addto) |
220 | |
221 | ! Advection |
222 | if(have_advection) call add_advection_element_cg(ele, test_function, t, & |
223 | velocity, grid_velocity, diffusivity, & |
224 | - density, olddensity, & |
225 | - dt_t, du_t, dug_t, drho_t, detwei, j_mat, matrix_addto, rhs_addto) |
226 | + density, olddensity, nvfrac, & |
227 | + dt_t, du_t, dug_t, drho_t, dnvfrac_t, detwei, j_mat, matrix_addto, rhs_addto) |
228 | |
229 | ! Absorption |
230 | if(have_absorption) call add_absorption_element_cg(ele, test_function, t, absorption, detwei, matrix_addto, rhs_addto) |
231 | @@ -784,8 +840,8 @@ |
232 | |
233 | ! Pressure |
234 | if(equation_type==FIELD_EQUATION_INTERNALENERGY) call add_pressurediv_element_cg(ele, test_function, t, & |
235 | - velocity, pressure, & |
236 | - du_t, detwei, rhs_addto) |
237 | + velocity, pressure, nvfrac, & |
238 | + du_t, dnvfrac_t, detwei, rhs_addto) |
239 | |
240 | ! Step 4: Insertion |
241 | |
242 | @@ -795,12 +851,13 @@ |
243 | |
244 | end subroutine assemble_advection_diffusion_element_cg |
245 | |
246 | - subroutine add_mass_element_cg(ele, test_function, t, density, olddensity, porosity_theta, detwei, detwei_old, detwei_new, matrix_addto, rhs_addto) |
247 | + subroutine add_mass_element_cg(ele, test_function, t, density, olddensity, porosity_theta, nvfrac, detwei, detwei_old, detwei_new, matrix_addto, rhs_addto) |
248 | integer, intent(in) :: ele |
249 | type(element_type), intent(in) :: test_function |
250 | type(scalar_field), intent(in) :: t |
251 | type(scalar_field), intent(in) :: density, olddensity |
252 | type(scalar_field), intent(in) :: porosity_theta |
253 | + type(scalar_field), intent(in) :: nvfrac |
254 | real, dimension(ele_ngi(t, ele)), intent(in) :: detwei, detwei_old, detwei_new |
255 | real, dimension(ele_loc(t, ele), ele_loc(t, ele)), intent(inout) :: matrix_addto |
256 | real, dimension(ele_loc(t, ele)), intent(inout) :: rhs_addto |
257 | @@ -827,6 +884,8 @@ |
258 | else |
259 | if (include_porosity) then |
260 | mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei*density_at_quad*porosity_theta_at_quad) |
261 | + else if(multiphase) then |
262 | + mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei*density_at_quad*ele_val_at_quad(nvfrac, ele)) |
263 | else |
264 | mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei*density_at_quad) |
265 | end if |
266 | @@ -877,8 +936,8 @@ |
267 | |
268 | subroutine add_advection_element_cg(ele, test_function, t, & |
269 | velocity, grid_velocity, diffusivity, & |
270 | - density, olddensity, & |
271 | - dt_t, du_t, dug_t, drho_t, detwei, j_mat, matrix_addto, rhs_addto) |
272 | + density, olddensity, nvfrac, & |
273 | + dt_t, du_t, dug_t, drho_t, dnvfrac_t, detwei, j_mat, matrix_addto, rhs_addto) |
274 | integer, intent(in) :: ele |
275 | type(element_type), intent(in) :: test_function |
276 | type(scalar_field), intent(in) :: t |
277 | @@ -886,10 +945,12 @@ |
278 | type(vector_field), pointer :: grid_velocity |
279 | type(tensor_field), intent(in) :: diffusivity |
280 | type(scalar_field), intent(in) :: density, olddensity |
281 | + type(scalar_field), intent(in) :: nvfrac |
282 | real, dimension(ele_loc(t, ele), ele_ngi(t, ele), mesh_dim(t)), intent(in) :: dt_t |
283 | real, dimension(ele_loc(velocity, ele), ele_ngi(velocity, ele), mesh_dim(t)) :: du_t |
284 | real, dimension(:, :, :) :: dug_t |
285 | real, dimension(ele_loc(density, ele), ele_ngi(density, ele), mesh_dim(density)), intent(in) :: drho_t |
286 | + real, dimension(:, :, :), intent(in) :: dnvfrac_t |
287 | real, dimension(ele_ngi(t, ele)), intent(in) :: detwei |
288 | real, dimension(mesh_dim(t), mesh_dim(t), ele_ngi(t, ele)), intent(in) :: j_mat |
289 | real, dimension(ele_loc(t, ele), ele_loc(t, ele)), intent(inout) :: matrix_addto |
290 | @@ -904,6 +965,10 @@ |
291 | real, dimension(velocity%dim, ele_ngi(density, ele)) :: densitygrad_at_quad |
292 | real, dimension(ele_ngi(density, ele)) :: udotgradrho_at_quad |
293 | |
294 | + real, dimension(ele_ngi(t, ele)) :: nvfrac_at_quad |
295 | + real, dimension(velocity%dim, ele_ngi(t, ele)) :: nvfracgrad_at_quad |
296 | + real, dimension(ele_ngi(t, ele)) :: udotgradnvfrac_at_quad |
297 | + |
298 | assert(have_advection) |
299 | |
300 | t_shape => ele_shape(t, ele) |
301 | @@ -922,6 +987,10 @@ |
302 | densitygrad_at_quad = density_theta*ele_grad_at_quad(density, ele, drho_t) & |
303 | +(1.-density_theta)*ele_grad_at_quad(olddensity, ele, drho_t) |
304 | udotgradrho_at_quad = sum(densitygrad_at_quad*velocity_at_quad, 1) |
305 | + |
306 | + if(multiphase) then |
307 | + nvfrac_at_quad = ele_val_at_quad(nvfrac, ele) |
308 | + end if |
309 | end select |
310 | |
311 | if(integrate_advection_by_parts) then |
312 | @@ -931,12 +1000,26 @@ |
313 | ! / / |
314 | select case(equation_type) |
315 | case(FIELD_EQUATION_INTERNALENERGY) |
316 | - advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei*density_at_quad) |
317 | + if(multiphase) then |
318 | + advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei*density_at_quad*nvfrac_at_quad) |
319 | + else |
320 | + advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei*density_at_quad) |
321 | + end if |
322 | + |
323 | if(abs(1.0 - beta) > epsilon(0.0)) then |
324 | velocity_div_at_quad = ele_div_at_quad(velocity, ele, du_t) |
325 | - advection_mat = advection_mat & |
326 | - - (1.0-beta) * shape_shape(test_function, t_shape, (velocity_div_at_quad*density_at_quad & |
327 | - +udotgradrho_at_quad)* detwei) |
328 | + |
329 | + if(multiphase) then |
330 | + |
331 | + nvfracgrad_at_quad = ele_grad_at_quad(nvfrac, ele, dnvfrac_t) |
332 | + udotgradnvfrac_at_quad = sum(nvfracgrad_at_quad*velocity_at_quad, 1) |
333 | + |
334 | + advection_mat = advection_mat - (1.0-beta) * ( shape_shape(test_function, t_shape, detwei*velocity_div_at_quad*density_at_quad*nvfrac_at_quad) & |
335 | + + shape_shape(test_function, t_shape, detwei*nvfrac_at_quad*udotgradrho_at_quad) & |
336 | + + shape_shape(test_function, t_shape, detwei*density_at_quad*udotgradnvfrac_at_quad) ) |
337 | + else |
338 | + advection_mat = advection_mat - (1.0-beta) * shape_shape(test_function, t_shape, (velocity_div_at_quad*density_at_quad + udotgradrho_at_quad)*detwei) |
339 | + end if |
340 | end if |
341 | case default |
342 | advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei) |
343 | @@ -953,12 +1036,31 @@ |
344 | ! / / |
345 | select case(equation_type) |
346 | case(FIELD_EQUATION_INTERNALENERGY) |
347 | - advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei*density_at_quad) |
348 | + |
349 | + if(multiphase) then |
350 | + ! vfrac*rho*nu*grad(internalenergy) |
351 | + advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei*density_at_quad*nvfrac_at_quad) |
352 | + else |
353 | + advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei*density_at_quad) |
354 | + end if |
355 | + |
356 | if(abs(beta) > epsilon(0.0)) then |
357 | velocity_div_at_quad = ele_div_at_quad(velocity, ele, du_t) |
358 | - advection_mat = advection_mat & |
359 | - + beta*shape_shape(test_function, t_shape, (velocity_div_at_quad*density_at_quad & |
360 | - +udotgradrho_at_quad)*detwei) |
361 | + |
362 | + if(multiphase) then |
363 | + ! advection_mat + internalenergy*div(vfrac*rho*nu) |
364 | + ! Split up div(vfrac*rho*nu) = vfrac*rho*div(nu) + nu*grad(vfrac*rho) = vfrac*rho*div(nu) + nu*(vfrac*grad(rho) + rho*grad(nvfrac)) |
365 | + |
366 | + nvfracgrad_at_quad = ele_grad_at_quad(nvfrac, ele, dnvfrac_t) |
367 | + udotgradnvfrac_at_quad = sum(nvfracgrad_at_quad*velocity_at_quad, 1) |
368 | + |
369 | + advection_mat = advection_mat + beta * ( shape_shape(test_function, t_shape, detwei*velocity_div_at_quad*density_at_quad*nvfrac_at_quad) & |
370 | + + shape_shape(test_function, t_shape, detwei*nvfrac_at_quad*udotgradrho_at_quad) & |
371 | + + shape_shape(test_function, t_shape, detwei*density_at_quad*udotgradnvfrac_at_quad) ) |
372 | + else |
373 | + advection_mat = advection_mat + beta*shape_shape(test_function, t_shape, (velocity_div_at_quad*density_at_quad & |
374 | + +udotgradrho_at_quad)*detwei) |
375 | + end if |
376 | end if |
377 | case default |
378 | advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei) |
379 | @@ -1060,26 +1162,39 @@ |
380 | |
381 | end subroutine add_diffusivity_element_cg |
382 | |
383 | - subroutine add_pressurediv_element_cg(ele, test_function, t, velocity, pressure, du_t, detwei, rhs_addto) |
384 | + subroutine add_pressurediv_element_cg(ele, test_function, t, velocity, pressure, nvfrac, du_t, dnvfrac_t, detwei, rhs_addto) |
385 | |
386 | integer, intent(in) :: ele |
387 | type(element_type), intent(in) :: test_function |
388 | type(scalar_field), intent(in) :: t |
389 | type(vector_field), intent(in) :: velocity |
390 | type(scalar_field), intent(in) :: pressure |
391 | - real, dimension(ele_loc(velocity, ele), ele_ngi(velocity, ele), mesh_dim(t)) :: du_t |
392 | + type(scalar_field), intent(in) :: nvfrac |
393 | + real, dimension(ele_loc(velocity, ele), ele_ngi(velocity, ele), mesh_dim(t)), intent(in) :: du_t |
394 | + real, dimension(:, :, :), intent(in) :: dnvfrac_t |
395 | real, dimension(ele_ngi(t, ele)), intent(in) :: detwei |
396 | real, dimension(ele_loc(t, ele)), intent(inout) :: rhs_addto |
397 | |
398 | + real, dimension(velocity%dim, ele_ngi(t, ele)) :: nvfracgrad_at_quad |
399 | + real, dimension(ele_ngi(t, ele)) :: udotgradnvfrac_at_quad |
400 | + |
401 | assert(equation_type==FIELD_EQUATION_INTERNALENERGY) |
402 | assert(ele_ngi(pressure, ele)==ele_ngi(t, ele)) |
403 | |
404 | - rhs_addto = rhs_addto - & |
405 | - shape_rhs(test_function, ele_div_at_quad(velocity, ele, du_t) * ele_val_at_quad(pressure, ele) * detwei) |
406 | + if(multiphase) then |
407 | + ! -p * (div(vfrac*nu)) |
408 | + nvfracgrad_at_quad = ele_grad_at_quad(nvfrac, ele, dnvfrac_t) |
409 | + udotgradnvfrac_at_quad = sum(nvfracgrad_at_quad*ele_val_at_quad(velocity, ele), 1) |
410 | + |
411 | + rhs_addto = rhs_addto - ( shape_rhs(test_function, ele_div_at_quad(velocity, ele, du_t) * ele_val_at_quad(nvfrac, ele) * ele_val_at_quad(pressure, ele) * detwei) & |
412 | + + shape_rhs(test_function, udotgradnvfrac_at_quad * ele_val_at_quad(pressure, ele) * detwei) ) |
413 | + else |
414 | + rhs_addto = rhs_addto - shape_rhs(test_function, ele_div_at_quad(velocity, ele, du_t) * ele_val_at_quad(pressure, ele) * detwei) |
415 | + end if |
416 | |
417 | end subroutine add_pressurediv_element_cg |
418 | |
419 | - subroutine assemble_advection_diffusion_face_cg(face, bc_type, t, t_bc, t_bc_2, matrix, rhs, positions, velocity, grid_velocity, density, olddensity) |
420 | + subroutine assemble_advection_diffusion_face_cg(face, bc_type, t, t_bc, t_bc_2, matrix, rhs, positions, velocity, grid_velocity, density, olddensity, nvfrac) |
421 | integer, intent(in) :: face |
422 | integer, intent(in) :: bc_type |
423 | type(scalar_field), intent(in) :: t |
424 | @@ -1092,6 +1207,7 @@ |
425 | type(vector_field), pointer :: grid_velocity |
426 | type(scalar_field), intent(in) :: density |
427 | type(scalar_field), intent(in) :: olddensity |
428 | + type(scalar_field), intent(in) :: nvfrac |
429 | |
430 | integer, dimension(face_loc(t, face)) :: face_nodes |
431 | real, dimension(face_ngi(t, face)) :: detwei |
432 | @@ -1125,7 +1241,7 @@ |
433 | |
434 | ! Advection |
435 | if(have_advection .and. integrate_advection_by_parts) & |
436 | - call add_advection_face_cg(face, bc_type, t, t_bc, velocity, grid_velocity, density, olddensity, detwei, normal, matrix_addto, rhs_addto) |
437 | + call add_advection_face_cg(face, bc_type, t, t_bc, velocity, grid_velocity, density, olddensity, nvfrac, detwei, normal, matrix_addto, rhs_addto) |
438 | |
439 | ! Diffusivity |
440 | if(have_diffusivity) call add_diffusivity_face_cg(face, bc_type, t, t_bc, t_bc_2, detwei, matrix_addto, rhs_addto) |
441 | @@ -1138,7 +1254,7 @@ |
442 | |
443 | end subroutine assemble_advection_diffusion_face_cg |
444 | |
445 | - subroutine add_advection_face_cg(face, bc_type, t, t_bc, velocity, grid_velocity, density, olddensity, detwei, normal, matrix_addto, rhs_addto) |
446 | + subroutine add_advection_face_cg(face, bc_type, t, t_bc, velocity, grid_velocity, density, olddensity, nvfrac, detwei, normal, matrix_addto, rhs_addto) |
447 | integer, intent(in) :: face |
448 | integer, intent(in) :: bc_type |
449 | type(scalar_field), intent(in) :: t |
450 | @@ -1147,6 +1263,7 @@ |
451 | type(vector_field), pointer :: grid_velocity |
452 | type(scalar_field), intent(in) :: density |
453 | type(scalar_field), intent(in) :: olddensity |
454 | + type(scalar_field), intent(in) :: nvfrac |
455 | real, dimension(face_ngi(t, face)), intent(in) :: detwei |
456 | real, dimension(mesh_dim(t), face_ngi(t, face)), intent(in) :: normal |
457 | real, dimension(face_loc(t, face), face_loc(t, face)), intent(inout) :: matrix_addto |
458 | @@ -1171,8 +1288,12 @@ |
459 | case(FIELD_EQUATION_INTERNALENERGY) |
460 | density_at_quad = density_theta*face_val_at_quad(density, face) & |
461 | +(1.0-density_theta)*face_val_at_quad(olddensity, face) |
462 | - |
463 | - advection_mat = shape_shape(t_shape, t_shape, detwei * sum(velocity_at_quad * normal, 1) * density_at_quad) |
464 | + |
465 | + if(multiphase) then |
466 | + advection_mat = shape_shape(t_shape, t_shape, detwei * sum(velocity_at_quad * normal, 1) * density_at_quad * face_val_at_quad(nvfrac, face)) |
467 | + else |
468 | + advection_mat = shape_shape(t_shape, t_shape, detwei * sum(velocity_at_quad * normal, 1) * density_at_quad) |
469 | + end if |
470 | case default |
471 | |
472 | advection_mat = shape_shape(t_shape, t_shape, detwei * sum(velocity_at_quad * normal, 1)) |
473 | |
474 | === modified file 'assemble/Compressible_Projection.F90' |
475 | --- assemble/Compressible_Projection.F90 2012-01-20 15:54:44 +0000 |
476 | +++ assemble/Compressible_Projection.F90 2012-08-09 00:46:20 +0000 |
477 | @@ -39,6 +39,7 @@ |
478 | use fefields, only: compute_cv_mass |
479 | use state_fields_module |
480 | use upwind_stabilisation |
481 | + use multiphase_module |
482 | implicit none |
483 | |
484 | ! Buffer for output messages. |
485 | @@ -56,7 +57,10 @@ |
486 | integer :: nu_bar_scheme |
487 | real :: nu_bar_scale |
488 | |
489 | -contains |
490 | + !! Are we running a multiphase flow simulation? |
491 | + logical :: multiphase |
492 | + |
493 | + contains |
494 | |
495 | subroutine assemble_compressible_projection_cv(state, cmc, rhs, dt, theta_pg, theta_divergence, cmcget) |
496 | |
497 | @@ -367,11 +371,12 @@ |
498 | |
499 | end subroutine assemble_mmat_compressible_projection_cv |
500 | |
501 | - subroutine assemble_compressible_projection_cg(state, cmc, rhs, dt, theta_pg, theta_divergence, cmcget) |
502 | + subroutine assemble_compressible_projection_cg(state, istate, cmc, rhs, dt, theta_pg, theta_divergence, cmcget) |
503 | |
504 | ! inputs: |
505 | ! bucket full of fields |
506 | type(state_type), dimension(:), intent(inout) :: state |
507 | + integer, intent(in) :: istate |
508 | |
509 | type(csr_matrix), intent(inout) :: cmc |
510 | type(scalar_field), intent(inout) :: rhs |
511 | @@ -380,17 +385,25 @@ |
512 | real, intent(in) :: theta_pg, theta_divergence |
513 | logical, intent(in) :: cmcget |
514 | |
515 | - if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then |
516 | - |
517 | - call assemble_1mat_compressible_projection_cg(state(1), cmc, rhs, dt, & |
518 | - theta_pg, theta_divergence, cmcget) |
519 | - |
520 | + if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then |
521 | + multiphase = .true. |
522 | + call assemble_1mat_compressible_projection_cg(state(istate), cmc, rhs, dt, & |
523 | + theta_pg, theta_divergence, cmcget) |
524 | else |
525 | - |
526 | - FLExit("Multimaterial compressible continuous_galerkin pressure not possible.") |
527 | - |
528 | - end if |
529 | - |
530 | + multiphase = .false. |
531 | + |
532 | + if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then |
533 | + |
534 | + call assemble_1mat_compressible_projection_cg(state(1), cmc, rhs, dt, & |
535 | + theta_pg, theta_divergence, cmcget) |
536 | + |
537 | + else |
538 | + |
539 | + FLExit("Multimaterial compressible continuous_galerkin pressure not possible.") |
540 | + |
541 | + end if |
542 | + |
543 | + end if |
544 | |
545 | end subroutine assemble_compressible_projection_cg |
546 | |
547 | @@ -441,6 +454,11 @@ |
548 | logical :: have_absorption, have_source |
549 | integer :: stat |
550 | |
551 | + !! Multiphase variables |
552 | + ! Volume fraction fields |
553 | + type(scalar_field), pointer :: vfrac |
554 | + type(scalar_field) :: nvfrac |
555 | + |
556 | ! ============================================================= |
557 | ! Subroutine to construct the matrix CT_m (a.k.a. C1/2/3T). |
558 | ! ============================================================= |
559 | @@ -470,6 +488,15 @@ |
560 | |
561 | velocity=>extract_vector_field(state, "Velocity") |
562 | nonlinearvelocity=>extract_vector_field(state, "NonlinearVelocity") ! maybe this should be updated after the velocity solve? |
563 | + |
564 | + ! Get the non-linear PhaseVolumeFraction field if multiphase |
565 | + if(multiphase) then |
566 | + vfrac => extract_scalar_field(state, "PhaseVolumeFraction") |
567 | + call allocate(nvfrac, vfrac%mesh, "NonlinearPhaseVolumeFraction") |
568 | + call zero(nvfrac) |
569 | + call get_nonlinear_volume_fraction(state, nvfrac) |
570 | + ewrite_minmax(nvfrac) |
571 | + end if |
572 | |
573 | pressure => extract_scalar_field(state, "Pressure") |
574 | |
575 | @@ -552,7 +579,11 @@ |
576 | ! Important note: with SUPG the test function derivatives have not been |
577 | ! modified. |
578 | |
579 | - ele_mat = (1./(dt*dt*theta_divergence*theta_pg))*shape_shape(test_shape, test_shape_ptr, detwei*drhodp_at_quad) |
580 | + if(multiphase) then |
581 | + ele_mat = (1./(dt*dt*theta_divergence*theta_pg))*shape_shape(test_shape, test_shape_ptr, detwei*ele_val_at_quad(nvfrac, ele)*drhodp_at_quad) |
582 | + else |
583 | + ele_mat = (1./(dt*dt*theta_divergence*theta_pg))*shape_shape(test_shape, test_shape_ptr, detwei*drhodp_at_quad) |
584 | + end if |
585 | ! / |
586 | ! rhs = |test_shape* & |
587 | ! / |
588 | @@ -560,8 +591,14 @@ |
589 | ! +(absorption)*(drhodp*theta*(eospressure - (pressure + atmospheric_pressure)) |
590 | ! - theta*density - (1-theta)*olddensity) |
591 | ! +source)dV |
592 | - ele_rhs = (1./dt)*shape_rhs(test_shape, detwei*((drhodp_at_quad*(eosp_at_quad - p_at_quad)) & |
593 | - +(olddensity_at_quad - density_at_quad))) |
594 | + |
595 | + if(multiphase) then |
596 | + ele_rhs = (1./dt)*shape_rhs(test_shape, detwei*(ele_val_at_quad(nvfrac, ele))*((drhodp_at_quad*(eosp_at_quad - p_at_quad)) & |
597 | + +(olddensity_at_quad - density_at_quad))) |
598 | + else |
599 | + ele_rhs = (1./dt)*shape_rhs(test_shape, detwei*((drhodp_at_quad*(eosp_at_quad - p_at_quad)) & |
600 | + +(olddensity_at_quad - density_at_quad))) |
601 | + end if |
602 | |
603 | if(have_source) then |
604 | ele_rhs = ele_rhs + shape_rhs(test_shape, detwei*ele_val_at_quad(source, ele)) |
605 | @@ -587,6 +624,10 @@ |
606 | |
607 | call deallocate(drhodp) |
608 | call deallocate(eospressure) |
609 | + |
610 | + if(multiphase) then |
611 | + call deallocate(nvfrac) |
612 | + end if |
613 | |
614 | end if |
615 | |
616 | @@ -597,17 +638,29 @@ |
617 | type(state_type), dimension(:), intent(inout) :: state |
618 | |
619 | type(scalar_field), pointer :: density |
620 | - |
621 | - if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then |
622 | - |
623 | - density=>extract_scalar_field(state(1),'Density') |
624 | - |
625 | - if(have_option(trim(density%option_path)//"/prognostic")) then |
626 | - |
627 | - call compressible_eos(state(1), density=density) |
628 | - |
629 | - end if |
630 | - |
631 | + |
632 | + integer :: istate |
633 | + |
634 | + if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then |
635 | + do istate=1,size(state) |
636 | + density=>extract_scalar_field(state(istate),'Density') |
637 | + |
638 | + if(have_option(trim(density%option_path)//"/prognostic")) then |
639 | + call compressible_eos(state(istate), density=density) |
640 | + end if |
641 | + end do |
642 | + else |
643 | + if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then |
644 | + |
645 | + density=>extract_scalar_field(state(1),'Density') |
646 | + |
647 | + if(have_option(trim(density%option_path)//"/prognostic")) then |
648 | + |
649 | + call compressible_eos(state(1), density=density) |
650 | + |
651 | + end if |
652 | + |
653 | + end if |
654 | end if |
655 | |
656 | end subroutine update_compressible_density |
657 | @@ -620,8 +673,8 @@ |
658 | do i=0, option_count("/material_phase")-1 |
659 | prognostic_pressure_path="/material_phase"//int2str(i)//"/scalar_field::Pressure/prognostic" |
660 | if (have_option(trim(prognostic_pressure_path)//"/spatial_discretisation/discontinuous_galerkin") & |
661 | - .and. have_option(trim(prognostic_pressure_path)//"/scheme/use_compressible_projection")) then |
662 | - FLExit("With a DG pressure you cannot have use_compressible_projection") |
663 | + .and. have_option(trim(prognostic_pressure_path)//"/scheme/use_compressible_projection_method")) then |
664 | + FLExit("With a DG pressure you cannot have use_compressible_projection_method") |
665 | end if |
666 | end do |
667 | |
668 | |
669 | === modified file 'assemble/Diagnostic_Fields_Matrices.F90' |
670 | --- assemble/Diagnostic_Fields_Matrices.F90 2012-05-27 16:09:33 +0000 |
671 | +++ assemble/Diagnostic_Fields_Matrices.F90 2012-08-09 00:46:20 +0000 |
672 | @@ -40,7 +40,7 @@ |
673 | use spud |
674 | use parallel_fields, only: zero_non_owned |
675 | use divergence_matrix_cv, only: assemble_divergence_matrix_cv |
676 | - use divergence_matrix_cg, only: assemble_divergence_matrix_cg |
677 | + use divergence_matrix_cg, only: assemble_divergence_matrix_cg, assemble_compressible_divergence_matrix_cg |
678 | use gradient_matrix_cg, only: assemble_gradient_matrix_cg |
679 | use parallel_tools |
680 | use sparsity_patterns, only: make_sparsity |
681 | @@ -58,7 +58,8 @@ |
682 | |
683 | public :: calculate_divergence_cv, calculate_divergence_fe, & |
684 | calculate_div_t_cv, calculate_div_t_fe, & |
685 | - calculate_grad_fe, calculate_sum_velocity_divergence |
686 | + calculate_grad_fe, calculate_sum_velocity_divergence, & |
687 | + calculate_compressible_continuity_residual |
688 | |
689 | contains |
690 | |
691 | @@ -356,9 +357,9 @@ |
692 | |
693 | ! Allocate memory for matrices and sparsity patterns |
694 | call allocate(ctfield, sum_velocity_divergence%mesh, name="CTField") |
695 | - call zero (ctfield) |
696 | + call zero(ctfield) |
697 | call allocate(temp, sum_velocity_divergence%mesh, name="Temp") |
698 | - |
699 | + |
700 | ! Require a mass matrix type if not CV tested equation. |
701 | if (.not. test_with_cv_dual) then |
702 | mass_sparsity=make_sparsity(sum_velocity_divergence%mesh, sum_velocity_divergence%mesh, "MassSparsity") |
703 | @@ -382,27 +383,29 @@ |
704 | |
705 | ! Allocate sparsity patterns, C^T matrix and C^T RHS for current state |
706 | divergence_sparsity=make_sparsity(sum_velocity_divergence%mesh, u%mesh, "DivergenceSparsity") |
707 | + |
708 | allocate(ct_m) |
709 | call allocate(ct_m, divergence_sparsity, (/1, u%dim/), name="DivergenceMatrix" ) |
710 | call allocate(ct_rhs, sum_velocity_divergence%mesh, name="CTRHS") |
711 | - |
712 | + |
713 | ! Reassemble C^T matrix here |
714 | if (test_with_cv_dual) then |
715 | call assemble_divergence_matrix_cv(ct_m, state(i), ct_rhs=ct_rhs, & |
716 | test_mesh=sum_velocity_divergence%mesh, field=u) |
717 | else |
718 | - if(i==1) then ! Construct the mass matrix (just do this once) |
719 | + if(i==1) then ! Construct the mass matrix (just do this once) |
720 | call assemble_divergence_matrix_cg(ct_m, state(i), ct_rhs=ct_rhs, & |
721 | - test_mesh=sum_velocity_divergence%mesh, field=u, & |
722 | - option_path=sum_velocity_divergence%option_path, div_mass=mass) |
723 | + test_mesh=sum_velocity_divergence%mesh, field=u, & |
724 | + option_path=sum_velocity_divergence%option_path, div_mass=mass) |
725 | else |
726 | call assemble_divergence_matrix_cg(ct_m, state(i), ct_rhs=ct_rhs, & |
727 | - test_mesh=sum_velocity_divergence%mesh, field=u, & |
728 | - option_path=sum_velocity_divergence%option_path) |
729 | - end if |
730 | + test_mesh=sum_velocity_divergence%mesh, field=u, & |
731 | + option_path=sum_velocity_divergence%option_path) |
732 | + end if |
733 | + |
734 | end if |
735 | |
736 | - ! Construct the linear system of equations to solve for \sum{div(vfrac*u)} |
737 | + ! Construct the linear system of equations |
738 | call zero(temp) |
739 | call mult(temp, ct_m, u) |
740 | call addto(temp, ct_rhs, -1.0) |
741 | @@ -417,7 +420,8 @@ |
742 | |
743 | end do |
744 | |
745 | - ! Solve for sum_velocity_divergence ( = \sum{div(vfrac*u)} ) |
746 | + ! Solve for sum_velocity_divergence |
747 | + ! ( = \sum{div(vfrac*u)} for incompressible multiphase flows ) |
748 | if (test_with_cv_dual) then |
749 | ! get the cv mass matrix |
750 | cv_mass => get_cv_mass(state(1), sum_velocity_divergence%mesh) |
751 | @@ -442,4 +446,148 @@ |
752 | |
753 | end subroutine calculate_sum_velocity_divergence |
754 | |
755 | + |
756 | + subroutine calculate_compressible_continuity_residual(state, compressible_continuity_residual) |
757 | + !!< Calculates the residual of the continity equation used in compressible multiphase flow simulations: |
758 | + !!< vfrac_c*d(rho_c)/dt + div(rho_c*vfrac_c*u_c) + \sum_i{ rho_c*div(vfrac_i*u_i) } |
759 | + |
760 | + type(state_type), dimension(:), intent(inout) :: state |
761 | + type(scalar_field), pointer :: compressible_continuity_residual, density, olddensity, vfrac |
762 | + type(scalar_field) :: drhodt |
763 | + |
764 | + ! Local variables |
765 | + type(vector_field), pointer :: u, x |
766 | + integer :: i, stat, ele |
767 | + logical :: dg |
768 | + |
769 | + type(csr_sparsity) :: divergence_sparsity |
770 | + type(block_csr_matrix), pointer :: ct_m |
771 | + |
772 | + type(csr_sparsity) :: mass_sparsity |
773 | + type(csr_matrix) :: mass |
774 | + type(scalar_field) :: ctfield, ct_rhs, temp |
775 | + |
776 | + type(element_type) :: test_function |
777 | + type(element_type), pointer :: compressible_continuity_residual_shape |
778 | + integer, dimension(:), pointer :: compressible_continuity_residual_nodes |
779 | + real, dimension(:), allocatable :: detwei |
780 | + real, dimension(:), allocatable :: drhodt_addto |
781 | + |
782 | + real :: dt |
783 | + |
784 | + ewrite(1,*) 'Entering calculate_compressible_continuity_residual' |
785 | + |
786 | + ! Allocate memory for matrices and sparsity patterns |
787 | + call allocate(ctfield, compressible_continuity_residual%mesh, name="CTField") |
788 | + call zero(ctfield) |
789 | + call allocate(temp, compressible_continuity_residual%mesh, name="Temp") |
790 | + |
791 | + mass_sparsity=make_sparsity(compressible_continuity_residual%mesh, compressible_continuity_residual%mesh, "MassSparsity") |
792 | + call allocate(mass, mass_sparsity, name="MassMatrix") |
793 | + call zero(mass) |
794 | + |
795 | + call get_option("/timestepping/timestep", dt) |
796 | + |
797 | + ! Sum up over the div's |
798 | + do i = 1, size(state) |
799 | + u => extract_vector_field(state(i), "Velocity", stat) |
800 | + |
801 | + ! If there's no velocity then cycle |
802 | + if(stat/=0) cycle |
803 | + ! If this is an aliased velocity then cycle |
804 | + if(aliased(u)) cycle |
805 | + ! If the velocity isn't prognostic then cycle |
806 | + if(.not.have_option(trim(u%option_path)//"/prognostic")) cycle |
807 | + |
808 | + ! If velocity field is prognostic, begin calculations below |
809 | + x => extract_vector_field(state(i), "Coordinate") |
810 | + |
811 | + ! Allocate sparsity patterns, C^T matrix and C^T RHS for current state |
812 | + divergence_sparsity=make_sparsity(compressible_continuity_residual%mesh, u%mesh, "DivergenceSparsity") |
813 | + allocate(ct_m) |
814 | + call allocate(ct_m, divergence_sparsity, (/1, u%dim/), name="DivergenceMatrix") |
815 | + call allocate(ct_rhs, compressible_continuity_residual%mesh, name="CTRHS") |
816 | + |
817 | + ! Reassemble C^T matrix here |
818 | + if(i==1) then ! Construct the mass matrix (just do this once) |
819 | + call assemble_compressible_divergence_matrix_cg(ct_m, state, i, ct_rhs, div_mass=mass) |
820 | + else |
821 | + call assemble_compressible_divergence_matrix_cg(ct_m, state, i, ct_rhs) |
822 | + end if |
823 | + |
824 | + if(have_option("/material_phase::"//trim(state(i)%name)//"/equation_of_state/compressible")) then |
825 | + ! Get the time derivative term for the compressible phase's density, vfrac_c * d(rho_c)/dt |
826 | + |
827 | + call allocate(drhodt, compressible_continuity_residual%mesh, name="drhodt") |
828 | + |
829 | + density => extract_scalar_field(state(i), "Density", stat) |
830 | + olddensity => extract_scalar_field(state(i), "OldDensity", stat) |
831 | + vfrac => extract_scalar_field(state(i), "PhaseVolumeFraction") |
832 | + |
833 | + ! Assumes Density and OldDensity are on the same mesh as Pressure, |
834 | + ! as it should be according to the manual. |
835 | + call zero(drhodt) |
836 | + |
837 | + dg = continuity(compressible_continuity_residual) < 0 |
838 | + |
839 | + element_loop: do ele = 1, element_count(compressible_continuity_residual) |
840 | + |
841 | + if(.not.dg .or. (dg .and. element_owned(compressible_continuity_residual,ele))) then |
842 | + |
843 | + allocate(detwei(ele_ngi(compressible_continuity_residual, ele))) |
844 | + allocate(drhodt_addto(ele_loc(compressible_continuity_residual, ele))) |
845 | + |
846 | + compressible_continuity_residual_nodes => ele_nodes(compressible_continuity_residual, ele) |
847 | + compressible_continuity_residual_shape => ele_shape(compressible_continuity_residual, ele) |
848 | + test_function = compressible_continuity_residual_shape |
849 | + |
850 | + call transform_to_physical(x, ele, detwei=detwei) |
851 | + |
852 | + drhodt_addto = shape_rhs(test_function, detwei*(ele_val_at_quad(vfrac,ele)*(ele_val_at_quad(density,ele) - ele_val_at_quad(olddensity,ele))/dt)) |
853 | + |
854 | + call addto(drhodt, compressible_continuity_residual_nodes, drhodt_addto) |
855 | + |
856 | + deallocate(detwei) |
857 | + deallocate(drhodt_addto) |
858 | + end if |
859 | + |
860 | + end do element_loop |
861 | + |
862 | + call addto(ctfield, drhodt) |
863 | + |
864 | + call deallocate(drhodt) |
865 | + |
866 | + end if |
867 | + |
868 | + ! Construct the linear system of equations |
869 | + call zero(temp) |
870 | + call mult(temp, ct_m, u) |
871 | + call addto(temp, ct_rhs, -1.0) |
872 | + |
873 | + ! Now add it to the sum |
874 | + call addto(ctfield, temp) |
875 | + |
876 | + call deallocate(ct_m) |
877 | + deallocate(ct_m) |
878 | + call deallocate(ct_rhs) |
879 | + call deallocate(divergence_sparsity) |
880 | + |
881 | + end do |
882 | + |
883 | + ! Solve for compressible_continuity_residual |
884 | + call zero(compressible_continuity_residual) |
885 | + call petsc_solve(compressible_continuity_residual, mass, ctfield) |
886 | + ewrite_minmax(compressible_continuity_residual) |
887 | + |
888 | + ! Deallocate memory |
889 | + call deallocate(ctfield) |
890 | + call deallocate(temp) |
891 | + call deallocate(mass_sparsity) |
892 | + call deallocate(mass) |
893 | + |
894 | + ewrite(1,*) 'Exiting calculate_compressible_continuity_residual' |
895 | + |
896 | + end subroutine calculate_compressible_continuity_residual |
897 | + |
898 | + |
899 | end module diagnostic_fields_matrices |
900 | |
901 | === modified file 'assemble/Diagnostic_fields_wrapper.F90' |
902 | --- assemble/Diagnostic_fields_wrapper.F90 2012-02-18 15:45:48 +0000 |
903 | +++ assemble/Diagnostic_fields_wrapper.F90 2012-08-09 00:46:20 +0000 |
904 | @@ -585,6 +585,21 @@ |
905 | FLExit("The SumVelocityDivergence field is only used in multiphase simulations.") |
906 | end if |
907 | end if |
908 | + |
909 | + s_field => extract_scalar_field(state(i), "CompressibleContinuityResidual", stat) |
910 | + if(stat == 0) then |
911 | + ! Check that we are running a compressible multiphase simulation |
912 | + if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1 .and. option_count("/material_phase/equation_of_state/compressible") > 0) then |
913 | + diagnostic = have_option(trim(s_field%option_path)//"/diagnostic") |
914 | + if(diagnostic .and. .not.(aliased(s_field))) then |
915 | + if(recalculate(trim(s_field%option_path))) then |
916 | + call calculate_compressible_continuity_residual(state, s_field) |
917 | + end if |
918 | + end if |
919 | + else |
920 | + FLExit("The CompressibleContinuityResidual field is only used in compressible multiphase simulations.") |
921 | + end if |
922 | + end if |
923 | |
924 | ! end of fields that cannot be called through the generic |
925 | ! calculate_diagnostic_variable interface, i.e. - those that need things |
926 | |
927 | === modified file 'assemble/Divergence_Matrix_CG.F90' |
928 | --- assemble/Divergence_Matrix_CG.F90 2012-05-27 16:09:33 +0000 |
929 | +++ assemble/Divergence_Matrix_CG.F90 2012-08-09 00:46:20 +0000 |
930 | @@ -58,7 +58,10 @@ |
931 | integer :: nu_bar_scheme |
932 | real :: nu_bar_scale |
933 | |
934 | -contains |
935 | + !! Are we running a multiphase flow simulation? |
936 | + logical :: multiphase |
937 | + |
938 | + contains |
939 | |
940 | subroutine assemble_divergence_matrix_cg(CT_m, state, ct_rhs, & |
941 | test_mesh, field, option_path, & |
942 | @@ -119,7 +122,6 @@ |
943 | logical :: l_get_ct |
944 | |
945 | !! Multiphase variables |
946 | - logical :: multiphase |
947 | ! Volume fraction fields |
948 | type(scalar_field), pointer :: vfrac |
949 | type(scalar_field) :: nvfrac |
950 | @@ -355,40 +357,55 @@ |
951 | |
952 | end subroutine assemble_divergence_matrix_cg |
953 | |
954 | - subroutine assemble_compressible_divergence_matrix_cg(ctp_m, state, ct_rhs) |
955 | + subroutine assemble_compressible_divergence_matrix_cg(ctp_m, state, istate, ct_rhs, div_mass) |
956 | |
957 | ! inputs/outputs |
958 | ! bucket full of fields |
959 | type(state_type), dimension(:), intent(inout) :: state |
960 | - |
961 | + integer, intent(in) :: istate |
962 | ! the compressible divergence matrix |
963 | type(block_csr_matrix), pointer :: ctp_m |
964 | |
965 | type(scalar_field), intent(inout), optional :: ct_rhs |
966 | |
967 | - if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then |
968 | - |
969 | - call assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state(1), ct_rhs) |
970 | - |
971 | + type(csr_matrix), intent(inout), optional :: div_mass |
972 | + |
973 | + if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then |
974 | + multiphase = .true. |
975 | + call assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state, istate, ct_rhs, div_mass) |
976 | else |
977 | - |
978 | - FLExit("Multimaterial compressible continuous_galerkin pressure not possible.") |
979 | - |
980 | + multiphase = .false. |
981 | + |
982 | + if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then |
983 | + |
984 | + call assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state, 1, ct_rhs, div_mass) |
985 | + |
986 | + else |
987 | + |
988 | + FLExit("Multimaterial compressible continuous_galerkin pressure not possible.") |
989 | + |
990 | + end if |
991 | + |
992 | end if |
993 | + |
994 | + |
995 | |
996 | end subroutine assemble_compressible_divergence_matrix_cg |
997 | |
998 | - subroutine assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state, ct_rhs) |
999 | + subroutine assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state, istate, ct_rhs, div_mass) |
1000 | |
1001 | ! inputs/outputs |
1002 | ! bucket full of fields |
1003 | - type(state_type), intent(inout) :: state |
1004 | + type(state_type), dimension(:), intent(inout) :: state |
1005 | + integer, intent(in) :: istate |
1006 | |
1007 | ! the compressible divergence matrix |
1008 | type(block_csr_matrix), pointer :: ctp_m |
1009 | |
1010 | type(scalar_field), intent(inout), optional :: ct_rhs |
1011 | |
1012 | + type(csr_matrix), intent(inout), optional :: div_mass |
1013 | + |
1014 | ! local |
1015 | type(mesh_type), pointer :: test_mesh |
1016 | |
1017 | @@ -408,6 +425,8 @@ |
1018 | real, dimension(:), allocatable :: density_at_quad, olddensity_at_quad |
1019 | real, dimension(:,:), allocatable :: density_grad_at_quad, nlvelocity_at_quad |
1020 | |
1021 | + real, dimension(:,:), allocatable :: div_mass_mat |
1022 | + |
1023 | ! loop integers |
1024 | integer :: ele, sele, dim |
1025 | |
1026 | @@ -425,7 +444,16 @@ |
1027 | integer, dimension(:), allocatable :: density_bc_type |
1028 | type(scalar_field) :: density_bc |
1029 | |
1030 | - integer :: stat |
1031 | + integer :: i, stat |
1032 | + |
1033 | + !! Multiphase variables |
1034 | + ! Volume fraction fields |
1035 | + type(scalar_field), pointer :: vfrac |
1036 | + type(scalar_field) :: nvfrac |
1037 | + type(element_type), pointer :: nvfrac_shape |
1038 | + ! Transformed gradient function for the non-linear PhaseVolumeFraction. |
1039 | + real, dimension(:, :, :), allocatable :: dnvfrac_t |
1040 | + logical :: is_compressible_phase ! Is this the (single) compressible phase in a multiphase simulation? |
1041 | |
1042 | ! ============================================================= |
1043 | ! Subroutine to construct the matrix ctp_m (a.k.a. C1/2/3TP). |
1044 | @@ -433,16 +461,42 @@ |
1045 | |
1046 | ewrite(2,*) 'In assemble_1mat_compressible_divergence_matrix_cg' |
1047 | |
1048 | - coordinate=> extract_vector_field(state, "Coordinate") |
1049 | - |
1050 | - density => extract_scalar_field(state, "Density") |
1051 | - olddensity => extract_scalar_field(state, "OldDensity") |
1052 | - |
1053 | - pressure => extract_scalar_field(state, "Pressure") |
1054 | - |
1055 | - velocity=>extract_vector_field(state, "Velocity") |
1056 | - nonlinearvelocity=>extract_vector_field(state, "NonlinearVelocity") ! maybe this should be updated after the velocity solve? |
1057 | - |
1058 | + coordinate=> extract_vector_field(state(istate), "Coordinate") |
1059 | + |
1060 | + density => extract_scalar_field(state(istate), "Density") |
1061 | + olddensity => extract_scalar_field(state(istate), "OldDensity") |
1062 | + |
1063 | + if(have_option(trim(state(istate)%option_path)//"/equation_of_state/compressible")) then |
1064 | + is_compressible_phase = .true. |
1065 | + else |
1066 | + is_compressible_phase = .false. |
1067 | + |
1068 | + ! Find the compressible phase and extract it's density to form rho_c * div(vfrac_i*u_i), where _c and _i represent |
1069 | + ! the compressible and incompressible phases respectively. |
1070 | + do i = 1, size(state) |
1071 | + if(have_option(trim(state(i)%option_path)//"/equation_of_state/compressible")) then |
1072 | + density => extract_scalar_field(state(i), "Density") |
1073 | + olddensity => extract_scalar_field(state(i), "OldDensity") |
1074 | + end if |
1075 | + end do |
1076 | + end if |
1077 | + |
1078 | + pressure => extract_scalar_field(state(istate), "Pressure") |
1079 | + |
1080 | + velocity=>extract_vector_field(state(istate), "Velocity") |
1081 | + nonlinearvelocity=>extract_vector_field(state(istate), "NonlinearVelocity") ! maybe this should be updated after the velocity solve? |
1082 | + |
1083 | + ! Get the non-linear PhaseVolumeFraction field if multiphase |
1084 | + if(multiphase) then |
1085 | + vfrac => extract_scalar_field(state(istate), "PhaseVolumeFraction") |
1086 | + call allocate(nvfrac, vfrac%mesh, "NonlinearPhaseVolumeFraction") |
1087 | + call zero(nvfrac) |
1088 | + call get_nonlinear_volume_fraction(state(istate), nvfrac) |
1089 | + ewrite_minmax(nvfrac) |
1090 | + else |
1091 | + nullify(vfrac) |
1092 | + end if |
1093 | + |
1094 | integrate_by_parts=have_option(trim(complete_field_path(density%option_path, stat=stat))//& |
1095 | &"/spatial_discretisation/continuous_galerkin/advection_terms/integrate_advection_by_parts")& |
1096 | .or. have_option(trim(complete_field_path(velocity%option_path, stat=stat))//& |
1097 | @@ -474,7 +528,7 @@ |
1098 | end if |
1099 | |
1100 | call get_option(trim(complete_field_path(density%option_path, stat=stat))//& |
1101 | - &"/temporal_discretisation/theta", theta) |
1102 | + &"/temporal_discretisation/theta", theta) |
1103 | call get_option("/timestepping/timestep", dt) |
1104 | |
1105 | test_mesh => pressure%mesh |
1106 | @@ -494,8 +548,14 @@ |
1107 | olddensity_at_quad(ele_ngi(density, 1)), & |
1108 | nlvelocity_at_quad(nonlinearvelocity%dim, ele_ngi(nonlinearvelocity, 1)), & |
1109 | density_grad_at_quad(field%dim, ele_ngi(density,1)), & |
1110 | - j_mat(field%dim, field%dim, ele_ngi(density, 1))) |
1111 | + j_mat(field%dim, field%dim, ele_ngi(density, 1)), & |
1112 | + div_mass_mat(ele_loc(test_mesh, 1), ele_loc(test_mesh, 1))) |
1113 | |
1114 | + if(multiphase) then |
1115 | + ! We will need grad(nvfrac) if we are not integrating by parts below |
1116 | + allocate(dnvfrac_t(ele_loc(nvfrac,1), ele_ngi(nvfrac,1), field%dim)) |
1117 | + end if |
1118 | + |
1119 | do ele=1, element_count(test_mesh) |
1120 | |
1121 | test_nodes=>ele_nodes(test_mesh, ele) |
1122 | @@ -509,7 +569,7 @@ |
1123 | olddensity_at_quad = ele_val_at_quad(olddensity, ele) |
1124 | |
1125 | nlvelocity_at_quad = ele_val_at_quad(nonlinearvelocity, ele) |
1126 | - |
1127 | + |
1128 | if(any(stabilisation_scheme == (/STABILISATION_STREAMLINE_UPWIND, STABILISATION_SUPG/))) then |
1129 | call transform_to_physical(coordinate, ele, test_shape_ptr, dshape = dtest_t, & |
1130 | detwei = detwei, j = j_mat) |
1131 | @@ -517,7 +577,7 @@ |
1132 | call transform_to_physical(coordinate, ele, test_shape_ptr, dshape = dtest_t, detwei=detwei) |
1133 | end if |
1134 | |
1135 | - if(.not.integrate_by_parts) then |
1136 | + if(.not.integrate_by_parts .or. (multiphase .and. integrate_by_parts .and. .not.is_compressible_phase)) then |
1137 | ! transform the field (velocity) derivatives into physical space |
1138 | call transform_to_physical(coordinate, ele, field_shape, dshape=dfield_t) |
1139 | |
1140 | @@ -545,16 +605,53 @@ |
1141 | |
1142 | |
1143 | if(integrate_by_parts) then |
1144 | + |
1145 | ! if SUPG is fixed for P>1 then this dtest_t should be updated |
1146 | - ele_mat = -dshape_shape(dtest_t, field_shape, & |
1147 | - detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad)) |
1148 | + if(multiphase .and. .not.is_compressible_phase) then |
1149 | + density_grad_at_quad = theta*(ele_grad_at_quad(density, ele, ddensity_t))+& |
1150 | + (1-theta)*(ele_grad_at_quad(olddensity, ele, ddensity_t)) |
1151 | + |
1152 | + ele_mat = -dshape_shape(dtest_t, field_shape, detwei*ele_val_at_quad(nvfrac, ele)*(theta*density_at_quad + (1-theta)*olddensity_at_quad)) - shape_shape_vector(test_shape, field_shape, detwei*ele_val_at_quad(nvfrac, ele), density_grad_at_quad) |
1153 | + |
1154 | + else if(multiphase .and. is_compressible_phase) then |
1155 | + ele_mat = -dshape_shape(dtest_t, field_shape, detwei*ele_val_at_quad(nvfrac, ele)*(theta*density_at_quad + (1-theta)*olddensity_at_quad)) |
1156 | + |
1157 | + else |
1158 | + |
1159 | + ele_mat = -dshape_shape(dtest_t, field_shape, detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad)) |
1160 | + end if |
1161 | else |
1162 | density_grad_at_quad = theta*(ele_grad_at_quad(density, ele, ddensity_t))+& |
1163 | (1-theta)*(ele_grad_at_quad(olddensity, ele, ddensity_t)) |
1164 | |
1165 | - ele_mat = shape_dshape(test_shape, dfield_t, & |
1166 | + if(multiphase) then |
1167 | + |
1168 | + ! If the field and nvfrac meshes are different, then we need to |
1169 | + ! compute the derivatives of the nvfrac shape functions. |
1170 | + if(.not.(nvfrac%mesh == field%mesh)) then |
1171 | + nvfrac_shape => ele_shape(nvfrac%mesh, ele) |
1172 | + call transform_to_physical(coordinate, ele, nvfrac_shape, dshape=dnvfrac_t) |
1173 | + else |
1174 | + dnvfrac_t = dfield_t |
1175 | + end if |
1176 | + |
1177 | + ! Split up the divergence term div(rho*vfrac*u) = rho*div(u*vfrac) + vfrac*u*grad(rho) |
1178 | + ! = (rho*vfrac*div(u) + rho*u*grad(vfrac)) + u*vfrac*grad(rho) |
1179 | + |
1180 | + ! First assemble rho*div(u*vfrac). This is the incompressible phase's divergence matrix. |
1181 | + ele_mat = shape_dshape(test_shape, dfield_t, detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad)*ele_val_at_quad(nvfrac, ele)) + shape_shape_vector(test_shape, field_shape, detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad), ele_grad_at_quad(nvfrac, ele, dnvfrac_t)) |
1182 | + |
1183 | + ! If the phase is compressible, then we now complete the assembly of div(rho*vfrac*u) below. |
1184 | + if(is_compressible_phase) then |
1185 | + ele_mat = ele_mat + shape_shape_vector(test_shape, field_shape, detwei*ele_val_at_quad(nvfrac, ele), density_grad_at_quad) |
1186 | + end if |
1187 | + |
1188 | + else |
1189 | + ele_mat = shape_dshape(test_shape, dfield_t, & |
1190 | detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad)) + & |
1191 | shape_shape_vector(test_shape, field_shape, detwei, density_grad_at_quad) |
1192 | + end if |
1193 | + |
1194 | end if |
1195 | |
1196 | ! Stabilisation does not return the right shape for this operator! |
1197 | @@ -568,6 +665,11 @@ |
1198 | do dim = 1, field%dim |
1199 | call addto(ctp_m, 1, dim, test_nodes, field_nodes, ele_mat(dim,:,:)) |
1200 | end do |
1201 | + |
1202 | + if(present(div_mass)) then |
1203 | + div_mass_mat = shape_shape(test_shape, test_shape, detwei) |
1204 | + call addto(div_mass, test_nodes, test_nodes, div_mass_mat) |
1205 | + end if |
1206 | |
1207 | call deallocate(test_shape) |
1208 | |
1209 | @@ -618,8 +720,14 @@ |
1210 | & detwei_f=detwei_bdy,& |
1211 | & normal=normal_bdy) |
1212 | |
1213 | - ele_mat_bdy = shape_shape_vector(test_shape_ptr, field_shape, & |
1214 | - detwei_bdy*density_bdy, normal_bdy) |
1215 | + |
1216 | + if(multiphase) then |
1217 | + ele_mat_bdy = shape_shape_vector(test_shape_ptr, field_shape, & |
1218 | + detwei_bdy*density_bdy*face_val_at_quad(nvfrac, sele), normal_bdy) |
1219 | + else |
1220 | + ele_mat_bdy = shape_shape_vector(test_shape_ptr, field_shape, & |
1221 | + detwei_bdy*density_bdy, normal_bdy) |
1222 | + end if |
1223 | |
1224 | do dim = 1, field%dim |
1225 | if((field_bc_type(dim, sele)==1).and.present(ct_rhs)) then |
1226 | @@ -641,6 +749,11 @@ |
1227 | deallocate(test_nodes_bdy, field_nodes_bdy) |
1228 | |
1229 | end if |
1230 | + |
1231 | + if(multiphase) then |
1232 | + deallocate(dnvfrac_t) |
1233 | + call deallocate(nvfrac) |
1234 | + end if |
1235 | |
1236 | end subroutine assemble_1mat_compressible_divergence_matrix_cg |
1237 | |
1238 | |
1239 | === modified file 'assemble/Makefile.dependencies' |
1240 | --- assemble/Makefile.dependencies 2012-08-02 13:03:52 +0000 |
1241 | +++ assemble/Makefile.dependencies 2012-08-09 00:46:20 +0000 |
1242 | @@ -79,12 +79,12 @@ |
1243 | ../include/boundary_conditions_from_options.mod ../include/colouring.mod \ |
1244 | ../include/elements.mod ../include/fdebug.h ../include/field_options.mod \ |
1245 | ../include/fields.mod ../include/fldebug.mod \ |
1246 | - ../include/global_parameters.mod ../include/petsc_solve_state_module.mod \ |
1247 | - ../include/porous_media.mod ../include/profiler.mod \ |
1248 | - ../include/quadrature.mod ../include/sparse_tools.mod \ |
1249 | - ../include/sparse_tools_petsc.mod ../include/sparsity_patterns_meshes.mod \ |
1250 | - ../include/spud.mod ../include/state_module.mod \ |
1251 | - ../include/upwind_stabilisation.mod |
1252 | + ../include/global_parameters.mod ../include/multiphase_module.mod \ |
1253 | + ../include/petsc_solve_state_module.mod ../include/porous_media.mod \ |
1254 | + ../include/profiler.mod ../include/quadrature.mod \ |
1255 | + ../include/sparse_tools.mod ../include/sparse_tools_petsc.mod \ |
1256 | + ../include/sparsity_patterns_meshes.mod ../include/spud.mod \ |
1257 | + ../include/state_module.mod ../include/upwind_stabilisation.mod |
1258 | |
1259 | ../include/advection_diffusion_dg.mod: Advection_Diffusion_DG.o |
1260 | @true |
1261 | @@ -154,10 +154,10 @@ |
1262 | Compressible_Projection.F90 ../include/equation_of_state.mod \ |
1263 | ../include/fdebug.h ../include/fefields.mod ../include/field_options.mod \ |
1264 | ../include/fields.mod ../include/fldebug.mod \ |
1265 | - ../include/global_parameters.mod ../include/sparse_matrices_fields.mod \ |
1266 | - ../include/sparse_tools.mod ../include/spud.mod \ |
1267 | - ../include/state_fields_module.mod ../include/state_module.mod \ |
1268 | - ../include/upwind_stabilisation.mod |
1269 | + ../include/global_parameters.mod ../include/multiphase_module.mod \ |
1270 | + ../include/sparse_matrices_fields.mod ../include/sparse_tools.mod \ |
1271 | + ../include/spud.mod ../include/state_fields_module.mod \ |
1272 | + ../include/state_module.mod ../include/upwind_stabilisation.mod |
1273 | |
1274 | ../include/coriolis_module.mod: Coriolis.o |
1275 | @true |
1276 | |
1277 | === modified file 'assemble/Momentum_Equation.F90' |
1278 | --- assemble/Momentum_Equation.F90 2012-04-16 14:55:55 +0000 |
1279 | +++ assemble/Momentum_Equation.F90 2012-08-09 00:46:20 +0000 |
1280 | @@ -92,7 +92,7 @@ |
1281 | logical :: assemble_schur_auxiliary_matrix |
1282 | |
1283 | ! Do we want to use the compressible projection method? |
1284 | - logical :: use_compressible_projection |
1285 | + logical :: compressible_eos |
1286 | ! Are we doing a full Schur solve? |
1287 | logical :: full_schur |
1288 | ! Are we lumping mass or assuming consistent mass? |
1289 | @@ -282,7 +282,7 @@ |
1290 | multiphase = .false. |
1291 | end if |
1292 | ! Do we have fluid-particle drag (for multi-phase simulations)? |
1293 | - have_fp_drag = option_count("/material_phase/multiphase_properties/particle_diameter") > 0 |
1294 | + have_fp_drag = have_option("/multiphase_interaction/fluid_particle_drag") |
1295 | |
1296 | ! Get the pressure p^{n}, and get the assembly options for the divergence and CMC matrices |
1297 | ! find the first non-aliased pressure |
1298 | @@ -388,7 +388,7 @@ |
1299 | ! get the CV tested pressure gradient matrix (i.e. the divergence matrix) |
1300 | ! if required with a different unique name. Note there is no need |
1301 | ! to again decide reassemble_ct_m as ctp_m for this case is assembled when ct_m is. |
1302 | - if ((.not. use_compressible_projection) .and. cg_pressure_cv_test_continuity) then |
1303 | + if ((.not. compressible_eos) .and. cg_pressure_cv_test_continuity) then |
1304 | ctp_m(istate)%ptr => get_velocity_divergence_matrix(state(istate), ct_m_name = "CVTestedVelocityDivergenceMatrix") |
1305 | end if |
1306 | |
1307 | @@ -467,17 +467,16 @@ |
1308 | end select |
1309 | end if |
1310 | |
1311 | - if (has_boundary_condition(u, "free_surface") .or. use_compressible_projection) then |
1312 | - ! this needs fixing for multiphase theta_pg could in principle be chosen |
1313 | + if (has_boundary_condition(u, "free_surface") .or. compressible_eos) then |
1314 | + ! This needs fixing for multiphase. theta_pg could in principle be chosen |
1315 | ! per phase but then we need an array and we'd have to include theta_pg |
1316 | ! in cmc_m, i.e. solve for theta_div*dt*dp instead of theta_div*theta_pg*dt*dp |
1317 | ! theta_div can only be set once, but where and what is the default? |
1318 | - if (multiphase) then |
1319 | - FLExit("Multiphase does not work with a free surface or compressible flow.") |
1320 | + if (has_boundary_condition(u, "free_surface") .and. multiphase) then |
1321 | + FLExit("Multiphase does not work with a free surface.") |
1322 | end if |
1323 | |
1324 | - |
1325 | - ! With free surface or compressible-projection pressures are at integer |
1326 | + ! With free surface or compressible-projection, pressures are at integer |
1327 | ! time levels and we apply a theta-weighting to the pressure gradient term |
1328 | ! Also, obtain theta-weighting to be used in divergence term |
1329 | call get_option( trim(u%option_path)//'/prognostic/temporal_discretisation/theta', & |
1330 | @@ -491,6 +490,15 @@ |
1331 | ewrite(2,*) "theta_pg: ", theta_pg |
1332 | ewrite(2,*) "Velocity divergence is evaluated at n+theta_divergence" |
1333 | ewrite(2,*) "theta_divergence: ", theta_divergence |
1334 | + |
1335 | + ! Note: Compressible multiphase simulations work, but only when use_theta_pg and use_theta_divergence |
1336 | + ! are false. This needs improving - see comment above. |
1337 | + if(compressible_eos .and. multiphase .and. (use_theta_pg .or. use_theta_divergence)) then |
1338 | + ewrite(-1,*) "Currently, for compressible multiphase flow simulations, the" |
1339 | + ewrite(-1,*) "temporal_discretisation/theta and temporal_discretisation/theta_divergence values" |
1340 | + ewrite(-1,*) "for each Velocity field must be set to 1.0." |
1341 | + FLExit("Multiphase does not work when use_theta_pg or use_theta_divergence are true.") |
1342 | + end if |
1343 | else |
1344 | ! Pressures are, as usual, staggered in time with the velocities |
1345 | use_theta_pg=.false. |
1346 | @@ -514,6 +522,7 @@ |
1347 | ! Allocation of big_m |
1348 | if(dg(istate)) then |
1349 | call allocate_big_m_dg(state(istate), big_m(istate), u) |
1350 | + |
1351 | if(subcycle(istate)) then |
1352 | u_sparsity => get_csr_sparsity_firstorder(state, u%mesh, u%mesh) |
1353 | ! subcycle_m currently only contains advection, so diagonal=.true. |
1354 | @@ -531,8 +540,8 @@ |
1355 | ! Initialise the big_m, ct_m and ctp_m matrices |
1356 | call zero(big_m(istate)) |
1357 | if(reassemble_ct_m) then |
1358 | - call zero(ct_m(istate)%ptr) |
1359 | - if ((.not. use_compressible_projection) .and. cg_pressure_cv_test_continuity) then |
1360 | + call zero(ct_m(istate)%ptr) |
1361 | + if ((.not. compressible_eos) .and. cg_pressure_cv_test_continuity) then |
1362 | call zero(ctp_m(istate)%ptr) |
1363 | end if |
1364 | end if |
1365 | @@ -580,7 +589,7 @@ |
1366 | call subtract_geostrophic_pressure_gradient(mom_rhs(istate), state(istate)) |
1367 | end if |
1368 | else |
1369 | - ! This call will form the ct_rhs, which for use_compressible_projection |
1370 | + ! This call will form the ct_rhs, which for compressible_eos |
1371 | ! or cg_pressure_cv_test_continuity is formed for a second time later below. |
1372 | call construct_momentum_cg(u, p, density, x, & |
1373 | big_m(istate), mom_rhs(istate), ct_m(istate)%ptr, & |
1374 | @@ -635,7 +644,7 @@ |
1375 | |
1376 | call profiler_tic(p, "assembly") |
1377 | if(cv_pressure) then |
1378 | - ! This call will form the ct_rhs, which for use_compressible_projection |
1379 | + ! This call will form the ct_rhs, which for compressible_eos |
1380 | ! is formed for a second time later below. |
1381 | call assemble_divergence_matrix_cv(ct_m(istate)%ptr, state(istate), ct_rhs=ct_rhs(istate), & |
1382 | test_mesh=p%mesh, field=u, get_ct=reassemble_ct_m) |
1383 | @@ -643,7 +652,7 @@ |
1384 | |
1385 | ! Assemble divergence matrix C^T. |
1386 | ! At the moment cg does its own ct assembly. We might change this in the future. |
1387 | - ! This call will form the ct_rhs, which for use_compressible_projection |
1388 | + ! This call will form the ct_rhs, which for compressible_eos |
1389 | ! or cg_pressure_cv_test_continuity is formed for a second time later below. |
1390 | if(dg(istate) .and. .not. cv_pressure) then |
1391 | call assemble_divergence_matrix_cg(ct_m(istate)%ptr, state(istate), ct_rhs=ct_rhs(istate), & |
1392 | @@ -683,7 +692,11 @@ |
1393 | |
1394 | ! Set up the left C matrix in CMC |
1395 | |
1396 | - if(use_compressible_projection) then |
1397 | + if(compressible_eos) then |
1398 | + ! Note: If we are running a compressible multiphase simulation then the C^T matrix for each phase becomes: |
1399 | + ! rho*div(vfrac*u) for each incompressible phase |
1400 | + ! rho*div(vfrac*u) + vfrac*u*grad(rho) for the single compressible phase. |
1401 | + |
1402 | allocate(ctp_m(istate)%ptr) |
1403 | call allocate(ctp_m(istate)%ptr, ct_m(istate)%ptr%sparsity, (/1, u%dim/), name="CTP_m") |
1404 | ! NOTE that this is not optimal in that the ct_rhs |
1405 | @@ -691,7 +704,7 @@ |
1406 | if(cv_pressure) then |
1407 | call assemble_compressible_divergence_matrix_cv(ctp_m(istate)%ptr, state, ct_rhs(istate)) |
1408 | else |
1409 | - call assemble_compressible_divergence_matrix_cg(ctp_m(istate)%ptr, state, ct_rhs(istate)) |
1410 | + call assemble_compressible_divergence_matrix_cg(ctp_m(istate)%ptr, state, istate, ct_rhs(istate)) |
1411 | end if |
1412 | else |
1413 | ! Incompressible scenario |
1414 | @@ -708,7 +721,7 @@ |
1415 | end if |
1416 | end if |
1417 | |
1418 | - if (use_compressible_projection .or. cg_pressure_cv_test_continuity) then |
1419 | + if (compressible_eos .or. cg_pressure_cv_test_continuity) then |
1420 | if (have_rotated_bcs(u)) then |
1421 | if (dg(istate)) then |
1422 | call zero_non_owned(u) |
1423 | @@ -1041,7 +1054,7 @@ |
1424 | |
1425 | call profiler_toc(u, "assembly") |
1426 | |
1427 | - if(use_compressible_projection) then |
1428 | + if(compressible_eos) then |
1429 | call deallocate(ctp_m(istate)%ptr) |
1430 | deallocate(ctp_m(istate)%ptr) |
1431 | end if |
1432 | @@ -1252,15 +1265,12 @@ |
1433 | |
1434 | ewrite(1,*) 'Entering get_pressure_options' |
1435 | |
1436 | - |
1437 | ! Are we using a compressible projection? |
1438 | - use_compressible_projection = have_option(trim(p%option_path)//& |
1439 | - "/prognostic/scheme& |
1440 | - &/use_compressible_projection_method") |
1441 | + compressible_eos = option_count("/material_phase/equation_of_state/compressible") > 0 |
1442 | |
1443 | reassemble_all_cmc_m = have_option(trim(p%option_path)//& |
1444 | "/prognostic/scheme/update_discretised_equation") .or. & |
1445 | - use_compressible_projection |
1446 | + compressible_eos |
1447 | |
1448 | reassemble_all_ct_m = have_option(trim(p%option_path)//& |
1449 | "/prognostic/scheme/update_discretised_equation") |
1450 | @@ -1579,15 +1589,15 @@ |
1451 | call deallocate(kmk_rhs) |
1452 | |
1453 | cmc_m => extract_csr_matrix(state(istate), "PressurePoissonMatrix", stat) |
1454 | - |
1455 | - if(use_compressible_projection) then |
1456 | + |
1457 | + if(compressible_eos .and. have_option(trim(state(istate)%option_path)//'/equation_of_state/compressible')) then |
1458 | call allocate(compress_projec_rhs, p%mesh, "CompressibleProjectionRHS") |
1459 | |
1460 | if(cv_pressure) then |
1461 | call assemble_compressible_projection_cv(state, cmc_m, compress_projec_rhs, dt, & |
1462 | theta_pg, theta_divergence, reassemble_cmc_m) |
1463 | else |
1464 | - call assemble_compressible_projection_cg(state, cmc_m, compress_projec_rhs, dt, & |
1465 | + call assemble_compressible_projection_cg(state, istate, cmc_m, compress_projec_rhs, dt, & |
1466 | theta_pg, theta_divergence, reassemble_cmc_m) |
1467 | end if |
1468 | |
1469 | @@ -1720,7 +1730,7 @@ |
1470 | call addto(p, delta_p, scale=1.0/(theta_pg*dt)) |
1471 | ewrite_minmax(p) |
1472 | |
1473 | - if(use_compressible_projection) then |
1474 | + if(compressible_eos) then |
1475 | call update_compressible_density(state) |
1476 | end if |
1477 | |
1478 | |
1479 | === modified file 'assemble/Multiphase.F90' |
1480 | --- assemble/Multiphase.F90 2011-09-03 00:57:43 +0000 |
1481 | +++ assemble/Multiphase.F90 2012-08-09 00:46:20 +0000 |
1482 | @@ -46,7 +46,7 @@ |
1483 | private |
1484 | public :: get_phase_submaterials, get_nonlinear_volume_fraction, & |
1485 | calculate_diagnostic_phase_volume_fraction, & |
1486 | - add_fluid_particle_drag |
1487 | + add_fluid_particle_drag, add_heat_transfer |
1488 | |
1489 | contains |
1490 | |
1491 | @@ -268,7 +268,7 @@ |
1492 | integer :: ele |
1493 | type(element_type) :: test_function |
1494 | type(element_type), pointer :: u_shape |
1495 | - integer, dimension(:), pointer :: u_ele |
1496 | + integer, dimension(:), pointer :: u_nodes |
1497 | logical :: dg |
1498 | |
1499 | type(vector_field), pointer :: velocity_fluid |
1500 | @@ -281,14 +281,20 @@ |
1501 | integer :: i, dim |
1502 | logical :: not_found ! Error flag. Have we found the fluid phase? |
1503 | integer :: istate_fluid, istate_particle |
1504 | - |
1505 | + |
1506 | + ! Types of drag correlation |
1507 | + integer, parameter :: DRAG_CORRELATION_TYPE_STOKES = 1, DRAG_CORRELATION_TYPE_WEN_YU = 2, DRAG_CORRELATION_TYPE_ERGUN = 3 |
1508 | |
1509 | ewrite(1, *) "Entering add_fluid_particle_drag" |
1510 | + |
1511 | + ! Let's check whether we actually have at least one particle phase. |
1512 | + if(option_count("/material_phase/multiphase_properties/particle_diameter") == 0) then |
1513 | + FLExit("Fluid-particle drag enabled but no particle_diameter has been specified for the particle phase(s).") |
1514 | + end if |
1515 | |
1516 | ! Get the timestepping options |
1517 | call get_option("/timestepping/timestep", dt) |
1518 | - call get_option(trim(u%option_path)//"/prognostic/temporal_discretisation/theta", & |
1519 | - theta) |
1520 | + call get_option(trim(u%option_path)//"/prognostic/temporal_discretisation/theta", theta) |
1521 | |
1522 | ! For the big_m matrix. Controls whether the off diagonal entries are used |
1523 | block_mask = .false. |
1524 | @@ -300,15 +306,13 @@ |
1525 | dg = continuity(u) < 0 |
1526 | |
1527 | ! Is this phase a particle phase? |
1528 | - is_particle_phase = have_option("/material_phase["//int2str(istate-1)//& |
1529 | - &"]/multiphase_properties/particle_diameter") |
1530 | + is_particle_phase = have_option(trim(state(istate)%option_path)//"/multiphase_properties/particle_diameter") |
1531 | |
1532 | ! Retrieve the index of the fluid phase in the state array. |
1533 | not_found = .true. |
1534 | if(is_particle_phase) then |
1535 | do i = 1, size(state) |
1536 | - if(.not.have_option("/material_phase["//int2str(i-1)//& |
1537 | - &"]/multiphase_properties/particle_diameter")) then |
1538 | + if(.not.have_option(trim(state(i)%option_path)//"/multiphase_properties/particle_diameter")) then |
1539 | |
1540 | velocity_fluid => extract_vector_field(state(i), "Velocity") |
1541 | ! Aliased material_phases will also not have a particle_diameter, |
1542 | @@ -359,6 +363,8 @@ |
1543 | type(tensor_field), pointer :: viscosity_fluid |
1544 | type(scalar_field) :: nvfrac_fluid, nvfrac_particle |
1545 | real :: d ! Particle diameter |
1546 | + character(len=OPTION_PATH_LEN) :: drag_correlation_name |
1547 | + integer :: drag_correlation |
1548 | |
1549 | ! Get the necessary fields to calculate the drag force |
1550 | velocity_fluid => extract_vector_field(state(istate_fluid), "Velocity") |
1551 | @@ -371,8 +377,7 @@ |
1552 | density_particle => extract_scalar_field(state(istate_particle), "Density") |
1553 | viscosity_fluid => extract_tensor_field(state(istate_fluid), "Viscosity") |
1554 | |
1555 | - call get_option("/material_phase["//int2str(istate_particle-1)//& |
1556 | - &"]/multiphase_properties/particle_diameter", d) |
1557 | + call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter", d) |
1558 | |
1559 | ! Calculate the non-linear approximation to the PhaseVolumeFractions |
1560 | call allocate(nvfrac_fluid, vfrac_fluid%mesh, "NonlinearPhaseVolumeFraction") |
1561 | @@ -387,13 +392,25 @@ |
1562 | nu_particle => extract_vector_field(state(istate_particle), "NonlinearVelocity") |
1563 | oldu_fluid => extract_vector_field(state(istate_fluid), "OldVelocity") |
1564 | oldu_particle => extract_vector_field(state(istate_particle), "OldVelocity") |
1565 | - |
1566 | + |
1567 | + call get_option("/multiphase_interaction/fluid_particle_drag/drag_correlation/name", drag_correlation_name) |
1568 | + select case(trim(drag_correlation_name)) |
1569 | + case("stokes") |
1570 | + drag_correlation = DRAG_CORRELATION_TYPE_STOKES |
1571 | + case("wen_yu") |
1572 | + drag_correlation = DRAG_CORRELATION_TYPE_WEN_YU |
1573 | + case("ergun") |
1574 | + drag_correlation = DRAG_CORRELATION_TYPE_ERGUN |
1575 | + case("default") |
1576 | + FLAbort("Unknown correlation for fluid-particle drag") |
1577 | + end select |
1578 | + |
1579 | ! ----- Volume integrals over elements ------------- |
1580 | call profiler_tic(u, "element_loop") |
1581 | element_loop: do ele = 1, element_count(u) |
1582 | |
1583 | if(.not.dg .or. (dg .and. element_owned(u,ele))) then |
1584 | - u_ele=>ele_nodes(u, ele) |
1585 | + u_nodes => ele_nodes(u, ele) |
1586 | u_shape => ele_shape(u, ele) |
1587 | test_function = u_shape |
1588 | |
1589 | @@ -403,7 +420,7 @@ |
1590 | density_fluid, density_particle, & |
1591 | nu_fluid, nu_particle, & |
1592 | oldu_fluid, oldu_particle, & |
1593 | - viscosity_fluid, d) |
1594 | + viscosity_fluid, d, drag_correlation) |
1595 | end if |
1596 | |
1597 | end do element_loop |
1598 | @@ -421,7 +438,7 @@ |
1599 | density_fluid, density_particle, & |
1600 | nu_fluid, nu_particle, & |
1601 | oldu_fluid, oldu_particle, & |
1602 | - viscosity_fluid, d) |
1603 | + viscosity_fluid, d, drag_correlation) |
1604 | |
1605 | integer, intent(in) :: ele |
1606 | type(element_type), intent(in) :: test_function |
1607 | @@ -436,6 +453,7 @@ |
1608 | type(vector_field), intent(in) :: oldu_fluid, oldu_particle |
1609 | type(tensor_field), intent(in) :: viscosity_fluid |
1610 | real, intent(in) :: d ! Particle diameter |
1611 | + integer, intent(in) :: drag_correlation |
1612 | |
1613 | ! Local variables |
1614 | real, dimension(ele_ngi(u,ele)) :: vfrac_fluid_gi, vfrac_particle_gi |
1615 | @@ -447,14 +465,14 @@ |
1616 | |
1617 | real, dimension(ele_loc(u, ele), ele_ngi(u, ele), x%dim) :: du_t |
1618 | real, dimension(ele_ngi(u, ele)) :: detwei |
1619 | - real, dimension(u%dim, ele_loc(u,ele)) :: interaction_rhs_mat |
1620 | + real, dimension(u%dim, ele_loc(u,ele)) :: interaction_rhs |
1621 | real, dimension(ele_loc(u, ele), ele_loc(u, ele)) :: interaction_big_m_mat |
1622 | real, dimension(u%dim, ele_loc(u,ele)) :: rhs_addto |
1623 | real, dimension(u%dim, u%dim, ele_loc(u,ele), ele_loc(u,ele)) :: big_m_tensor_addto |
1624 | |
1625 | - real, dimension(ele_ngi(u,ele)) :: particle_re ! Particle Reynolds number |
1626 | - real, dimension(ele_ngi(u,ele)) :: drag_coefficient |
1627 | - real, dimension(ele_ngi(u,ele)) :: magnitude ! |v_f - v_p| |
1628 | + real, dimension(ele_ngi(u,ele)) :: particle_re_gi ! Particle Reynolds number |
1629 | + real, dimension(ele_ngi(u,ele)) :: drag_coefficient_gi |
1630 | + real, dimension(ele_ngi(u,ele)) :: magnitude_gi ! |v_f - v_p| |
1631 | |
1632 | real, dimension(ele_ngi(u,ele)) :: K |
1633 | real, dimension(ele_ngi(u,ele)) :: drag_force_big_m |
1634 | @@ -476,34 +494,55 @@ |
1635 | |
1636 | ! Compute the magnitude of the relative velocity |
1637 | do gi = 1, ele_ngi(u,ele) |
1638 | - magnitude(gi) = norm2(nu_fluid_gi(:,gi) - nu_particle_gi(:,gi)) |
1639 | + magnitude_gi(gi) = norm2(nu_fluid_gi(:,gi) - nu_particle_gi(:,gi)) |
1640 | end do |
1641 | |
1642 | ! Compute the particle Reynolds number |
1643 | ! (Assumes isotropic viscosity for now) |
1644 | - particle_re = (vfrac_fluid_gi*density_fluid_gi*magnitude*d) / viscosity_fluid_gi(1,1,:) |
1645 | + particle_re_gi = (vfrac_fluid_gi*density_fluid_gi*magnitude_gi*d) / viscosity_fluid_gi(1,1,:) |
1646 | |
1647 | ! Compute the drag coefficient |
1648 | - do gi = 1, ele_ngi(u,ele) |
1649 | - if(particle_re(gi) < 1000) then |
1650 | - drag_coefficient(gi) = (24.0/particle_re(gi)) |
1651 | - ! For the Wen & Yu (1966) drag term, this should be: |
1652 | - ! drag_coefficient(gi) = (24.0/particle_re(gi))*(1.0+0.15*particle_re(gi)**0.687) |
1653 | - else |
1654 | - drag_coefficient(gi) = 0.44 |
1655 | - end if |
1656 | - end do |
1657 | + select case(drag_correlation) |
1658 | + case(DRAG_CORRELATION_TYPE_STOKES) |
1659 | + ! Stokes drag correlation |
1660 | + do gi = 1, ele_ngi(u,ele) |
1661 | + if(particle_re_gi(gi) < 1000) then |
1662 | + drag_coefficient_gi(gi) = (24.0/particle_re_gi(gi)) |
1663 | + else |
1664 | + drag_coefficient_gi(gi) = 0.44 |
1665 | + end if |
1666 | + end do |
1667 | + |
1668 | + case(DRAG_CORRELATION_TYPE_WEN_YU) |
1669 | + ! Wen & Yu (1966) drag correlation |
1670 | + do gi = 1, ele_ngi(u,ele) |
1671 | + if(particle_re_gi(gi) < 1000) then |
1672 | + drag_coefficient_gi(gi) = (24.0/particle_re_gi(gi))*(1.0+0.15*particle_re_gi(gi)**0.687) |
1673 | + else |
1674 | + drag_coefficient_gi(gi) = 0.44 |
1675 | + end if |
1676 | + end do |
1677 | + |
1678 | + case(DRAG_CORRELATION_TYPE_ERGUN) |
1679 | + ! No drag coefficient is needed here. |
1680 | + end select |
1681 | |
1682 | - ! Don't let the drag_coefficient be NaN |
1683 | + ! Don't let the drag_coefficient_gi be NaN |
1684 | do gi = 1, ele_ngi(u,ele) |
1685 | - if(particle_re(gi) == 0) then |
1686 | - drag_coefficient(gi) = 1e16 |
1687 | + if(particle_re_gi(gi) == 0) then |
1688 | + drag_coefficient_gi(gi) = 1e16 |
1689 | end if |
1690 | end do |
1691 | |
1692 | - ! For the Wen & Yu (1966) drag term, this should be: |
1693 | - ! K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient*(vfrac_fluid_gi*density_fluid_gi*magnitude)/(d*vfrac_fluid_gi**2.7) |
1694 | - K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient*(vfrac_fluid_gi*density_fluid_gi*magnitude)/(d) |
1695 | + select case(drag_correlation) |
1696 | + case(DRAG_CORRELATION_TYPE_STOKES) |
1697 | + K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient_gi*(vfrac_fluid_gi*density_fluid_gi*magnitude_gi)/(d) |
1698 | + case(DRAG_CORRELATION_TYPE_WEN_YU) |
1699 | + ! Wen & Yu (1966) drag term |
1700 | + K = vfrac_particle_gi*(3.0/4.0)*drag_coefficient_gi*(vfrac_fluid_gi*density_fluid_gi*magnitude_gi)/(d*vfrac_fluid_gi**2.7) |
1701 | + case(DRAG_CORRELATION_TYPE_ERGUN) |
1702 | + K = 150.0*((vfrac_particle_gi**2)*viscosity_fluid_gi(1,1,:))/(vfrac_fluid_gi*(d**2)) + 1.75*(vfrac_particle_gi*density_fluid_gi*magnitude_gi/d) |
1703 | + end select |
1704 | |
1705 | if(is_particle_phase) then |
1706 | drag_force_big_m = -K |
1707 | @@ -519,7 +558,7 @@ |
1708 | |
1709 | ! Form the element interaction/drag matrix |
1710 | interaction_big_m_mat = shape_shape(test_function, u_shape, detwei*drag_force_big_m) |
1711 | - interaction_rhs_mat = shape_vector_rhs(test_function, drag_force_rhs, detwei) |
1712 | + interaction_rhs = shape_vector_rhs(test_function, drag_force_rhs, detwei) |
1713 | |
1714 | ! Add contribution |
1715 | big_m_tensor_addto = 0.0 |
1716 | @@ -531,16 +570,328 @@ |
1717 | end if |
1718 | do dim = 1, u%dim |
1719 | big_m_tensor_addto(dim, dim, :, :) = big_m_tensor_addto(dim, dim, :, :) - dt*theta*interaction_big_m_mat |
1720 | - rhs_addto(dim,:) = rhs_addto(dim,:) + matmul(interaction_big_m_mat, oldu_val(dim,:)) + interaction_rhs_mat(dim,:) |
1721 | + rhs_addto(dim,:) = rhs_addto(dim,:) + matmul(interaction_big_m_mat, oldu_val(dim,:)) + interaction_rhs(dim,:) |
1722 | end do |
1723 | |
1724 | ! Add the contribution to mom_rhs |
1725 | - call addto(mom_rhs, u_ele, rhs_addto) |
1726 | + call addto(mom_rhs, u_nodes, rhs_addto) |
1727 | ! Add to the big_m matrix |
1728 | - call addto(big_m, u_ele, u_ele, big_m_tensor_addto, block_mask=block_mask) |
1729 | + call addto(big_m, u_nodes, u_nodes, big_m_tensor_addto, block_mask=block_mask) |
1730 | |
1731 | end subroutine add_fluid_particle_drag_element |
1732 | |
1733 | end subroutine add_fluid_particle_drag |
1734 | |
1735 | + |
1736 | + !! Multiphase energy interaction term Q_i |
1737 | + !! to be added to the RHS of the internal energy equation. |
1738 | + subroutine add_heat_transfer(state, istate, internal_energy, matrix, rhs) |
1739 | + !!< This computes the inter-phase heat transfer term. |
1740 | + !!< Only between fluid and particle phase pairs. |
1741 | + !!< Uses the empirical correlation by Gunn (1978). |
1742 | + |
1743 | + type(state_type), dimension(:), intent(inout) :: state |
1744 | + integer, intent(in) :: istate |
1745 | + type(scalar_field), intent(in) :: internal_energy |
1746 | + type(csr_matrix), intent(inout) :: matrix |
1747 | + type(scalar_field), intent(inout) :: rhs |
1748 | + |
1749 | + ! Local variables |
1750 | + integer :: ele |
1751 | + type(element_type) :: test_function |
1752 | + type(element_type), pointer :: internal_energy_shape |
1753 | + integer, dimension(:), pointer :: internal_energy_nodes |
1754 | + logical :: dg |
1755 | + |
1756 | + type(vector_field), pointer :: x, velocity_fluid |
1757 | + |
1758 | + real :: dt, theta |
1759 | + |
1760 | + logical :: is_particle_phase |
1761 | + logical :: not_found ! Error flag. Have we found the fluid phase? |
1762 | + integer :: i, istate_fluid, istate_particle |
1763 | + |
1764 | + |
1765 | + ewrite(1, *) "Entering add_heat_transfer" |
1766 | + |
1767 | + ! Get the timestepping options |
1768 | + call get_option("/timestepping/timestep", dt) |
1769 | + call get_option(trim(internal_energy%option_path)//"/prognostic/temporal_discretisation/theta", & |
1770 | + theta) |
1771 | + |
1772 | + ! Get the coordinate field from state(istate) |
1773 | + x => extract_vector_field(state(istate), "Coordinate") |
1774 | + ewrite_minmax(x) |
1775 | + assert(x%dim == mesh_dim(internal_energy)) |
1776 | + assert(ele_count(x) == ele_count(internal_energy)) |
1777 | + |
1778 | + ! Are we using a discontinuous Galerkin discretisation? |
1779 | + dg = continuity(internal_energy) < 0 |
1780 | + |
1781 | + ! Is this phase a particle phase? |
1782 | + is_particle_phase = have_option(trim(state(istate)%option_path)//"/multiphase_properties/particle_diameter") |
1783 | + |
1784 | + ! Retrieve the index of the fluid phase in the state array. |
1785 | + not_found = .true. |
1786 | + if(is_particle_phase) then |
1787 | + do i = 1, size(state) |
1788 | + if(.not.have_option(trim(state(i)%option_path)//"/multiphase_properties/particle_diameter")) then |
1789 | + |
1790 | + velocity_fluid => extract_vector_field(state(i), "Velocity") |
1791 | + ! Aliased material_phases will also not have a particle_diameter, |
1792 | + ! so here we make sure that we don't count these as the fluid phase |
1793 | + if(.not.aliased(velocity_fluid)) then |
1794 | + istate_fluid = i |
1795 | + if(.not.not_found) then |
1796 | + FLExit("Heat transfer term does not currently support more than one fluid phase.") |
1797 | + end if |
1798 | + not_found = .false. |
1799 | + end if |
1800 | + |
1801 | + end if |
1802 | + end do |
1803 | + else |
1804 | + istate_fluid = istate |
1805 | + not_found = .false. |
1806 | + end if |
1807 | + |
1808 | + if(not_found) then |
1809 | + FLExit("No fluid phase found for the heat transfer term.") |
1810 | + end if |
1811 | + |
1812 | + ! If we have a fluid-particle pair, then assemble the heat transfer term |
1813 | + if(is_particle_phase) then |
1814 | + call assemble_heat_transfer(istate_fluid, istate) |
1815 | + else |
1816 | + state_loop: do i = 1, size(state) |
1817 | + if(i /= istate_fluid) then |
1818 | + call assemble_heat_transfer(istate_fluid, i) |
1819 | + end if |
1820 | + end do state_loop |
1821 | + end if |
1822 | + |
1823 | + ewrite(1, *) "Exiting add_heat_transfer" |
1824 | + |
1825 | + contains |
1826 | + |
1827 | + subroutine assemble_heat_transfer(istate_fluid, istate_particle) |
1828 | + |
1829 | + integer, intent(in) :: istate_fluid, istate_particle |
1830 | + |
1831 | + type(scalar_field), pointer :: vfrac_fluid, vfrac_particle |
1832 | + type(scalar_field), pointer :: density_fluid, density_particle |
1833 | + type(vector_field), pointer :: velocity_fluid, velocity_particle |
1834 | + type(scalar_field), pointer :: internal_energy_fluid, internal_energy_particle |
1835 | + type(scalar_field), pointer :: old_internal_energy_fluid, old_internal_energy_particle |
1836 | + type(vector_field), pointer :: nu_fluid, nu_particle ! Non-linear approximation to the Velocities |
1837 | + type(tensor_field), pointer :: viscosity_fluid |
1838 | + type(scalar_field) :: nvfrac_fluid, nvfrac_particle |
1839 | + real :: d ! Particle diameter |
1840 | + real :: k ! Effective gas conductivity |
1841 | + real :: C_fluid, C_particle ! Specific heat of the fluid and particle phases at constant volume |
1842 | + real :: gamma ! Ratio of specific heats for compressible phase |
1843 | + integer :: kstat, cstat_fluid, cstat_particle, gstat |
1844 | + |
1845 | + ! Get the necessary fields to calculate the heat transfer term |
1846 | + velocity_fluid => extract_vector_field(state(istate_fluid), "Velocity") |
1847 | + velocity_particle => extract_vector_field(state(istate_particle), "Velocity") |
1848 | + if(.not.aliased(velocity_particle)) then ! Don't count the aliased material_phases |
1849 | + |
1850 | + vfrac_fluid => extract_scalar_field(state(istate_fluid), "PhaseVolumeFraction") |
1851 | + vfrac_particle => extract_scalar_field(state(istate_particle), "PhaseVolumeFraction") |
1852 | + density_fluid => extract_scalar_field(state(istate_fluid), "Density") |
1853 | + density_particle => extract_scalar_field(state(istate_particle), "Density") |
1854 | + viscosity_fluid => extract_tensor_field(state(istate_fluid), "Viscosity") |
1855 | + |
1856 | + call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/particle_diameter", d) |
1857 | + |
1858 | + call get_option(trim(state(istate_fluid)%option_path)//"/multiphase_properties/effective_conductivity", k, kstat) |
1859 | + |
1860 | + call get_option(trim(state(istate_fluid)%option_path)//"/multiphase_properties/specific_heat", C_fluid, cstat_fluid) |
1861 | + |
1862 | + call get_option(trim(state(istate_particle)%option_path)//"/multiphase_properties/specific_heat", C_particle, cstat_particle) |
1863 | + |
1864 | + call get_option(trim(state(istate_fluid)%option_path)//"/equation_of_state/compressible/stiffened_gas/ratio_specific_heats", gamma, gstat) |
1865 | + |
1866 | + if(kstat /= 0) then |
1867 | + FLExit("For inter-phase heat transfer, an effective_conductivity needs to be specified for the fluid phase.") |
1868 | + end if |
1869 | + if(cstat_fluid /= 0 .or. cstat_particle /= 0) then |
1870 | + FLExit("For inter-phase heat transfer, a specific_heat needs to be specified for each phase.") |
1871 | + end if |
1872 | + if(gstat /= 0) then |
1873 | + FLExit("For inter-phase heat transfer, ratio_specific_heats needs to be specified for the compressible phase.") |
1874 | + end if |
1875 | + |
1876 | + ! Calculate the non-linear approximation to the PhaseVolumeFractions |
1877 | + call allocate(nvfrac_fluid, vfrac_fluid%mesh, "NonlinearPhaseVolumeFraction") |
1878 | + call allocate(nvfrac_particle, vfrac_particle%mesh, "NonlinearPhaseVolumeFraction") |
1879 | + call zero(nvfrac_fluid) |
1880 | + call zero(nvfrac_particle) |
1881 | + call get_nonlinear_volume_fraction(state(istate_fluid), nvfrac_fluid) |
1882 | + call get_nonlinear_volume_fraction(state(istate_particle), nvfrac_particle) |
1883 | + |
1884 | + ! Get the non-linear approximation to the Velocities |
1885 | + nu_fluid => extract_vector_field(state(istate_fluid), "NonlinearVelocity") |
1886 | + nu_particle => extract_vector_field(state(istate_particle), "NonlinearVelocity") |
1887 | + |
1888 | + ! Get the current and old internal energy fields |
1889 | + internal_energy_fluid => extract_scalar_field(state(istate_fluid), "InternalEnergy") |
1890 | + internal_energy_particle => extract_scalar_field(state(istate_particle), "InternalEnergy") |
1891 | + old_internal_energy_fluid => extract_scalar_field(state(istate_fluid), "OldInternalEnergy") |
1892 | + old_internal_energy_particle => extract_scalar_field(state(istate_particle), "OldInternalEnergy") |
1893 | + |
1894 | + ! ----- Volume integrals over elements ------------- |
1895 | + call profiler_tic(internal_energy, "element_loop") |
1896 | + element_loop: do ele = 1, element_count(internal_energy) |
1897 | + |
1898 | + if(.not.dg .or. (dg .and. element_owned(internal_energy,ele))) then |
1899 | + internal_energy_nodes => ele_nodes(internal_energy, ele) |
1900 | + internal_energy_shape => ele_shape(internal_energy, ele) |
1901 | + test_function = internal_energy_shape |
1902 | + |
1903 | + call add_heat_transfer_element(ele, test_function, internal_energy_shape, & |
1904 | + x, internal_energy, matrix, rhs, & |
1905 | + nvfrac_fluid, nvfrac_particle, & |
1906 | + density_fluid, density_particle, & |
1907 | + nu_fluid, nu_particle, & |
1908 | + internal_energy_fluid, & |
1909 | + internal_energy_particle, & |
1910 | + old_internal_energy_fluid, & |
1911 | + old_internal_energy_particle, & |
1912 | + viscosity_fluid, d, k, C_fluid, & |
1913 | + C_particle, gamma) |
1914 | + end if |
1915 | + |
1916 | + end do element_loop |
1917 | + call profiler_toc(internal_energy, "element_loop") |
1918 | + |
1919 | + call deallocate(nvfrac_fluid) |
1920 | + call deallocate(nvfrac_particle) |
1921 | + end if |
1922 | + |
1923 | + end subroutine assemble_heat_transfer |
1924 | + |
1925 | + subroutine add_heat_transfer_element(ele, test_function, internal_energy_shape, & |
1926 | + x, internal_energy, matrix, rhs, & |
1927 | + vfrac_fluid, vfrac_particle, & |
1928 | + density_fluid, density_particle, & |
1929 | + nu_fluid, nu_particle, & |
1930 | + internal_energy_fluid, & |
1931 | + internal_energy_particle, & |
1932 | + old_internal_energy_fluid, & |
1933 | + old_internal_energy_particle, & |
1934 | + viscosity_fluid, d, k, C_fluid, & |
1935 | + C_particle, gamma) |
1936 | + |
1937 | + integer, intent(in) :: ele |
1938 | + type(element_type), intent(in) :: test_function |
1939 | + type(element_type), intent(in) :: internal_energy_shape |
1940 | + type(vector_field), intent(in) :: x |
1941 | + type(scalar_field), intent(in) :: internal_energy |
1942 | + type(csr_matrix), intent(inout) :: matrix |
1943 | + type(scalar_field), intent(inout) :: rhs |
1944 | + |
1945 | + type(scalar_field), intent(in) :: vfrac_fluid, vfrac_particle |
1946 | + type(scalar_field), intent(in) :: density_fluid, density_particle |
1947 | + type(vector_field), intent(in) :: nu_fluid, nu_particle |
1948 | + type(scalar_field), intent(in) :: internal_energy_fluid, internal_energy_particle |
1949 | + type(scalar_field), intent(in) :: old_internal_energy_fluid, old_internal_energy_particle |
1950 | + type(tensor_field), intent(in) :: viscosity_fluid |
1951 | + real, intent(in) :: d, k, C_fluid, C_particle, gamma |
1952 | + |
1953 | + ! Local variables |
1954 | + real, dimension(ele_ngi(x,ele)) :: internal_energy_fluid_gi, internal_energy_particle_gi |
1955 | + real, dimension(ele_ngi(x,ele)) :: vfrac_fluid_gi, vfrac_particle_gi |
1956 | + real, dimension(ele_ngi(x,ele)) :: density_fluid_gi, density_particle_gi |
1957 | + real, dimension(x%dim, ele_ngi(x,ele)) :: nu_fluid_gi, nu_particle_gi |
1958 | + real, dimension(x%dim, x%dim, ele_ngi(x,ele)) :: viscosity_fluid_gi |
1959 | + |
1960 | + real, dimension(ele_loc(internal_energy,ele)) :: old_internal_energy_val |
1961 | + |
1962 | + real, dimension(ele_ngi(x,ele)) :: detwei |
1963 | + real, dimension(ele_ngi(x,ele)) :: coefficient_for_matrix, coefficient_for_rhs |
1964 | + real, dimension(ele_loc(internal_energy,ele)) :: rhs_addto |
1965 | + real, dimension(ele_loc(internal_energy,ele), ele_loc(internal_energy,ele)) :: matrix_addto |
1966 | + |
1967 | + real, dimension(ele_ngi(x,ele)) :: particle_Re ! Particle Reynolds number |
1968 | + real, dimension(ele_ngi(x,ele)) :: Pr ! Prandtl number |
1969 | + real, dimension(ele_ngi(x,ele)) :: particle_Nu ! Particle Nusselt number |
1970 | + real, dimension(ele_ngi(x,ele)) :: velocity_magnitude ! |v_f - v_p| |
1971 | + |
1972 | + real, dimension(ele_ngi(x,ele)) :: Q ! heat transfer term = Q*(T_p - T_f) |
1973 | + real, dimension(ele_loc(internal_energy,ele), ele_loc(internal_energy,ele)) :: heat_transfer_matrix |
1974 | + real, dimension(ele_loc(internal_energy,ele)) :: heat_transfer_rhs |
1975 | + |
1976 | + integer :: gi |
1977 | + |
1978 | + ! Compute detwei |
1979 | + call transform_to_physical(x, ele, detwei) |
1980 | + |
1981 | + ! Get the values of the necessary fields at the Gauss points |
1982 | + vfrac_fluid_gi = ele_val_at_quad(vfrac_fluid, ele) |
1983 | + vfrac_particle_gi = ele_val_at_quad(vfrac_particle, ele) |
1984 | + density_fluid_gi = ele_val_at_quad(density_fluid, ele) |
1985 | + density_particle_gi = ele_val_at_quad(density_particle, ele) |
1986 | + nu_fluid_gi = ele_val_at_quad(nu_fluid, ele) |
1987 | + nu_particle_gi = ele_val_at_quad(nu_particle, ele) |
1988 | + viscosity_fluid_gi = ele_val_at_quad(viscosity_fluid, ele) |
1989 | + internal_energy_fluid_gi = ele_val_at_quad(internal_energy_fluid, ele) |
1990 | + internal_energy_particle_gi = ele_val_at_quad(internal_energy_particle, ele) |
1991 | + |
1992 | + ! Compute the magnitude of the relative velocity |
1993 | + do gi = 1, ele_ngi(x,ele) |
1994 | + velocity_magnitude(gi) = norm2(nu_fluid_gi(:,gi) - nu_particle_gi(:,gi)) |
1995 | + end do |
1996 | + |
1997 | + ! Compute the particle Reynolds number |
1998 | + ! (Assumes isotropic viscosity for now) |
1999 | + particle_Re = (density_fluid_gi*velocity_magnitude*d) / viscosity_fluid_gi(1,1,:) |
2000 | + |
2001 | + ! Compute the Prandtl number |
2002 | + ! (Assumes isotropic viscosity for now) |
2003 | + ! Note: C_fluid (at constant volume) multiplied by gamma = C_fluid at constant pressure |
2004 | + Pr = C_fluid*gamma*viscosity_fluid_gi(1,1,:)/k |
2005 | + |
2006 | + particle_Nu = (7.0 - 10.0*vfrac_fluid_gi + 5.0*vfrac_fluid_gi**2)*(1.0 + 0.7*(particle_Re**0.2)*(Pr**(1.0/3.0))) + & |
2007 | + (1.33 - 2.4*vfrac_fluid_gi + 1.2*vfrac_fluid_gi**2)*(particle_Re**0.7)*(Pr**(1.0/3.0)) |
2008 | + |
2009 | + Q = (6.0*k*vfrac_particle_gi*particle_Nu)/(d**2) |
2010 | + |
2011 | + ! Note that the transfer term is defined in terms of temperatures (T_fluid and T_particle) |
2012 | + ! Let's convert the temperatures to internal energy (E) using E_i = C_i*T_i, |
2013 | + ! where C is the specific heat of phase i at constant volume. |
2014 | + if(is_particle_phase) then |
2015 | + coefficient_for_matrix = -Q/C_particle |
2016 | + coefficient_for_rhs = -Q*(-internal_energy_fluid_gi/C_fluid) |
2017 | + else |
2018 | + coefficient_for_matrix = -Q/C_fluid |
2019 | + coefficient_for_rhs = Q*(internal_energy_particle_gi/C_particle) |
2020 | + end if |
2021 | + |
2022 | + ! Form the element heat transfer matrix and RHS |
2023 | + heat_transfer_matrix = shape_shape(test_function, internal_energy_shape, detwei*coefficient_for_matrix) |
2024 | + heat_transfer_rhs = shape_rhs(test_function, coefficient_for_rhs*detwei) |
2025 | + |
2026 | + ! Add contribution |
2027 | + matrix_addto = 0.0 |
2028 | + rhs_addto = 0.0 |
2029 | + if(is_particle_phase) then |
2030 | + old_internal_energy_val = ele_val(old_internal_energy_particle, ele) |
2031 | + else |
2032 | + old_internal_energy_val = ele_val(old_internal_energy_fluid, ele) |
2033 | + end if |
2034 | + |
2035 | + matrix_addto = matrix_addto - dt*theta*heat_transfer_matrix |
2036 | + rhs_addto = rhs_addto + matmul(heat_transfer_matrix, old_internal_energy_val) + heat_transfer_rhs |
2037 | + |
2038 | + ! Add to the internal energy equation's RHS |
2039 | + call addto(rhs, internal_energy_nodes, rhs_addto) |
2040 | + ! Add to the matrix |
2041 | + call addto(matrix, internal_energy_nodes, internal_energy_nodes, matrix_addto) |
2042 | + |
2043 | + end subroutine add_heat_transfer_element |
2044 | + |
2045 | + end subroutine add_heat_transfer |
2046 | + |
2047 | end module multiphase_module |
2048 | |
2049 | === modified file 'main/Fluids.F90' |
2050 | --- main/Fluids.F90 2012-05-29 13:43:29 +0000 |
2051 | +++ main/Fluids.F90 2012-08-09 00:46:20 +0000 |
2052 | @@ -699,7 +699,7 @@ |
2053 | else if(have_option(trim(field_optionpath_list(it)) // & |
2054 | & "/prognostic/spatial_discretisation/continuous_galerkin")) then |
2055 | |
2056 | - call solve_field_equation_cg(field_name_list(it), state(field_state_list(it)), dt) |
2057 | + call solve_field_equation_cg(field_name_list(it), state, field_state_list(it), dt) |
2058 | else |
2059 | |
2060 | ewrite(2, *) "Not solving scalar field " // trim(field_name_list(it)) // " in state " // trim(state(field_state_list(it))%name) //" in an advdif-like subroutine." |
2061 | |
2062 | === modified file 'main/Shallow_Water.F90' |
2063 | --- main/Shallow_Water.F90 2011-07-22 12:39:08 +0000 |
2064 | +++ main/Shallow_Water.F90 2012-08-09 00:46:20 +0000 |
2065 | @@ -241,7 +241,7 @@ |
2066 | call project_cartesian_to_local(state(1), v_field) |
2067 | end if |
2068 | |
2069 | - call execute_timestep(state(1), dt) |
2070 | + call execute_timestep(state, 1, dt) |
2071 | |
2072 | call set_prescribed_field_values(state, exclude_interpolated=.true., & |
2073 | exclude_nonreprescribed=.true., time=current_time + dt) |
2074 | @@ -528,9 +528,10 @@ |
2075 | |
2076 | end subroutine insert_time_in_state |
2077 | |
2078 | - subroutine execute_timestep(state, dt) |
2079 | + subroutine execute_timestep(state, istate, dt) |
2080 | implicit none |
2081 | - type(state_type), intent(inout) :: state |
2082 | + type(state_type), dimension(:), intent(inout) :: state |
2083 | + integer, intent(in) :: istate |
2084 | real, intent(in) :: dt |
2085 | |
2086 | !! Layer thickness |
2087 | @@ -548,20 +549,20 @@ |
2088 | real :: energy |
2089 | logical :: have_source |
2090 | |
2091 | - !Pull the fields out of state |
2092 | - D=>extract_scalar_field(state, "LayerThickness") |
2093 | - U=>extract_vector_field(state, "LocalVelocity") |
2094 | + !Pull the fields out of state(istate) |
2095 | + D=>extract_scalar_field(state(istate), "LayerThickness") |
2096 | + U=>extract_vector_field(state(istate), "LocalVelocity") |
2097 | have_source = .false. |
2098 | - if (has_vector_field(state, "LocalVelocitySource")) then |
2099 | + if (has_vector_field(state(istate), "LocalVelocitySource")) then |
2100 | have_source = .true. |
2101 | - source => extract_vector_field(state, "LocalVelocitySource") |
2102 | + source => extract_vector_field(state(istate), "LocalVelocitySource") |
2103 | end if |
2104 | dim = U%dim |
2105 | |
2106 | call execute_timestep_setup(D,U,d_rhs,u_rhs,advecting_u, & |
2107 | old_u,old_d,delta_d,delta_u) |
2108 | |
2109 | - call insert(state,advecting_u,"NonlinearVelocity") |
2110 | + call insert(state(istate),advecting_u,"NonlinearVelocity") |
2111 | |
2112 | call set(advecting_u,u) |
2113 | call set(old_u,u) |
2114 | @@ -577,20 +578,20 @@ |
2115 | do d1 = 1, dim |
2116 | velocity_cpt = extract_scalar_field(u,d1) |
2117 | old_velocity_cpt = extract_scalar_field(old_u,d1) |
2118 | - call insert(state,velocity_cpt,"VelocityComponent") |
2119 | - call insert(state,old_velocity_cpt,"VelocityComponentOld") |
2120 | - call solve_advection_diffusion_dg("VelocityComponent", state) |
2121 | + call insert(state(istate),velocity_cpt,"VelocityComponent") |
2122 | + call insert(state(istate),old_velocity_cpt,"VelocityComponentOld") |
2123 | + call solve_advection_diffusion_dg("VelocityComponent", state(istate)) |
2124 | end do |
2125 | end if |
2126 | !pressure advection |
2127 | if(.not.exclude_pressure_advection) then |
2128 | - call solve_field_equation_cg("LayerThickness", state, dt, & |
2129 | + call solve_field_equation_cg("LayerThickness", state, 1, dt, & |
2130 | "NonlinearVelocity") |
2131 | end if |
2132 | |
2133 | if(hybridized) then |
2134 | call solve_hybridized_helmholtz(& |
2135 | - &state,& |
2136 | + &state(istate),& |
2137 | &U_out=U_rhs,D_out=D_rhs,& |
2138 | &compute_cartesian=.true.,& |
2139 | &check_continuity=.true.,output_dense=.false.) |
2140 | @@ -623,8 +624,8 @@ |
2141 | !Construct explicit parts of h rhs in wave equation |
2142 | call get_d_rhs(d_rhs,u_rhs,D,U,div_mat,big_mat,D0,dt,theta) |
2143 | |
2144 | - if (has_scalar_field(state, "LayerThicknessSource")) then |
2145 | - d_src => extract_scalar_field(state, "LayerThicknessSource") |
2146 | + if (has_scalar_field(state(istate), "LayerThicknessSource")) then |
2147 | + d_src => extract_scalar_field(state(istate), "LayerThicknessSource") |
2148 | call allocate(md_src, d_src%mesh, "MassMatrixTimesLayerThicknessSource") |
2149 | call zero(md_src) |
2150 | call mult(md_src, h_mass_mat, d_src) |
2151 | |
2152 | === modified file 'manual/bibliography.bib' |
2153 | --- manual/bibliography.bib 2011-11-25 12:32:35 +0000 |
2154 | +++ manual/bibliography.bib 2012-08-09 00:46:20 +0000 |
2155 | @@ -912,6 +912,16 @@ |
2156 | year = {2005} |
2157 | } |
2158 | |
2159 | +@Article{ ergun1952, |
2160 | + author = {Ergun, S.}, |
2161 | + journal = {Chemical Engineering Progress}, |
2162 | + number = {2}, |
2163 | + pages = {89-94}, |
2164 | + title = {Fluid flow through packed columns}, |
2165 | + volume = {48}, |
2166 | + year = {1952} |
2167 | +} |
2168 | + |
2169 | @Article{ erturk2005, |
2170 | author = {E. Erturk and T. C. Corke and C. Gokcol}, |
2171 | journal = {International Journal for Numerical Methods in Fluids}, |
2172 | @@ -1298,6 +1308,16 @@ |
2173 | year = {1994} |
2174 | } |
2175 | |
2176 | +@Article{ gunn1978, |
2177 | + author = {Gunn, D. J.}, |
2178 | + journal = {International Journal of Heat and Mass Transfer}, |
2179 | + number = {4}, |
2180 | + pages = {467-476}, |
2181 | + title = {Transfer of heat or mass to particles in fixed and fluidised beds}, |
2182 | + volume = {21}, |
2183 | + year = {1978} |
2184 | +} |
2185 | + |
2186 | @Book{ haidvogelbook, |
2187 | author = {D. B. Haidvogel and A. Beckmann}, |
2188 | publisher = {Imperial {C}ollege {P}ress}, |
2189 | @@ -2877,6 +2897,15 @@ |
2190 | year = 2008 |
2191 | } |
2192 | |
2193 | +@Article{ wen_yu_1966, |
2194 | + title = {Mechanics of Fluidization}, |
2195 | + author = {C.~Y. Wen and Y.~H. Yu}, |
2196 | + journal = {Chem. Eng. Prog. Symp. Ser.}, |
2197 | + volume = 62, |
2198 | + number = 100, |
2199 | + year = 1966 |
2200 | +} |
2201 | + |
2202 | @Book{ wilcox1998turbulence, |
2203 | author = {Wilcox, D.C.}, |
2204 | publisher = {DCW industries La Canada, CA}, |
2205 | |
2206 | === modified file 'manual/configuring_fluidity.tex' |
2207 | --- manual/configuring_fluidity.tex 2012-07-26 15:38:42 +0000 |
2208 | +++ manual/configuring_fluidity.tex 2012-08-09 00:46:20 +0000 |
2209 | @@ -719,7 +719,7 @@ |
2210 | |
2211 | When configuring ocean problems (or single material/single phase fluids |
2212 | problems), only one \option{/material\_phase} is required. Multi-material |
2213 | -and multi-phase problems will require one \option{/material\_phase} for each |
2214 | +and multiphase problems will require one \option{/material\_phase} for each |
2215 | phase or material in the problem. |
2216 | |
2217 | Note that you must give each of your material phases a name. |
2218 | @@ -1184,7 +1184,8 @@ |
2219 | \index{conjugate gradient} |
2220 | \index{multigrid} |
2221 | \index{PETSc} |
2222 | -As described in Section \ref{ND_Linear_solvers}, for the solution of large sparse linear systems, the so called iterative methods are usually employed. These methods avoid having to explicitly construct the inverse of the matrix, which is generally dense and therefore costly to compute (both in memory and computer time). FLUIDITY is linked to PETSc: a suite of data structures and routines for the scalable (parallel) solution of scientific applications modelled by partial differential equations. It employs the MPI standard for parallelism. FLUIDITY therefore supports any iterative method provided by the PETSc library (http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/KSP/KSPType.html --- available methods may depend on the PETSc library installed on your system). Examples include Conjugate Gradient (CG), GMRES and FGMRES (Flexible GMRES). Some options are explicitly listed under \option{solver/iterative\_method}, for example CG: \option{solver/iterative\_method::cg}, whereas others can be selected entering the name of the chosen scheme in \option{solver/iterative\_method}. |
2223 | +As described in Section \ref{ND_Linear_solvers}, for the solution of large sparse linear systems, the so called iterative methods are usually employed. These methods avoid having to explicitly construct the inverse of the matrix, which is generally dense and therefore costly to compute (both in memory and computer time). FLUIDITY is linked to PETSc: a suite of data structures and routines for the scalable (parallel) solution of scientific applications modelled by partial differential equations. It employs the MPI standard for parallelism. FLUIDITY therefore supports any iterative method provided by the PETSc library (http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/KSP/KSPType.html --- available methods may depend on the PETSc library installed on your system). Examples include Conjugate Gradient (CG), GMRES and FGMRES (Flexible GMRES). Some options are explicitly listed under \option{solver/iterative\_method}, for example CG: \option{solver/iterative\_method::cg}, whereas |
2224 | +others can be selected entering the name of the chosen scheme in \option{solver/iterative\_method}. |
2225 | |
2226 | \subsection{Preconditioner} |
2227 | \index{linear solvers!preconditioners} |
2228 | @@ -2313,7 +2314,7 @@ |
2229 | Models with one prognostic velocity field per \option{material\_phase} are referred to as multiple phase models. The use of these multiple velocity fields permits the inter-penetration and interaction between different phases. |
2230 | |
2231 | \subsubsection{Simulation requirements} |
2232 | -In a multi-phase simulation, each \option{material\_phase} requires: |
2233 | +In a multiphase simulation, each \option{material\_phase} requires: |
2234 | \begin{itemize} |
2235 | \item an equation of state, and |
2236 | \item a PhaseVolumeFraction scalar field. |
2237 | @@ -2328,13 +2329,15 @@ |
2238 | \alpha_N = 1 - \sum_{i=1}^{N-1}\alpha_i\textrm{.} |
2239 | \end{equation} |
2240 | |
2241 | -\subsubsection{Inter-phase interactions} |
2242 | -Inter-phase interactions can be included under \option{../multiphase\_properties} in Diamond. Currently, \fluidity\ only supports fluid-particle drag between the continuous phase and dispersed phase(s), given by: |
2243 | -\begin{equation}\label{eq:multiphase_drag_term} |
2244 | +\subsubsection{Inter-phase momentum interaction} |
2245 | +Currently, \fluidity\ only supports fluid-particle drag between the continuous phase and dispersed phase(s). This can be enabled by switching on the \option{/multiphase\_interaction/fluid\_particle\_drag} option in Diamond and specifying the drag correlation: Stokes, \cite{wen_yu_1966} or \cite{ergun1952}. |
2246 | + |
2247 | +The drag force using the Stokes drag correlation is given by: |
2248 | +\begin{equation}\label{eq:stokes_drag_force} |
2249 | \mathbf{F}_D = \frac{3\ \alpha_p\ C_D\ \alpha_f\ \rho_f\ |\mathbf{u}_f-\mathbf{u}_p|\ (\mathbf{u}_f-\mathbf{u}_p)}{4\ d}, |
2250 | \end{equation} |
2251 | where $f$ and $p$ denote the fluid (i.e. continuous) and particle (i.e. dispersed) phases respectively, and $d$ is the diameter of a single particle in the dispersed phase. The drag coefficient $C_D$ is defined as: |
2252 | -\begin{equation}\label{eq:drag_coefficient} |
2253 | +\begin{equation}\label{eq:stokes_drag_coefficient} |
2254 | C_D = \frac{24}{\mathrm{Re}}, |
2255 | \end{equation} |
2256 | with |
2257 | @@ -2343,17 +2346,52 @@ |
2258 | \end{equation} |
2259 | Note that $\mu_f$ denotes the isotropic viscosity of the fluid (i.e. continuous) phase. |
2260 | |
2261 | -Within a dispersed phase, entering a value for $d$ in the \option{../multiphase\_properties/particle\_diameter} option enables fluid-particle interaction between itself and the continuous phase. |
2262 | +With the drag correlation by \cite{wen_yu_1966}, $\mathbf{F}_D$ and $C_D$ become: |
2263 | +\begin{equation}\label{eq:wen_yu_drag_force} |
2264 | +\mathbf{F}_D = \frac{3\ \alpha_p\ C_D\ \alpha_f\ \rho_f\ |\mathbf{u}_f-\mathbf{u}_p|\ (\mathbf{u}_f-\mathbf{u}_p)}{4\ d\alpha_f^{2.7}}, |
2265 | +\end{equation} |
2266 | +\begin{equation}\label{eq:wen_yu_drag_coefficient} |
2267 | +C_D = \frac{24}{\mathrm{Re}}\left(1.0 + 0.15\mathrm{Re}^{0.687}\right). |
2268 | +\end{equation} |
2269 | +In contrast to the Stokes drag correlation, the \cite{wen_yu_1966} drag correlation is more suitable for larger particle Reynolds number flows. |
2270 | + |
2271 | +For dense multiphase flows with $\alpha_p > 0.2$, the drag correlation by \cite{ergun1952} is often the most appropriate: |
2272 | +\begin{equation}\label{eq:ergun_drag_force} |
2273 | +\mathbf{F}_D = \left(150\frac{\alpha_p^2\mu_f}{\alpha_f d_p^2} + 1.75\frac{\alpha_p\rho_f|\mathbf{u}_f-\mathbf{u}_p|}{d_p}\right)\left(\mathbf{u}_f-\mathbf{u}_p\right). |
2274 | +\end{equation} |
2275 | + |
2276 | +Note that within each dispersed phase, a value for $d$ must be specified in \option{../multiphase\_properties/particle\_diameter} regardless of which drag correlation is used. |
2277 | + |
2278 | +\subsubsection{Inter-phase energy interaction} |
2279 | +This subsection considers compressible multiphase flow simulations where each phase has its own InternalEnergy field. Users can choose to include the inter-phase heat transfer term by \cite{gunn1978}, denoted $Q$ in Section \ref{sec:multiphase_equations}, on the RHS of each internal energy equation: |
2280 | + |
2281 | +\begin{equation}\label{eq:gunn_heat_transfer_term} |
2282 | +Q = \frac{6 k \alpha_p \mathrm{Nu}_p}{d_p^2}\left(\frac{e_f}{C_f} - \frac{e_p}{C_p}\right), |
2283 | +\end{equation} |
2284 | +where |
2285 | +\begin{equation} |
2286 | +\mathrm{Nu}_p = \left(7 - 10\alpha_f + 5\alpha_f^2\right)\left(1 + 0.7\mathrm{Re}_p^{0.2}\mathrm{Pr}^{\frac{1}{3}}\right) + \left(1.33 - 2.4\alpha_f + 1.2\alpha_f^2\right)\mathrm{Re}_p^{0.7}\mathrm{Pr}^{\frac{1}{3}}, |
2287 | +\end{equation} |
2288 | +\begin{equation} |
2289 | + \mathrm{Pr} = \frac{C_f \gamma \mu_f}{k}, |
2290 | +\end{equation} |
2291 | +and |
2292 | +\begin{equation} |
2293 | + \mathrm{Re} = \frac{\rho_f d_p |\mathbf{u}_f-\mathbf{u}_p|}{\mu_f}, |
2294 | +\end{equation} |
2295 | +are the Nusselt, Prandtl and Reynolds number respectively. $C_i$ denotes the specific heat of phase i at constant volume, and $k$ denotes the effective conductivity of the fluid phase; these must be specified in \option{../multiphase\_properties/effective\_conductivity} and \option{../multiphase\_properties/specific\_heat}. Note that we have written the above in terms of internal energy (rather than temperature) by dividing $e_i$ by the specific heat at constant volume. |
2296 | + |
2297 | +To include this term, switch on the \option{/multiphase\_interaction/heat\_transfer} option in Diamond. |
2298 | |
2299 | \subsubsection{Current limitations} |
2300 | \begin{itemize} |
2301 | - \item Boussinesq and fully compressible multi-phase simulations are not yet supported. |
2302 | + \item Boussinesq multiphase simulations are not yet supported. |
2303 | \item The momentum equation for each \option{material\_phase} can only be discretised in non-conservative form. |
2304 | \item The stress term must be in tensor form, and only works with an isotropic viscosity. |
2305 | - \item \option{bassi\_rebay} and \option{compact\_discontinuous\_galerkin} are currently the only DG viscosity schemes available for multi-phase flow simulations. |
2306 | + \item \option{bassi\_rebay} and \option{compact\_discontinuous\_galerkin} are currently the only DG viscosity schemes available for multiphase flow simulations. |
2307 | \item Discontinuous \option{PhaseVolumeFraction} fields are not yet supported. |
2308 | - \item The Pressure field only supports \option{continuous\_galerkin} and \option{control\_volume} discretisations. |
2309 | - \item Prescribed velocity fields cannot yet be used in multi-phase simulations. |
2310 | + \item For compressible multiphase flow simulations, the Pressure field only supports a \option{continuous\_galerkin} discretisation. |
2311 | + \item Prescribed velocity fields cannot yet be used in multiphase simulations. |
2312 | \item Fluid-particle drag can currently only support one continuous (i.e. fluid) phase. |
2313 | \end{itemize} |
2314 | |
2315 | @@ -2396,7 +2434,8 @@ |
2316 | |
2317 | This section describes how to configure \fluidity to simulate single phase incompressible Darcy flow. |
2318 | |
2319 | -First, the \option{Porosity} and \option{Permeability} fields found under the option element \option{/porous\_media} require configuring. Both fields can be either \option{prescribed} or \option{diagnostic}. The \option{Permeability} field can be either a \option{scalar\_field} or a \option{vector\_field}. The \option{Porosity} field will only be used in the transport of scalar fields and the metric tensor. The \option{Permeability} field will only be used in forming the absorption coefficient in the Darcy velocity equation. The use of these two fields should be considered when selecting their associated mesh. If \option{Porosity} is to be included in the transport of scalar fields (and the metric tensor) that are discretised using control volumes, due to the way this method is coded, the mesh associated with the \option{Porosity} must be either element wise (meaning discontinuous order zero) or the same order as field to be transported. It is therefore recommended that the \option{Porosity} field always be associated with a mesh that is discontinuous order zero. It is recommended that the \option{Permeability} field also be associated with a mesh that is discontinuous order zero, due to the recommended Darcy velocity - pressure element pair (being p0p1cv or p0p1 with the continuity tested with the p1 control volume dual space) |
2320 | +First, the \option{Porosity} and \option{Permeability} fields found under the option element \option{/porous\_media} require configuring. Both fields can be either \option{prescribed} or \option{diagnostic}. The \option{Permeability} field can be either a \option{scalar\_field} or a \option{vector\_field}. The \option{Porosity} field will only be used in the transport of scalar fields and the metric tensor. The \option{Permeability} field will only be used in forming the absorption coefficient in the Darcy velocity equation. The use of these two fields should be considered when selecting their associated mesh. If \option{Porosity} is to be included in the transport of scalar fields (and the metric tensor) that are discretised using control volumes, due to the way this method is coded, the mesh associated with the \option{Porosity} must be either element wise (meaning discontinuous order zero) or the same order as field to be transported. It is therefore recommended that the \option{Porosity} field always be |
2321 | +associated with a mesh that is discontinuous order zero. It is recommended that the \option{Permeability} field also be associated with a mesh that is discontinuous order zero, due to the recommended Darcy velocity - pressure element pair (being p0p1cv or p0p1 with the continuity tested with the p1 control volume dual space) |
2322 | |
2323 | Second, the prognostic \option{Velocity} vector field of the (only) phase is now taken to be the Darcy velocity \ref{sec:porous_media_darcy_flow_equations}. There is no option to select Darcy flow but this can be achieved via a careful selection of the \option{Velocity} field options tree given by: |
2324 | \begin{itemize} |
2325 | |
2326 | === modified file 'manual/model_equations.tex' |
2327 | --- manual/model_equations.tex 2012-07-26 15:38:42 +0000 |
2328 | +++ manual/model_equations.tex 2012-08-09 00:46:20 +0000 |
2329 | @@ -727,7 +727,7 @@ |
2330 | |
2331 | \subsection{Multi-material simulations} |
2332 | \index{multi-material flow} |
2333 | -The ability to differentiate between regions with distinct material properties is of fundamental importance in the modelling of many physical systems. Two different approaches exist for achieving this: the multi-material approach and the multi-phase approach. The multi-material approach is implemented within \fluidity, and the multi-phase approach (discussed in the next subsection) is currently under development. |
2334 | +The ability to differentiate between regions with distinct material properties is of fundamental importance in the modelling of many physical systems. Two different approaches exist for achieving this: the multi-material approach and the multiphase approach. The multi-material approach is implemented within \fluidity, and the multiphase approach (discussed in the next subsection) is currently under development. |
2335 | In situations where the model can resolve physical mixing of immiscible materials, or where there is no mixing, only one velocity field (and hence one momentum equation) is required to describe the flow of all materials. The \emph{multi-material} approach, considers all materials to be immiscible materials separated by a sharp interface. |
2336 | |
2337 | In a multi-material approach, the various forms of the conservation equations can be solved for multiple material flows if the point-wise mass density $\rho$ in the equations is defined as the bulk density at each point. If the flow comprises $n$ materials and the volume fraction of the $i^{th}$ material is denoted $\phi_i$ then the bulk density is given by: |
2338 | @@ -749,14 +749,17 @@ |
2339 | \end{equation} |
2340 | where the volume fraction field at time zero must be specified. |
2341 | |
2342 | -\subsection{Multi-phase simulations} |
2343 | +\subsection{Multiphase simulations} |
2344 | \label{sec:multiphase_equations} |
2345 | -\index{multi-phase flow} |
2346 | -Multi-phase flows are defined by \cite{prosperettiEtAl2007} to be flows in which two or more phases of matter (solid, liquid, gas, etc) are simultaneously present and are allowed to inter-penetrate. Simple examples include the flow of a fizzy drink which is composed of a liquid and a finite number of gas bubbles, the transportation of solid sediment particles in a river, and the flow of blood cells around the human body. |
2347 | +\index{multiphase flow} |
2348 | +Multiphase flows are defined by \cite{prosperettiEtAl2007} to be flows in which two or more phases of matter (solid, liquid, gas, etc) are simultaneously present and are allowed to inter-penetrate. Simple examples include the flow of a fizzy drink which is composed of a liquid and a finite number of gas bubbles, the transportation of solid sediment particles in a river, and the flow of blood cells around the human body. |
2349 | |
2350 | Further to the above definition, each phase is classed as either \textit{continuous} or \textit{dispersed}, where a continuous phase is a connected liquid or gas substance in which dispersed phases (comprising a finite number of solid particles, liquid droplets and/or gas bubbles) may be immersed \citep{croweEtAl1998}. |
2351 | |
2352 | -To enable the mixing and inter-penetration of phases, a separate velocity field (and hence a separate momentum equation) is assigned to each one and solved for. Extra terms are then included to account for inter-phase interactions. Furthermore, the model currently assumes no mass transfer between phases, incompressible flow, and a common pressure field $p$ so that only one continuity equation is used. Thus, the continuity equation and momentum equation for phase $i$ (based on the derivation in \cite{ishii1975}, written in non-conservative form) are: |
2353 | +To enable the mixing and inter-penetration of phases, a separate velocity field (and hence a separate momentum equation) is assigned to each one and solved for. Extra terms are then included to account for inter-phase interactions. Furthermore, the model currently assumes no mass transfer between phases, incompressible flow, and a common pressure field $p$ so that only one continuity equation is used. The form of the governing equations depends on whether we are dealing with incompressible or compressible flow. |
2354 | + |
2355 | +\subsubsection{Incompressible flow} |
2356 | +For an incompressible multiphase flow, the continuity equation and momentum equation for phase $i$ (based on the derivation in \cite{ishii1975}, written in non-conservative form) are: |
2357 | \begin{equation} |
2358 | \sum_{i=1}^N\nabla\cdot\left(\alpha_i\mathbf{u}_i\right) = 0, |
2359 | \end{equation} |
2360 | @@ -769,6 +772,19 @@ |
2361 | \end{equation} |
2362 | where $\mathbf{u}_i$, $\rho_i$, $\mu_i$ and $\alpha_i$ are the velocity, density, isotropic viscosity and volume fraction of phase $i$ respectively, and $\mathbf{F}_i$ represents the forces imposed on phase $i$ by the other $N-1$ phases. Details of momentum transfer terms are given in Chapter \ref{chap:configuration}. |
2363 | |
2364 | +\subsubsection{Compressible flow} |
2365 | +Fluidity currently only supports the case where the continuous phase is compressible. All other phases (the dispersed phases) are incompressible. The momentum equation is the same as the one given above, but the continuity equation becomes: |
2366 | +\begin{equation} |
2367 | +\alpha_1\frac{\partial\rho_1}{\partial t} + \nabla\cdot\left(\alpha_1\rho_1\mathbf{u}_1\right) + \rho_1\left(\sum_{i=2}^N\nabla\cdot\left(\alpha_i\mathbf{u}_i\right)\right) = 0, |
2368 | +\end{equation} |
2369 | +where the compressible phase has an index of $i=1$ here. |
2370 | + |
2371 | +An internal energy equation may also need to be solved depending on the equation of state used for the compressible phase. This is given by: |
2372 | +\begin{equation} |
2373 | +\alpha_i\rho_i\frac{\partial e_i}{\partial t} + \alpha_i\rho_i\mathbf{u}_i\cdot\nabla e_i = -p\nabla\cdot\left(\alpha_i\mathbf{u}_i\right) + Q_i |
2374 | +\end{equation} |
2375 | +where $e_i$ is the internal energy and $Q_i$ is an inter-phase energy transfer term described in Chapter \ref{chap:configuration}. |
2376 | + |
2377 | \subsection{Porous Media Darcy Flow} |
2378 | \label{sec:porous_media_darcy_flow_equations} |
2379 | \index{porous media darcy flow} |
2380 | |
2381 | === modified file 'manual/visualisation_and_diagnostics.tex' |
2382 | --- manual/visualisation_and_diagnostics.tex 2012-03-06 10:47:48 +0000 |
2383 | +++ manual/visualisation_and_diagnostics.tex 2012-08-09 00:46:20 +0000 |
2384 | @@ -66,6 +66,7 @@ |
2385 | Limitations: Requires the diagnostic \option{scalar\_field::IsopycnalCoordinate} and (therefore) is not parallelised. |
2386 | \item[BulkMaterialPressure:]Calculates the bulk material pressure based on the MaterialDensity and MaterialVolumeFraction (and MaterialInternalEnergy if appropriate) for the equation of state of all materials. |
2387 | \item[CFLNumber:]CFLNumber as defined on the co-ordinate mesh. It is calculated as $\triangle t \vec{u} \mat{J}^{-1}$ where $\triangle t$ is the timestep, $\vec{u}$ the velocity and $\mat{J}$ the Jacobian. |
2388 | +\item[CompressibleContinuity:]Computes the residual of the compressible multiphase continuity equation. Used only in compressible multiphase flow simulations. |
2389 | \item[ControlVolumeCFLNumber:]CFL Number as defined on a control volume mesh. It is calculated as $\triangle t \frac{1}{V} \int _{cv_{face}} \vec{u}$ where $\triangle t$ is the timestep and $\vec{u}$ the velocity. The integral is taken over the faces of the control volume and $V$ is the volume of the control volume. |
2390 | \item[ControlVolumeDivergence:]Divergence of the velocity field where the divergence operator is defined using the control volume $\mathrm{C}^\mathrm{T}$ matrix. This assumes that the test space is discontinuous control volumes. |
2391 | \item[CVMaterialDensityCFLNumber:]Courant Number as defined on a control volume mesh and incorporating the MaterialDensity. Requires a MaterialDensity field! |
2392 | |
2393 | === modified file 'schemas/fluidity_options.rnc' |
2394 | --- schemas/fluidity_options.rnc 2012-06-27 14:47:10 +0000 |
2395 | +++ schemas/fluidity_options.rnc 2012-08-09 00:46:20 +0000 |
2396 | @@ -15,6 +15,7 @@ |
2397 | include "embedded_models.rnc" |
2398 | include "flredecomp.rnc" |
2399 | include "porous_media.rnc" |
2400 | +include "multiphase_interaction.rnc" |
2401 | include "equation_of_state.rnc" |
2402 | |
2403 | start = |
2404 | @@ -1387,6 +1388,17 @@ |
2405 | comment, |
2406 | element particle_diameter { |
2407 | real |
2408 | + }?, |
2409 | + ## If this is the fluid phase, |
2410 | + ## the effective conductivity is required |
2411 | + ## for the inter-phase heat transfer term. |
2412 | + element effective_conductivity { |
2413 | + real |
2414 | + }?, |
2415 | + ## The specific heat (at constant volume) is required |
2416 | + ## for the inter-phase heat transfer term. |
2417 | + element specific_heat { |
2418 | + real |
2419 | }? |
2420 | }?, |
2421 | sediment? |
2422 | @@ -1570,7 +1582,8 @@ |
2423 | reduced_model?, |
2424 | porous_media_model?, |
2425 | embedded_models?, |
2426 | - flredecomp? |
2427 | + flredecomp?, |
2428 | + multiphase_interaction? |
2429 | } |
2430 | ) |
2431 | |
2432 | @@ -4681,7 +4694,27 @@ |
2433 | } |
2434 | } |
2435 | ) |
2436 | - } |
2437 | + }| |
2438 | + |
2439 | + ## CompressibleContinuityResidual |
2440 | + ## |
2441 | + ## Computes the residual of the continuity equation used in compressible multiphase flow simulations |
2442 | + ## i.e. vfrac_c*d(rho_c)/dt + div(rho_c*vfrac_c*u_c) + \sum_i{ rho_c*div(vfrac_i*u_i) } |
2443 | + ## where _c and _i denote the compressible and incompressible phases respectively. |
2444 | + element scalar_field { |
2445 | + attribute rank { "0" }, |
2446 | + attribute name { "CompressibleContinuityResidual" }, |
2447 | + ( |
2448 | + element diagnostic { |
2449 | + velocity_mesh_choice, |
2450 | + internal_algorithm, |
2451 | + diagnostic_scalar_field_no_adapt, |
2452 | + element solver { |
2453 | + linear_solver_options_sym |
2454 | + } |
2455 | + } |
2456 | + ) |
2457 | + } |
2458 | |
2459 | # Insert new diagnostic scalar fields here using the template: |
2460 | # element scalar_field { |
2461 | |
2462 | === modified file 'schemas/fluidity_options.rng' |
2463 | --- schemas/fluidity_options.rng 2012-08-02 13:03:52 +0000 |
2464 | +++ schemas/fluidity_options.rng 2012-08-09 00:46:20 +0000 |
2465 | @@ -16,6 +16,7 @@ |
2466 | <include href="embedded_models.rng"/> |
2467 | <include href="flredecomp.rng"/> |
2468 | <include href="porous_media.rng"/> |
2469 | + <include href="multiphase_interaction.rng"/> |
2470 | <include href="equation_of_state.rng"/> |
2471 | <start> |
2472 | <element name="fluidity_options"> |
2473 | @@ -1852,6 +1853,21 @@ |
2474 | <ref name="real"/> |
2475 | </element> |
2476 | </optional> |
2477 | + <optional> |
2478 | + <element name="effective_conductivity"> |
2479 | + <a:documentation>If this is the fluid phase, |
2480 | +the effective conductivity is required |
2481 | +for the inter-phase heat transfer term.</a:documentation> |
2482 | + <ref name="real"/> |
2483 | + </element> |
2484 | + </optional> |
2485 | + <optional> |
2486 | + <element name="specific_heat"> |
2487 | + <a:documentation>The specific heat (at constant volume) is required |
2488 | +for the inter-phase heat transfer term.</a:documentation> |
2489 | + <ref name="real"/> |
2490 | + </element> |
2491 | + </optional> |
2492 | </element> |
2493 | </optional> |
2494 | <optional> |
2495 | @@ -2185,6 +2201,9 @@ |
2496 | <optional> |
2497 | <ref name="flredecomp"/> |
2498 | </optional> |
2499 | + <optional> |
2500 | + <ref name="multiphase_interaction"/> |
2501 | + </optional> |
2502 | </element> |
2503 | </start> |
2504 | <define name="sediment"> |
2505 | @@ -5935,6 +5954,27 @@ |
2506 | </element> |
2507 | </element> |
2508 | </element> |
2509 | + <element name="scalar_field"> |
2510 | + <a:documentation>CompressibleContinuityResidual |
2511 | + |
2512 | +Computes the residual of the continuity equation used in compressible multiphase flow simulations |
2513 | +i.e. vfrac_c*d(rho_c)/dt + div(rho_c*vfrac_c*u_c) + \sum_i{ rho_c*div(vfrac_i*u_i) } |
2514 | +where _c and _i denote the compressible and incompressible phases respectively.</a:documentation> |
2515 | + <attribute name="rank"> |
2516 | + <value>0</value> |
2517 | + </attribute> |
2518 | + <attribute name="name"> |
2519 | + <value>CompressibleContinuityResidual</value> |
2520 | + </attribute> |
2521 | + <element name="diagnostic"> |
2522 | + <ref name="velocity_mesh_choice"/> |
2523 | + <ref name="internal_algorithm"/> |
2524 | + <ref name="diagnostic_scalar_field_no_adapt"/> |
2525 | + <element name="solver"> |
2526 | + <ref name="linear_solver_options_sym"/> |
2527 | + </element> |
2528 | + </element> |
2529 | + </element> |
2530 | </choice> |
2531 | <!-- |
2532 | Insert new diagnostic scalar fields here using the template: |
2533 | |
2534 | === added file 'schemas/multiphase_interaction.rnc' |
2535 | --- schemas/multiphase_interaction.rnc 1970-01-01 00:00:00 +0000 |
2536 | +++ schemas/multiphase_interaction.rnc 2012-08-09 00:46:20 +0000 |
2537 | @@ -0,0 +1,58 @@ |
2538 | +# Phase interaction options for Fluidity's multiphase flow model. |
2539 | + |
2540 | +multiphase_interaction = |
2541 | + ( |
2542 | + ## Phase interaction options for Fluidity's multiphase flow model. |
2543 | + element multiphase_interaction { |
2544 | + comment, |
2545 | + fluid_particle_drag?, |
2546 | + heat_transfer? |
2547 | + } |
2548 | + ) |
2549 | + |
2550 | +fluid_particle_drag = |
2551 | + ( |
2552 | + ## Fluid-particle drag term. |
2553 | + element fluid_particle_drag { |
2554 | + drag_correlation |
2555 | + } |
2556 | + ) |
2557 | + |
2558 | +drag_correlation = |
2559 | + ( |
2560 | + ## Stokes drag correlation: 24/Re_p |
2561 | + ## where Re_p is the particle Reynolds number. |
2562 | + element drag_correlation { |
2563 | + attribute name { "stokes" }, |
2564 | + comment |
2565 | + }| |
2566 | + ## Fluid-particle drag term by Wen & Yu (1966). |
2567 | + element drag_correlation { |
2568 | + attribute name { "wen_yu" }, |
2569 | + comment |
2570 | + }| |
2571 | + ## Fluid-particle drag term by Ergun (1952). |
2572 | + element drag_correlation { |
2573 | + attribute name { "ergun" }, |
2574 | + comment |
2575 | + } |
2576 | + ) |
2577 | + |
2578 | +heat_transfer = |
2579 | + ( |
2580 | + ## Heat transfer term for the |
2581 | + ## multiphase internal energy equation. |
2582 | + ## Note: Only for fluid-particle phase pairs. |
2583 | + element heat_transfer { |
2584 | + heat_transfer_coefficient |
2585 | + } |
2586 | + ) |
2587 | + |
2588 | +heat_transfer_coefficient = |
2589 | + ( |
2590 | + ## Heat transfer coefficient by Gunn (1978). |
2591 | + element heat_transfer_coefficient { |
2592 | + attribute name { "gunn" }, |
2593 | + comment |
2594 | + } |
2595 | + ) |
2596 | \ No newline at end of file |
2597 | |
2598 | === added file 'schemas/multiphase_interaction.rng' |
2599 | --- schemas/multiphase_interaction.rng 1970-01-01 00:00:00 +0000 |
2600 | +++ schemas/multiphase_interaction.rng 2012-08-09 00:46:20 +0000 |
2601 | @@ -0,0 +1,65 @@ |
2602 | +<?xml version="1.0" encoding="UTF-8"?> |
2603 | +<grammar xmlns:a="http://relaxng.org/ns/compatibility/annotations/1.0" xmlns="http://relaxng.org/ns/structure/1.0"> |
2604 | + <!-- Phase interaction options for Fluidity's multiphase flow model. --> |
2605 | + <define name="multiphase_interaction"> |
2606 | + <element name="multiphase_interaction"> |
2607 | + <a:documentation>Phase interaction options for Fluidity's multiphase flow model.</a:documentation> |
2608 | + <ref name="comment"/> |
2609 | + <optional> |
2610 | + <ref name="fluid_particle_drag"/> |
2611 | + </optional> |
2612 | + <optional> |
2613 | + <ref name="heat_transfer"/> |
2614 | + </optional> |
2615 | + </element> |
2616 | + </define> |
2617 | + <define name="fluid_particle_drag"> |
2618 | + <element name="fluid_particle_drag"> |
2619 | + <a:documentation>Fluid-particle drag term.</a:documentation> |
2620 | + <ref name="drag_correlation"/> |
2621 | + </element> |
2622 | + </define> |
2623 | + <define name="drag_correlation"> |
2624 | + <choice> |
2625 | + <element name="drag_correlation"> |
2626 | + <a:documentation>Stokes drag correlation: 24/Re_p |
2627 | +where Re_p is the particle Reynolds number.</a:documentation> |
2628 | + <attribute name="name"> |
2629 | + <value>stokes</value> |
2630 | + </attribute> |
2631 | + <ref name="comment"/> |
2632 | + </element> |
2633 | + <element name="drag_correlation"> |
2634 | + <a:documentation>Fluid-particle drag term by Wen & Yu (1966).</a:documentation> |
2635 | + <attribute name="name"> |
2636 | + <value>wen_yu</value> |
2637 | + </attribute> |
2638 | + <ref name="comment"/> |
2639 | + </element> |
2640 | + <element name="drag_correlation"> |
2641 | + <a:documentation>Fluid-particle drag term by Ergun (1952).</a:documentation> |
2642 | + <attribute name="name"> |
2643 | + <value>ergun</value> |
2644 | + </attribute> |
2645 | + <ref name="comment"/> |
2646 | + </element> |
2647 | + </choice> |
2648 | + </define> |
2649 | + <define name="heat_transfer"> |
2650 | + <element name="heat_transfer"> |
2651 | + <a:documentation>Heat transfer term for the |
2652 | +multiphase internal energy equation. |
2653 | +Note: Only for fluid-particle phase pairs.</a:documentation> |
2654 | + <ref name="heat_transfer_coefficient"/> |
2655 | + </element> |
2656 | + </define> |
2657 | + <define name="heat_transfer_coefficient"> |
2658 | + <element name="heat_transfer_coefficient"> |
2659 | + <a:documentation>Heat transfer coefficient by Gunn (1978).</a:documentation> |
2660 | + <attribute name="name"> |
2661 | + <value>gunn</value> |
2662 | + </attribute> |
2663 | + <ref name="comment"/> |
2664 | + </element> |
2665 | + </define> |
2666 | +</grammar> |
2667 | |
2668 | === added directory 'tests/flux_bc' |
2669 | === added file 'tests/flux_bc/Makefile' |
2670 | --- tests/flux_bc/Makefile 1970-01-01 00:00:00 +0000 |
2671 | +++ tests/flux_bc/Makefile 2012-08-09 00:46:20 +0000 |
2672 | @@ -0,0 +1,22 @@ |
2673 | +preprocess: |
2674 | + @echo **********Creating meshes |
2675 | + gmsh -2 -o flux_bc_2d.msh src/flux_bc_2d.geo |
2676 | + gmsh -3 -o flux_bc_3d.msh src/flux_bc_3d.geo |
2677 | + ../../bin/gmsh2triangle --2d flux_bc_2d.msh |
2678 | + ../../bin/gmsh2triangle flux_bc_3d.msh |
2679 | + |
2680 | +run: |
2681 | + @echo **********Running 2D simulation |
2682 | + ../../bin/fluidity flux_bc_2d.flml |
2683 | + @echo **********Running 3D simulation |
2684 | + ../../bin/fluidity flux_bc_3d.flml |
2685 | + |
2686 | +input: clean preprocess |
2687 | + |
2688 | +clean: |
2689 | + rm -f *.stat *.steady_state* |
2690 | + rm -f *.d.* *.vtu |
2691 | + rm -f *.msh |
2692 | + rm -f *.ele *.edge *.node *.poly *.face |
2693 | + rm -f matrixdump* *checkpoint* *.log* *.err* |
2694 | + |
2695 | |
2696 | === added file 'tests/flux_bc/flux_bc.xml' |
2697 | --- tests/flux_bc/flux_bc.xml 1970-01-01 00:00:00 +0000 |
2698 | +++ tests/flux_bc/flux_bc.xml 2012-08-09 00:46:20 +0000 |
2699 | @@ -0,0 +1,44 @@ |
2700 | +<?xml version="1.0" encoding="UTF-8" ?> |
2701 | +<!DOCTYPE testproblem SYSTEM "regressiontest.dtd"> |
2702 | + |
2703 | +<testproblem> |
2704 | + |
2705 | + <name>flux_bc</name> |
2706 | + <owner userid="ctj10"/> |
2707 | + <tags>flml</tags> |
2708 | + |
2709 | + <problem_definition length="medium" nprocs="1"> |
2710 | + <command_line>make run</command_line> |
2711 | + <!-- This checks the integral of the PhaseVolumeFraction field at t = 1 second. --> |
2712 | + <!-- The volumetric flux (r) is set to 0.5 m^3/s/m^2, and the flux boundary condition --> |
2713 | + <!-- enforces d/dt \int_{volume} vfrac = \int_{surface} r --> |
2714 | + </problem_definition> |
2715 | + |
2716 | + <variables> |
2717 | + <variable name="vfrac_integral_2d" language="python"> |
2718 | +from fluidity_tools import stat_parser |
2719 | +s = stat_parser("flux_bc_2d.stat") |
2720 | +vfrac_integral_2d = s["Tephra"]["PhaseVolumeFraction"]["integral"] |
2721 | + </variable> |
2722 | + |
2723 | + <variable name="vfrac_integral_3d" language="python"> |
2724 | +from fluidity_tools import stat_parser |
2725 | +s = stat_parser("flux_bc_3d.stat") |
2726 | +vfrac_integral_3d = s["Tephra"]["PhaseVolumeFraction"]["integral"] |
2727 | + </variable> |
2728 | + </variables> |
2729 | + |
2730 | + <pass_tests> |
2731 | + <test name="Integral of PhaseVolumeFraction after 1 second = 0.15 in 2D simulation." language="python"> |
2732 | +assert abs(vfrac_integral_2d[10] - 0.15) < 1.0e-9 |
2733 | + </test> |
2734 | + |
2735 | + <test name="Integral of PhaseVolumeFraction after 1 second = 0.045 in 3D simulation." language="python"> |
2736 | +assert abs(vfrac_integral_3d[10] - 0.045) < 1.0e-9 |
2737 | + </test> |
2738 | + </pass_tests> |
2739 | + |
2740 | + <warn_tests> |
2741 | + </warn_tests> |
2742 | + |
2743 | +</testproblem> |
2744 | |
2745 | === added file 'tests/flux_bc/flux_bc_2d.flml' |
2746 | --- tests/flux_bc/flux_bc_2d.flml 1970-01-01 00:00:00 +0000 |
2747 | +++ tests/flux_bc/flux_bc_2d.flml 2012-08-09 00:46:20 +0000 |
2748 | @@ -0,0 +1,396 @@ |
2749 | +<?xml version='1.0' encoding='utf-8'?> |
2750 | +<fluidity_options> |
2751 | + <simulation_name> |
2752 | + <string_value lines="1">flux_bc_2d</string_value> |
2753 | + </simulation_name> |
2754 | + <problem_type> |
2755 | + <string_value lines="1">fluids</string_value> |
2756 | + </problem_type> |
2757 | + <geometry> |
2758 | + <dimension> |
2759 | + <integer_value rank="0">2</integer_value> |
2760 | + </dimension> |
2761 | + <mesh name="CoordinateMesh"> |
2762 | + <from_file file_name="flux_bc_2d"> |
2763 | + <format name="triangle"/> |
2764 | + <stat> |
2765 | + <include_in_stat/> |
2766 | + </stat> |
2767 | + </from_file> |
2768 | + </mesh> |
2769 | + <mesh name="VelocityMesh"> |
2770 | + <from_mesh> |
2771 | + <mesh name="CoordinateMesh"/> |
2772 | + <mesh_shape> |
2773 | + <polynomial_degree> |
2774 | + <integer_value rank="0">1</integer_value> |
2775 | + </polynomial_degree> |
2776 | + </mesh_shape> |
2777 | + <mesh_continuity> |
2778 | + <string_value>discontinuous</string_value> |
2779 | + </mesh_continuity> |
2780 | + <stat> |
2781 | + <exclude_from_stat/> |
2782 | + </stat> |
2783 | + </from_mesh> |
2784 | + </mesh> |
2785 | + <mesh name="PressureMesh"> |
2786 | + <from_mesh> |
2787 | + <mesh name="CoordinateMesh"/> |
2788 | + <mesh_shape> |
2789 | + <polynomial_degree> |
2790 | + <integer_value rank="0">2</integer_value> |
2791 | + </polynomial_degree> |
2792 | + </mesh_shape> |
2793 | + <stat> |
2794 | + <exclude_from_stat/> |
2795 | + </stat> |
2796 | + </from_mesh> |
2797 | + </mesh> |
2798 | + <quadrature> |
2799 | + <degree> |
2800 | + <integer_value rank="0">4</integer_value> |
2801 | + </degree> |
2802 | + </quadrature> |
2803 | + </geometry> |
2804 | + <io> |
2805 | + <dump_format> |
2806 | + <string_value>vtk</string_value> |
2807 | + </dump_format> |
2808 | + <dump_period> |
2809 | + <constant> |
2810 | + <real_value rank="0">1.0</real_value> |
2811 | + </constant> |
2812 | + </dump_period> |
2813 | + <output_mesh name="VelocityMesh"/> |
2814 | + <stat> |
2815 | + <output_at_start/> |
2816 | + </stat> |
2817 | + </io> |
2818 | + <timestepping> |
2819 | + <current_time> |
2820 | + <real_value rank="0">0.0</real_value> |
2821 | + </current_time> |
2822 | + <timestep> |
2823 | + <real_value rank="0">0.1</real_value> |
2824 | + </timestep> |
2825 | + <finish_time> |
2826 | + <real_value rank="0">1.0</real_value> |
2827 | + </finish_time> |
2828 | + <nonlinear_iterations> |
2829 | + <integer_value rank="0">2</integer_value> |
2830 | + <tolerance> |
2831 | + <real_value rank="0">1.0e-12</real_value> |
2832 | + <infinity_norm/> |
2833 | + </tolerance> |
2834 | + </nonlinear_iterations> |
2835 | + </timestepping> |
2836 | + <physical_parameters> |
2837 | + <gravity> |
2838 | + <magnitude> |
2839 | + <real_value rank="0">9.8</real_value> |
2840 | + </magnitude> |
2841 | + <vector_field name="GravityDirection" rank="1"> |
2842 | + <prescribed> |
2843 | + <mesh name="CoordinateMesh"/> |
2844 | + <value name="WholeMesh"> |
2845 | + <constant> |
2846 | + <real_value shape="2" dim1="dim" rank="1">0.0 -1.0</real_value> |
2847 | + </constant> |
2848 | + </value> |
2849 | + <output/> |
2850 | + <stat> |
2851 | + <include_in_stat/> |
2852 | + </stat> |
2853 | + <detectors> |
2854 | + <exclude_from_detectors/> |
2855 | + </detectors> |
2856 | + </prescribed> |
2857 | + </vector_field> |
2858 | + </gravity> |
2859 | + </physical_parameters> |
2860 | + <material_phase name="Tephra"> |
2861 | + <equation_of_state> |
2862 | + <fluids> |
2863 | + <linear> |
2864 | + <reference_density> |
2865 | + <real_value rank="0">2340</real_value> |
2866 | + </reference_density> |
2867 | + </linear> |
2868 | + </fluids> |
2869 | + </equation_of_state> |
2870 | + <scalar_field name="Pressure" rank="0"> |
2871 | + <prognostic> |
2872 | + <mesh name="PressureMesh"/> |
2873 | + <spatial_discretisation> |
2874 | + <continuous_galerkin> |
2875 | + <remove_stabilisation_term/> |
2876 | + </continuous_galerkin> |
2877 | + </spatial_discretisation> |
2878 | + <reference_node> |
2879 | + <integer_value rank="0">1</integer_value> |
2880 | + </reference_node> |
2881 | + <scheme> |
2882 | + <poisson_pressure_solution> |
2883 | + <string_value lines="1">never</string_value> |
2884 | + </poisson_pressure_solution> |
2885 | + <use_projection_method/> |
2886 | + </scheme> |
2887 | + <solver> |
2888 | + <iterative_method name="cg"/> |
2889 | + <preconditioner name="sor"/> |
2890 | + <relative_error> |
2891 | + <real_value rank="0">1.0e-7</real_value> |
2892 | + </relative_error> |
2893 | + <max_iterations> |
2894 | + <integer_value rank="0">1000</integer_value> |
2895 | + </max_iterations> |
2896 | + <never_ignore_solver_failures/> |
2897 | + <diagnostics> |
2898 | + <monitors/> |
2899 | + </diagnostics> |
2900 | + </solver> |
2901 | + <output/> |
2902 | + <stat/> |
2903 | + <convergence> |
2904 | + <include_in_convergence/> |
2905 | + </convergence> |
2906 | + <detectors> |
2907 | + <exclude_from_detectors/> |
2908 | + </detectors> |
2909 | + <steady_state> |
2910 | + <include_in_steady_state/> |
2911 | + </steady_state> |
2912 | + <no_interpolation/> |
2913 | + </prognostic> |
2914 | + </scalar_field> |
2915 | + <scalar_field name="Density" rank="0"> |
2916 | + <diagnostic> |
2917 | + <algorithm name="Internal" material_phase_support="multiple"/> |
2918 | + <mesh name="CoordinateMesh"/> |
2919 | + <output/> |
2920 | + <stat/> |
2921 | + <convergence> |
2922 | + <include_in_convergence/> |
2923 | + </convergence> |
2924 | + <detectors> |
2925 | + <include_in_detectors/> |
2926 | + </detectors> |
2927 | + <steady_state> |
2928 | + <include_in_steady_state/> |
2929 | + </steady_state> |
2930 | + </diagnostic> |
2931 | + </scalar_field> |
2932 | + <vector_field name="Velocity" rank="1"> |
2933 | + <prognostic> |
2934 | + <mesh name="VelocityMesh"/> |
2935 | + <equation name="LinearMomentum"/> |
2936 | + <spatial_discretisation> |
2937 | + <discontinuous_galerkin> |
2938 | + <viscosity_scheme> |
2939 | + <bassi_rebay/> |
2940 | + </viscosity_scheme> |
2941 | + <advection_scheme> |
2942 | + <upwind/> |
2943 | + <integrate_advection_by_parts> |
2944 | + <twice/> |
2945 | + </integrate_advection_by_parts> |
2946 | + </advection_scheme> |
2947 | + </discontinuous_galerkin> |
2948 | + <conservative_advection> |
2949 | + <real_value rank="0">0.0</real_value> |
2950 | + </conservative_advection> |
2951 | + </spatial_discretisation> |
2952 | + <temporal_discretisation> |
2953 | + <theta> |
2954 | + <real_value rank="0">1.0</real_value> |
2955 | + </theta> |
2956 | + <relaxation> |
2957 | + <real_value rank="0">0.5</real_value> |
2958 | + </relaxation> |
2959 | + </temporal_discretisation> |
2960 | + <solver> |
2961 | + <iterative_method name="gmres"> |
2962 | + <restart> |
2963 | + <integer_value rank="0">30</integer_value> |
2964 | + </restart> |
2965 | + </iterative_method> |
2966 | + <preconditioner name="sor"/> |
2967 | + <relative_error> |
2968 | + <real_value rank="0">1.0e-7</real_value> |
2969 | + </relative_error> |
2970 | + <max_iterations> |
2971 | + <integer_value rank="0">1000</integer_value> |
2972 | + </max_iterations> |
2973 | + <never_ignore_solver_failures/> |
2974 | + <diagnostics> |
2975 | + <monitors/> |
2976 | + </diagnostics> |
2977 | + </solver> |
2978 | + <initial_condition name="WholeMesh"> |
2979 | + <constant> |
2980 | + <real_value shape="2" dim1="dim" rank="1">0.0 0.0</real_value> |
2981 | + </constant> |
2982 | + </initial_condition> |
2983 | + <boundary_conditions name="Sides"> |
2984 | + <surface_ids> |
2985 | + <integer_value shape="1" rank="1">666</integer_value> |
2986 | + </surface_ids> |
2987 | + <type name="dirichlet"> |
2988 | + <apply_weakly/> |
2989 | + <align_bc_with_cartesian> |
2990 | + <x_component> |
2991 | + <constant> |
2992 | + <real_value rank="0">0.0</real_value> |
2993 | + </constant> |
2994 | + </x_component> |
2995 | + </align_bc_with_cartesian> |
2996 | + </type> |
2997 | + </boundary_conditions> |
2998 | + <boundary_conditions name="Bottom"> |
2999 | + <surface_ids> |
3000 | + <integer_value shape="1" rank="1">444</integer_value> |
3001 | + </surface_ids> |
3002 | + <type name="dirichlet"> |
3003 | + <apply_weakly/> |
3004 | + <align_bc_with_cartesian> |
3005 | + <y_component> |
3006 | + <constant> |
3007 | + <real_value rank="0">0.0</real_value> |
3008 | + </constant> |
3009 | + </y_component> |
3010 | + </align_bc_with_cartesian> |
3011 | + </type> |
3012 | + </boundary_conditions> |
3013 | + <boundary_conditions name="Top"> |
3014 | + <surface_ids> |
3015 | + <integer_value shape="1" rank="1">333</integer_value> |
3016 | + </surface_ids> |
3017 | + <type name="dirichlet"> |
3018 | + <apply_weakly/> |
3019 | + <align_bc_with_cartesian> |
3020 | + <y_component> |
3021 | + <constant> |
3022 | + <real_value rank="0">0.0</real_value> |
3023 | + </constant> |
3024 | + </y_component> |
3025 | + </align_bc_with_cartesian> |
3026 | + </type> |
3027 | + </boundary_conditions> |
3028 | + <tensor_field name="Viscosity" rank="2"> |
3029 | + <prescribed> |
3030 | + <value name="WholeMesh"> |
3031 | + <isotropic> |
3032 | + <constant> |
3033 | + <real_value rank="0">0.001</real_value> |
3034 | + </constant> |
3035 | + </isotropic> |
3036 | + </value> |
3037 | + <output/> |
3038 | + </prescribed> |
3039 | + </tensor_field> |
3040 | + <output/> |
3041 | + <stat> |
3042 | + <include_in_stat/> |
3043 | + <previous_time_step> |
3044 | + <exclude_from_stat/> |
3045 | + </previous_time_step> |
3046 | + <nonlinear_field> |
3047 | + <exclude_from_stat/> |
3048 | + </nonlinear_field> |
3049 | + </stat> |
3050 | + <convergence> |
3051 | + <include_in_convergence/> |
3052 | + </convergence> |
3053 | + <detectors> |
3054 | + <include_in_detectors/> |
3055 | + </detectors> |
3056 | + <steady_state> |
3057 | + <include_in_steady_state/> |
3058 | + </steady_state> |
3059 | + <consistent_interpolation/> |
3060 | + </prognostic> |
3061 | + </vector_field> |
3062 | + <scalar_field name="PhaseVolumeFraction" rank="0"> |
3063 | + <prognostic> |
3064 | + <mesh name="CoordinateMesh"/> |
3065 | + <equation name="AdvectionDiffusion"/> |
3066 | + <spatial_discretisation> |
3067 | + <control_volumes> |
3068 | + <face_value name="FiniteElement"> |
3069 | + <limit_face_value> |
3070 | + <limiter name="Sweby"/> |
3071 | + </limit_face_value> |
3072 | + </face_value> |
3073 | + <diffusion_scheme name="BassiRebay"/> |
3074 | + </control_volumes> |
3075 | + <conservative_advection> |
3076 | + <real_value rank="0">1.0</real_value> |
3077 | + </conservative_advection> |
3078 | + </spatial_discretisation> |
3079 | + <temporal_discretisation> |
3080 | + <theta> |
3081 | + <real_value rank="0">1.0</real_value> |
3082 | + </theta> |
3083 | + </temporal_discretisation> |
3084 | + <solver> |
3085 | + <iterative_method name="gmres"> |
3086 | + <restart> |
3087 | + <integer_value rank="0">30</integer_value> |
3088 | + </restart> |
3089 | + </iterative_method> |
3090 | + <preconditioner name="sor"/> |
3091 | + <relative_error> |
3092 | + <real_value rank="0">1.0e-7</real_value> |
3093 | + </relative_error> |
3094 | + <max_iterations> |
3095 | + <integer_value rank="0">1000</integer_value> |
3096 | + </max_iterations> |
3097 | + <never_ignore_solver_failures/> |
3098 | + <diagnostics> |
3099 | + <monitors/> |
3100 | + </diagnostics> |
3101 | + </solver> |
3102 | + <initial_condition name="WholeMesh"> |
3103 | + <constant> |
3104 | + <real_value rank="0">0.0</real_value> |
3105 | + </constant> |
3106 | + </initial_condition> |
3107 | + <boundary_conditions name="Top"> |
3108 | + <surface_ids> |
3109 | + <integer_value shape="1" rank="1">333</integer_value> |
3110 | + </surface_ids> |
3111 | + <type name="flux"> |
3112 | + <constant> |
3113 | + <real_value rank="0">0.5</real_value> |
3114 | + </constant> |
3115 | + </type> |
3116 | + </boundary_conditions> |
3117 | + <boundary_conditions name="Bottom"> |
3118 | + <surface_ids> |
3119 | + <integer_value shape="1" rank="1">444</integer_value> |
3120 | + </surface_ids> |
3121 | + <type name="zero_flux"/> |
3122 | + </boundary_conditions> |
3123 | + <boundary_conditions name="Sides"> |
3124 | + <surface_ids> |
3125 | + <integer_value shape="1" rank="1">666</integer_value> |
3126 | + </surface_ids> |
3127 | + <type name="zero_flux"/> |
3128 | + </boundary_conditions> |
3129 | + <output/> |
3130 | + <stat/> |
3131 | + <convergence> |
3132 | + <include_in_convergence/> |
3133 | + </convergence> |
3134 | + <detectors> |
3135 | + <include_in_detectors/> |
3136 | + </detectors> |
3137 | + <steady_state> |
3138 | + <include_in_steady_state/> |
3139 | + </steady_state> |
3140 | + <consistent_interpolation/> |
3141 | + </prognostic> |
3142 | + </scalar_field> |
3143 | + </material_phase> |
3144 | +</fluidity_options> |
3145 | |
3146 | === added file 'tests/flux_bc/flux_bc_3d.flml' |
3147 | --- tests/flux_bc/flux_bc_3d.flml 1970-01-01 00:00:00 +0000 |
3148 | +++ tests/flux_bc/flux_bc_3d.flml 2012-08-09 00:46:20 +0000 |
3149 | @@ -0,0 +1,405 @@ |
3150 | +<?xml version='1.0' encoding='utf-8'?> |
3151 | +<fluidity_options> |
3152 | + <simulation_name> |
3153 | + <string_value lines="1">flux_bc_3d</string_value> |
3154 | + </simulation_name> |
3155 | + <problem_type> |
3156 | + <string_value lines="1">fluids</string_value> |
3157 | + </problem_type> |
3158 | + <geometry> |
3159 | + <dimension> |
3160 | + <integer_value rank="0">3</integer_value> |
3161 | + </dimension> |
3162 | + <mesh name="CoordinateMesh"> |
3163 | + <from_file file_name="flux_bc_3d"> |
3164 | + <format name="triangle"/> |
3165 | + <stat> |
3166 | + <include_in_stat/> |
3167 | + </stat> |
3168 | + </from_file> |
3169 | + </mesh> |
3170 | + <mesh name="VelocityMesh"> |
3171 | + <from_mesh> |
3172 | + <mesh name="CoordinateMesh"/> |
3173 | + <mesh_shape> |
3174 | + <polynomial_degree> |
3175 | + <integer_value rank="0">1</integer_value> |
3176 | + </polynomial_degree> |
3177 | + </mesh_shape> |
3178 | + <mesh_continuity> |
3179 | + <string_value>discontinuous</string_value> |
3180 | + </mesh_continuity> |
3181 | + <stat> |
3182 | + <exclude_from_stat/> |
3183 | + </stat> |
3184 | + </from_mesh> |
3185 | + </mesh> |
3186 | + <mesh name="PressureMesh"> |
3187 | + <from_mesh> |
3188 | + <mesh name="CoordinateMesh"/> |
3189 | + <mesh_shape> |
3190 | + <polynomial_degree> |
3191 | + <integer_value rank="0">2</integer_value> |
3192 | + </polynomial_degree> |
3193 | + </mesh_shape> |
3194 | + <stat> |
3195 | + <exclude_from_stat/> |
3196 | + </stat> |
3197 | + </from_mesh> |
3198 | + </mesh> |
3199 | + <quadrature> |
3200 | + <degree> |
3201 | + <integer_value rank="0">4</integer_value> |
3202 | + </degree> |
3203 | + </quadrature> |
3204 | + </geometry> |
3205 | + <io> |
3206 | + <dump_format> |
3207 | + <string_value>vtk</string_value> |
3208 | + </dump_format> |
3209 | + <dump_period> |
3210 | + <constant> |
3211 | + <real_value rank="0">1.0</real_value> |
3212 | + </constant> |
3213 | + </dump_period> |
3214 | + <output_mesh name="VelocityMesh"/> |
3215 | + <stat> |
3216 | + <output_at_start/> |
3217 | + </stat> |
3218 | + </io> |
3219 | + <timestepping> |
3220 | + <current_time> |
3221 | + <real_value rank="0">0.0</real_value> |
3222 | + </current_time> |
3223 | + <timestep> |
3224 | + <real_value rank="0">0.1</real_value> |
3225 | + </timestep> |
3226 | + <finish_time> |
3227 | + <real_value rank="0">1.0</real_value> |
3228 | + </finish_time> |
3229 | + <nonlinear_iterations> |
3230 | + <integer_value rank="0">2</integer_value> |
3231 | + <tolerance> |
3232 | + <real_value rank="0">1.0e-12</real_value> |
3233 | + <infinity_norm/> |
3234 | + </tolerance> |
3235 | + </nonlinear_iterations> |
3236 | + </timestepping> |
3237 | + <physical_parameters> |
3238 | + <gravity> |
3239 | + <magnitude> |
3240 | + <real_value rank="0">9.8</real_value> |
3241 | + </magnitude> |
3242 | + <vector_field name="GravityDirection" rank="1"> |
3243 | + <prescribed> |
3244 | + <mesh name="CoordinateMesh"/> |
3245 | + <value name="WholeMesh"> |
3246 | + <constant> |
3247 | + <real_value shape="3" dim1="dim" rank="1">0.0 -1.0 0.0</real_value> |
3248 | + </constant> |
3249 | + </value> |
3250 | + <output/> |
3251 | + <stat> |
3252 | + <include_in_stat/> |
3253 | + </stat> |
3254 | + <detectors> |
3255 | + <exclude_from_detectors/> |
3256 | + </detectors> |
3257 | + </prescribed> |
3258 | + </vector_field> |
3259 | + </gravity> |
3260 | + </physical_parameters> |
3261 | + <material_phase name="Tephra"> |
3262 | + <equation_of_state> |
3263 | + <fluids> |
3264 | + <linear> |
3265 | + <reference_density> |
3266 | + <real_value rank="0">2340.0</real_value> |
3267 | + </reference_density> |
3268 | + </linear> |
3269 | + </fluids> |
3270 | + </equation_of_state> |
3271 | + <scalar_field name="Pressure" rank="0"> |
3272 | + <prognostic> |
3273 | + <mesh name="PressureMesh"/> |
3274 | + <spatial_discretisation> |
3275 | + <continuous_galerkin> |
3276 | + <remove_stabilisation_term/> |
3277 | + </continuous_galerkin> |
3278 | + </spatial_discretisation> |
3279 | + <reference_node> |
3280 | + <integer_value rank="0">1</integer_value> |
3281 | + </reference_node> |
3282 | + <scheme> |
3283 | + <poisson_pressure_solution> |
3284 | + <string_value lines="1">never</string_value> |
3285 | + </poisson_pressure_solution> |
3286 | + <use_projection_method/> |
3287 | + </scheme> |
3288 | + <solver> |
3289 | + <iterative_method name="cg"/> |
3290 | + <preconditioner name="sor"/> |
3291 | + <relative_error> |
3292 | + <real_value rank="0">1.0e-7</real_value> |
3293 | + </relative_error> |
3294 | + <max_iterations> |
3295 | + <integer_value rank="0">1000</integer_value> |
3296 | + </max_iterations> |
3297 | + <never_ignore_solver_failures/> |
3298 | + <diagnostics> |
3299 | + <monitors/> |
3300 | + </diagnostics> |
3301 | + </solver> |
3302 | + <output/> |
3303 | + <stat/> |
3304 | + <convergence> |
3305 | + <include_in_convergence/> |
3306 | + </convergence> |
3307 | + <detectors> |
3308 | + <exclude_from_detectors/> |
3309 | + </detectors> |
3310 | + <steady_state> |
3311 | + <include_in_steady_state/> |
3312 | + </steady_state> |
3313 | + <no_interpolation/> |
3314 | + </prognostic> |
3315 | + </scalar_field> |
3316 | + <scalar_field name="Density" rank="0"> |
3317 | + <diagnostic> |
3318 | + <algorithm name="Internal" material_phase_support="multiple"/> |
3319 | + <mesh name="CoordinateMesh"/> |
3320 | + <output/> |
3321 | + <stat/> |
3322 | + <convergence> |
3323 | + <include_in_convergence/> |
3324 | + </convergence> |
3325 | + <detectors> |
3326 | + <include_in_detectors/> |
3327 | + </detectors> |
3328 | + <steady_state> |
3329 | + <include_in_steady_state/> |
3330 | + </steady_state> |
3331 | + </diagnostic> |
3332 | + </scalar_field> |
3333 | + <vector_field name="Velocity" rank="1"> |
3334 | + <prognostic> |
3335 | + <mesh name="VelocityMesh"/> |
3336 | + <equation name="LinearMomentum"/> |
3337 | + <spatial_discretisation> |
3338 | + <discontinuous_galerkin> |
3339 | + <viscosity_scheme> |
3340 | + <bassi_rebay/> |
3341 | + </viscosity_scheme> |
3342 | + <advection_scheme> |
3343 | + <upwind/> |
3344 | + <integrate_advection_by_parts> |
3345 | + <twice/> |
3346 | + </integrate_advection_by_parts> |
3347 | + </advection_scheme> |
3348 | + </discontinuous_galerkin> |
3349 | + <conservative_advection> |
3350 | + <real_value rank="0">0.0</real_value> |
3351 | + </conservative_advection> |
3352 | + </spatial_discretisation> |
3353 | + <temporal_discretisation> |
3354 | + <theta> |
3355 | + <real_value rank="0">1.0</real_value> |
3356 | + </theta> |
3357 | + <relaxation> |
3358 | + <real_value rank="0">0.5</real_value> |
3359 | + </relaxation> |
3360 | + </temporal_discretisation> |
3361 | + <solver> |
3362 | + <iterative_method name="gmres"> |
3363 | + <restart> |
3364 | + <integer_value rank="0">30</integer_value> |
3365 | + </restart> |
3366 | + </iterative_method> |
3367 | + <preconditioner name="sor"/> |
3368 | + <relative_error> |
3369 | + <real_value rank="0">1.0e-7</real_value> |
3370 | + </relative_error> |
3371 | + <max_iterations> |
3372 | + <integer_value rank="0">1000</integer_value> |
3373 | + </max_iterations> |
3374 | + <never_ignore_solver_failures/> |
3375 | + <diagnostics> |
3376 | + <monitors/> |
3377 | + </diagnostics> |
3378 | + </solver> |
3379 | + <initial_condition name="WholeMesh"> |
3380 | + <constant> |
3381 | + <real_value shape="3" dim1="dim" rank="1">0.0 0.0 0.0</real_value> |
3382 | + </constant> |
3383 | + </initial_condition> |
3384 | + <boundary_conditions name="Top"> |
3385 | + <surface_ids> |
3386 | + <integer_value shape="1" rank="1">444</integer_value> |
3387 | + </surface_ids> |
3388 | + <type name="dirichlet"> |
3389 | + <apply_weakly/> |
3390 | + <align_bc_with_cartesian> |
3391 | + <y_component> |
3392 | + <constant> |
3393 | + <real_value rank="0">0</real_value> |
3394 | + </constant> |
3395 | + </y_component> |
3396 | + </align_bc_with_cartesian> |
3397 | + </type> |
3398 | + </boundary_conditions> |
3399 | + <boundary_conditions name="Sides_X"> |
3400 | + <surface_ids> |
3401 | + <integer_value shape="1" rank="1">777</integer_value> |
3402 | + </surface_ids> |
3403 | + <type name="dirichlet"> |
3404 | + <apply_weakly/> |
3405 | + <align_bc_with_cartesian> |
3406 | + <x_component> |
3407 | + <constant> |
3408 | + <real_value rank="0">0.0</real_value> |
3409 | + </constant> |
3410 | + </x_component> |
3411 | + </align_bc_with_cartesian> |
3412 | + </type> |
3413 | + </boundary_conditions> |
3414 | + <boundary_conditions name="Sides_Z"> |
3415 | + <surface_ids> |
3416 | + <integer_value shape="1" rank="1">666</integer_value> |
3417 | + </surface_ids> |
3418 | + <type name="dirichlet"> |
3419 | + <apply_weakly/> |
3420 | + <align_bc_with_cartesian> |
3421 | + <z_component> |
3422 | + <constant> |
3423 | + <real_value rank="0">0.0</real_value> |
3424 | + </constant> |
3425 | + </z_component> |
3426 | + </align_bc_with_cartesian> |
3427 | + </type> |
3428 | + </boundary_conditions> |
3429 | + <boundary_conditions name="Bottom"> |
3430 | + <surface_ids> |
3431 | + <integer_value shape="1" rank="1">333</integer_value> |
3432 | + </surface_ids> |
3433 | + <type name="dirichlet"> |
3434 | + <apply_weakly/> |
3435 | + <align_bc_with_cartesian> |
3436 | + <y_component> |
3437 | + <constant> |
3438 | + <real_value rank="0">0.0</real_value> |
3439 | + </constant> |
3440 | + </y_component> |
3441 | + </align_bc_with_cartesian> |
3442 | + </type> |
3443 | + </boundary_conditions> |
3444 | + <tensor_field name="Viscosity" rank="2"> |
3445 | + <prescribed> |
3446 | + <value name="WholeMesh"> |
3447 | + <isotropic> |
3448 | + <constant> |
3449 | + <real_value rank="0">0.001</real_value> |
3450 | + </constant> |
3451 | + </isotropic> |
3452 | + </value> |
3453 | + <output/> |
3454 | + </prescribed> |
3455 | + </tensor_field> |
3456 | + <output/> |
3457 | + <stat> |
3458 | + <include_in_stat/> |
3459 | + <previous_time_step> |
3460 | + <exclude_from_stat/> |
3461 | + </previous_time_step> |
3462 | + <nonlinear_field> |
3463 | + <exclude_from_stat/> |
3464 | + </nonlinear_field> |
3465 | + </stat> |
3466 | + <convergence> |
3467 | + <include_in_convergence/> |
3468 | + </convergence> |
3469 | + <detectors> |
3470 | + <include_in_detectors/> |
3471 | + </detectors> |
3472 | + <steady_state> |
3473 | + <include_in_steady_state/> |
3474 | + </steady_state> |
3475 | + <consistent_interpolation/> |
3476 | + </prognostic> |
3477 | + </vector_field> |
3478 | + <scalar_field name="PhaseVolumeFraction" rank="0"> |
3479 | + <prognostic> |
3480 | + <mesh name="CoordinateMesh"/> |
3481 | + <equation name="AdvectionDiffusion"/> |
3482 | + <spatial_discretisation> |
3483 | + <control_volumes> |
3484 | + <face_value name="FiniteElement"> |
3485 | + <limit_face_value> |
3486 | + <limiter name="Sweby"/> |
3487 | + </limit_face_value> |
3488 | + </face_value> |
3489 | + <diffusion_scheme name="BassiRebay"/> |
3490 | + </control_volumes> |
3491 | + <conservative_advection> |
3492 | + <real_value rank="0">1.0</real_value> |
3493 | + </conservative_advection> |
3494 | + </spatial_discretisation> |
3495 | + <temporal_discretisation> |
3496 | + <theta> |
3497 | + <real_value rank="0">1.0</real_value> |
3498 | + </theta> |
3499 | + </temporal_discretisation> |
3500 | + <solver> |
3501 | + <iterative_method name="gmres"> |
3502 | + <restart> |
3503 | + <integer_value rank="0">30</integer_value> |
3504 | + </restart> |
3505 | + </iterative_method> |
3506 | + <preconditioner name="sor"/> |
3507 | + <relative_error> |
3508 | + <real_value rank="0">1.0e-7</real_value> |
3509 | + </relative_error> |
3510 | + <max_iterations> |
3511 | + <integer_value rank="0">1000</integer_value> |
3512 | + </max_iterations> |
3513 | + <never_ignore_solver_failures/> |
3514 | + <diagnostics> |
3515 | + <monitors/> |
3516 | + </diagnostics> |
3517 | + </solver> |
3518 | + <initial_condition name="WholeMesh"> |
3519 | + <constant> |
3520 | + <real_value rank="0">0.0</real_value> |
3521 | + </constant> |
3522 | + </initial_condition> |
3523 | + <boundary_conditions name="Top"> |
3524 | + <surface_ids> |
3525 | + <integer_value shape="1" rank="1">444</integer_value> |
3526 | + </surface_ids> |
3527 | + <type name="flux"> |
3528 | + <constant> |
3529 | + <real_value rank="0">0.5</real_value> |
3530 | + </constant> |
3531 | + </type> |
3532 | + </boundary_conditions> |
3533 | + <boundary_conditions name="Others"> |
3534 | + <surface_ids> |
3535 | + <integer_value shape="3" rank="1">333 666 777</integer_value> |
3536 | + </surface_ids> |
3537 | + <type name="zero_flux"/> |
3538 | + </boundary_conditions> |
3539 | + <output/> |
3540 | + <stat/> |
3541 | + <convergence> |
3542 | + <include_in_convergence/> |
3543 | + </convergence> |
3544 | + <detectors> |
3545 | + <include_in_detectors/> |
3546 | + </detectors> |
3547 | + <steady_state> |
3548 | + <include_in_steady_state/> |
3549 | + </steady_state> |
3550 | + <consistent_interpolation/> |
3551 | + </prognostic> |
3552 | + </scalar_field> |
3553 | + </material_phase> |
3554 | +</fluidity_options> |
3555 | |
3556 | === added directory 'tests/flux_bc/src' |
3557 | === added file 'tests/flux_bc/src/flux_bc_2d.geo' |
3558 | --- tests/flux_bc/src/flux_bc_2d.geo 1970-01-01 00:00:00 +0000 |
3559 | +++ tests/flux_bc/src/flux_bc_2d.geo 2012-08-09 00:46:20 +0000 |
3560 | @@ -0,0 +1,19 @@ |
3561 | +Point(1) = {0.0, 0.0, 0.0, 0.05}; |
3562 | +Point(2) = {0.0, 1.0, 0.0, 0.05}; |
3563 | +Point(3) = {0.3, 1.0, 0.0, 0.05}; |
3564 | +Point(4) = {0.3, 0.0, 0.0, 0.05}; |
3565 | +Line(1) = {3,4}; |
3566 | +Line(2) = {4,1}; |
3567 | +Line(3) = {1,2}; |
3568 | +Line(4) = {2,3}; |
3569 | +Line Loop(5) = {1,2,3,4}; |
3570 | +Plane Surface(6) = {5}; |
3571 | +// Top of the box |
3572 | +Physical Line(333) = {4}; |
3573 | +// Box sides |
3574 | +Physical Line(666) = {3,1}; |
3575 | +// Box bottom |
3576 | +Physical Line(444) = {2}; |
3577 | +// This is just to ensure all the interior |
3578 | +// elements get written out. |
3579 | +Physical Surface(10) = {6}; |
3580 | |
3581 | === added file 'tests/flux_bc/src/flux_bc_3d.geo' |
3582 | --- tests/flux_bc/src/flux_bc_3d.geo 1970-01-01 00:00:00 +0000 |
3583 | +++ tests/flux_bc/src/flux_bc_3d.geo 2012-08-09 00:46:20 +0000 |
3584 | @@ -0,0 +1,69 @@ |
3585 | +Point(1) = {0.0, 0.0, 0.0, 0.05}; |
3586 | +Point(2) = {0.3, 0.0, 0.0, 0.05}; |
3587 | +Point(3) = {0.3, 0.0, 0.3, 0.05}; |
3588 | +Point(4) = {0.0, 0.0, 0.3, 0.05}; |
3589 | +Point(5) = {0.0, 1.0, 0.0, 0.05}; |
3590 | +Point(6) = {0.3, 1.0, 0.0, 0.05}; |
3591 | +Point(7) = {0.3, 1.0, 0.3, 0.05}; |
3592 | +Point(8) = {0.0, 1.0, 0.3, 0.05}; |
3593 | + |
3594 | +Line(1) = {1,2}; |
3595 | +Line(2) = {2,3}; |
3596 | +Line(3) = {3,4}; |
3597 | +Line(4) = {4,1}; |
3598 | + |
3599 | +Line(5) = {5,6}; |
3600 | +Line(6) = {6,7}; |
3601 | +Line(7) = {7,8}; |
3602 | +Line(8) = {8,5}; |
3603 | + |
3604 | +Line(9) = {5,1}; |
3605 | +Line(10) = {6,2}; |
3606 | +Line(11) = {7,3}; |
3607 | +Line(12) = {8,4}; |
3608 | + |
3609 | +// Bottom surface |
3610 | +Line Loop(5) = {1,2,3,4}; |
3611 | + |
3612 | +// Side surface, nearest to the x axis |
3613 | +Line Loop(6) = {1,-10,-5,9}; |
3614 | + |
3615 | +// Side surface, furthest from the x axis |
3616 | +Line Loop(7) = {-3,-11,7,12}; |
3617 | + |
3618 | +// Side surface, nearest to the z axis |
3619 | +Line Loop(8) = {-4,-12,8,9}; |
3620 | + |
3621 | +// Side surface, furthest from the z axis |
3622 | +Line Loop(9) = {2,-11,-6,10}; |
3623 | + |
3624 | +// Top surface |
3625 | +Line Loop(10) = {5,6,7,8}; |
3626 | + |
3627 | +Plane Surface(20) = {5}; |
3628 | +Plane Surface(21) = {6}; |
3629 | +Plane Surface(22) = {7}; |
3630 | +Plane Surface(23) = {8}; |
3631 | +Plane Surface(24) = {9}; |
3632 | +Plane Surface(25) = {10}; |
3633 | + |
3634 | +// Bottom of the box |
3635 | +Physical Line(333) = {5}; |
3636 | +// Sides of the box perpendicular to the z axis |
3637 | +Physical Line(666) = {6,7}; |
3638 | +// Sides of the box perpendicular to the x axis |
3639 | +Physical Line(777) = {8,9}; |
3640 | +// Top of the box |
3641 | +Physical Line(444) = {10}; |
3642 | + |
3643 | +// This is just to ensure all the interior |
3644 | +// elements get written out. |
3645 | +Physical Surface(333) = {20}; |
3646 | +Physical Surface(666) = {21,22}; |
3647 | +Physical Surface(777) = {23,24}; |
3648 | +Physical Surface(444) = {25}; |
3649 | + |
3650 | +Surface Loop(800) = {22, 20, 21, 24, 25, 23}; |
3651 | +Volume(900) = {800}; |
3652 | +Physical Volume(999) = {900}; |
3653 | + |
3654 | |
3655 | === added directory 'tests/mphase_dusty_gas_shock_tube' |
3656 | === added file 'tests/mphase_dusty_gas_shock_tube/Makefile' |
3657 | --- tests/mphase_dusty_gas_shock_tube/Makefile 1970-01-01 00:00:00 +0000 |
3658 | +++ tests/mphase_dusty_gas_shock_tube/Makefile 2012-08-09 00:46:20 +0000 |
3659 | @@ -0,0 +1,18 @@ |
3660 | +preprocess: |
3661 | + @echo **********Creating 1D mesh |
3662 | + ../../bin/interval --dx=0.00271002710027 -- -1.35501355014 1.35501355014 line |
3663 | + |
3664 | +run: |
3665 | + @echo **********Running simulation |
3666 | + ../../bin/fluidity -v2 -l mphase_dusty_gas_shock_tube.flml |
3667 | + ../../bin/fluidity -v2 -l single_phase_frozen_flow_test.flml |
3668 | + |
3669 | +input: clean preprocess |
3670 | + |
3671 | +clean: |
3672 | + rm -f *.stat *.steady_state* |
3673 | + rm -f *.d.* *.vtu |
3674 | + rm -f *.msh |
3675 | + rm -f *.ele *.edge *.node *.poly *.bound |
3676 | + rm -f matrixdump* *.log* *.err* |
3677 | + |
3678 | |
3679 | === added file 'tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml' |
3680 | --- tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml 1970-01-01 00:00:00 +0000 |
3681 | +++ tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.flml 2012-08-09 00:46:20 +0000 |
3682 | @@ -0,0 +1,733 @@ |
3683 | +<?xml version='1.0' encoding='utf-8'?> |
3684 | +<fluidity_options> |
3685 | + <simulation_name> |
3686 | + <string_value lines="1">mphase_dusty_gas_shock_tube</string_value> |
3687 | + </simulation_name> |
3688 | + <problem_type> |
3689 | + <string_value lines="1">multiphase</string_value> |
3690 | + </problem_type> |
3691 | + <geometry> |
3692 | + <dimension> |
3693 | + <integer_value rank="0">1</integer_value> |
3694 | + </dimension> |
3695 | + <mesh name="CoordinateMesh"> |
3696 | + <from_file file_name="line"> |
3697 | + <format name="triangle"/> |
3698 | + <stat> |
3699 | + <include_in_stat/> |
3700 | + </stat> |
3701 | + </from_file> |
3702 | + </mesh> |
3703 | + <mesh name="VelocityMesh"> |
3704 | + <from_mesh> |
3705 | + <mesh name="CoordinateMesh"/> |
3706 | + <mesh_shape> |
3707 | + <polynomial_degree> |
3708 | + <integer_value rank="0">1</integer_value> |
3709 | + </polynomial_degree> |
3710 | + </mesh_shape> |
3711 | + <stat> |
3712 | + <exclude_from_stat/> |
3713 | + </stat> |
3714 | + </from_mesh> |
3715 | + </mesh> |
3716 | + <mesh name="PressureMesh"> |
3717 | + <from_mesh> |
3718 | + <mesh name="CoordinateMesh"/> |
3719 | + <mesh_shape> |
3720 | + <polynomial_degree> |
3721 | + <integer_value rank="0">1</integer_value> |
3722 | + </polynomial_degree> |
3723 | + </mesh_shape> |
3724 | + <stat> |
3725 | + <exclude_from_stat/> |
3726 | + </stat> |
3727 | + </from_mesh> |
3728 | + </mesh> |
3729 | + <mesh name="DensityMesh"> |
3730 | + <from_mesh> |
3731 | + <mesh name="CoordinateMesh"/> |
3732 | + <mesh_shape> |
3733 | + <polynomial_degree> |
3734 | + <integer_value rank="0">1</integer_value> |
3735 | + </polynomial_degree> |
3736 | + </mesh_shape> |
3737 | + <stat> |
3738 | + <exclude_from_stat/> |
3739 | + </stat> |
3740 | + </from_mesh> |
3741 | + </mesh> |
3742 | + <quadrature> |
3743 | + <degree> |
3744 | + <integer_value rank="0">5</integer_value> |
3745 | + </degree> |
3746 | + </quadrature> |
3747 | + </geometry> |
3748 | + <io> |
3749 | + <dump_format> |
3750 | + <string_value>vtk</string_value> |
3751 | + </dump_format> |
3752 | + <dump_period> |
3753 | + <constant> |
3754 | + <real_value rank="0">0.0</real_value> |
3755 | + </constant> |
3756 | + </dump_period> |
3757 | + <output_mesh name="VelocityMesh"/> |
3758 | + <stat> |
3759 | + <output_at_start/> |
3760 | + </stat> |
3761 | + </io> |
3762 | + <timestepping> |
3763 | + <current_time> |
3764 | + <real_value rank="0">0.0</real_value> |
3765 | + </current_time> |
3766 | + <timestep> |
3767 | + <real_value rank="0">0.000001</real_value> |
3768 | + </timestep> |
3769 | + <finish_time> |
3770 | + <real_value rank="0">0.00037772998222</real_value> |
3771 | + <comment>This corresponds to a non-dimensional time of tau = 4. Multiplying this by the reference time of l/u_ref gives 0.00037772998222 seconds. |
3772 | + |
3773 | +l = (4/3)*(rho_p/rho_g)*d |
3774 | +u_ref = sqrt(p/rho_g)</comment> |
3775 | + </finish_time> |
3776 | + <nonlinear_iterations> |
3777 | + <integer_value rank="0">3</integer_value> |
3778 | + <tolerance> |
3779 | + <real_value rank="0">1.0e-9</real_value> |
3780 | + <infinity_norm/> |
3781 | + </tolerance> |
3782 | + </nonlinear_iterations> |
3783 | + </timestepping> |
3784 | + <material_phase name="Air"> |
3785 | + <equation_of_state> |
3786 | + <compressible> |
3787 | + <stiffened_gas> |
3788 | + <ratio_specific_heats> |
3789 | + <real_value rank="0">1.4</real_value> |
3790 | + </ratio_specific_heats> |
3791 | + </stiffened_gas> |
3792 | + </compressible> |
3793 | + </equation_of_state> |
3794 | + <scalar_field name="Pressure" rank="0"> |
3795 | + <prognostic> |
3796 | + <mesh name="PressureMesh"/> |
3797 | + <spatial_discretisation> |
3798 | + <continuous_galerkin> |
3799 | + <remove_stabilisation_term/> |
3800 | + </continuous_galerkin> |
3801 | + </spatial_discretisation> |
3802 | + <scheme> |
3803 | + <poisson_pressure_solution> |
3804 | + <string_value lines="1">never</string_value> |
3805 | + </poisson_pressure_solution> |
3806 | + <use_compressible_projection_method/> |
3807 | + </scheme> |
3808 | + <solver> |
3809 | + <iterative_method name="preonly"/> |
3810 | + <preconditioner name="lu"/> |
3811 | + <relative_error> |
3812 | + <real_value rank="0">1e-7</real_value> |
3813 | + </relative_error> |
3814 | + <max_iterations> |
3815 | + <integer_value rank="0">10000</integer_value> |
3816 | + </max_iterations> |
3817 | + <never_ignore_solver_failures/> |
3818 | + <diagnostics> |
3819 | + <monitors/> |
3820 | + </diagnostics> |
3821 | + </solver> |
3822 | + <output/> |
3823 | + <stat/> |
3824 | + <convergence> |
3825 | + <include_in_convergence/> |
3826 | + </convergence> |
3827 | + <detectors> |
3828 | + <exclude_from_detectors/> |
3829 | + </detectors> |
3830 | + <steady_state> |
3831 | + <include_in_steady_state/> |
3832 | + </steady_state> |
3833 | + <no_interpolation/> |
3834 | + </prognostic> |
3835 | + </scalar_field> |
3836 | + <scalar_field name="Density" rank="0"> |
3837 | + <prognostic> |
3838 | + <mesh name="PressureMesh"/> |
3839 | + <spatial_discretisation> |
3840 | + <continuous_galerkin> |
3841 | + <stabilisation> |
3842 | + <no_stabilisation/> |
3843 | + </stabilisation> |
3844 | + <advection_terms/> |
3845 | + <mass_terms/> |
3846 | + </continuous_galerkin> |
3847 | + <conservative_advection> |
3848 | + <real_value rank="0">0.0</real_value> |
3849 | + </conservative_advection> |
3850 | + </spatial_discretisation> |
3851 | + <temporal_discretisation> |
3852 | + <theta> |
3853 | + <real_value rank="0">1.0</real_value> |
3854 | + </theta> |
3855 | + </temporal_discretisation> |
3856 | + <initial_condition name="WholeMesh"> |
3857 | + <python> |
3858 | + <string_value lines="20" type="code" language="python">def val(X,t): |
3859 | + if(X[0] <= 0): |
3860 | + return 12.3 |
3861 | + else: |
3862 | + return 1.23</string_value> |
3863 | + </python> |
3864 | + </initial_condition> |
3865 | + <output/> |
3866 | + <stat/> |
3867 | + <convergence> |
3868 | + <include_in_convergence/> |
3869 | + </convergence> |
3870 | + <detectors> |
3871 | + <include_in_detectors/> |
3872 | + </detectors> |
3873 | + <steady_state> |
3874 | + <include_in_steady_state/> |
3875 | + </steady_state> |
3876 | + <consistent_interpolation/> |
3877 | + </prognostic> |
3878 | + </scalar_field> |
3879 | + <vector_field name="Velocity" rank="1"> |
3880 | + <prognostic> |
3881 | + <mesh name="VelocityMesh"/> |
3882 | + <equation name="LinearMomentum"/> |
3883 | + <spatial_discretisation> |
3884 | + <continuous_galerkin> |
3885 | + <stabilisation> |
3886 | + <streamline_upwind> |
3887 | + <nu_bar_optimal/> |
3888 | + <nu_scale name="0.5"> |
3889 | + <real_value shape="1" rank="0">0.5</real_value> |
3890 | + </nu_scale> |
3891 | + </streamline_upwind> |
3892 | + </stabilisation> |
3893 | + <mass_terms> |
3894 | + <lump_mass_matrix/> |
3895 | + </mass_terms> |
3896 | + <advection_terms/> |
3897 | + <stress_terms> |
3898 | + <tensor_form/> |
3899 | + </stress_terms> |
3900 | + </continuous_galerkin> |
3901 | + <conservative_advection> |
3902 | + <real_value rank="0">0.0</real_value> |
3903 | + </conservative_advection> |
3904 | + </spatial_discretisation> |
3905 | + <temporal_discretisation> |
3906 | + <theta> |
3907 | + <real_value rank="0">1.0</real_value> |
3908 | + </theta> |
3909 | + <relaxation> |
3910 | + <real_value rank="0">0.5</real_value> |
3911 | + </relaxation> |
3912 | + </temporal_discretisation> |
3913 | + <solver> |
3914 | + <iterative_method name="preonly"/> |
3915 | + <preconditioner name="lu"/> |
3916 | + <relative_error> |
3917 | + <real_value rank="0">1e-7</real_value> |
3918 | + </relative_error> |
3919 | + <max_iterations> |
3920 | + <integer_value rank="0">10000</integer_value> |
3921 | + </max_iterations> |
3922 | + <never_ignore_solver_failures/> |
3923 | + <diagnostics> |
3924 | + <monitors/> |
3925 | + </diagnostics> |
3926 | + </solver> |
3927 | + <initial_condition name="WholeMesh"> |
3928 | + <constant> |
3929 | + <real_value shape="1" dim1="dim" rank="1">0.0</real_value> |
3930 | + </constant> |
3931 | + </initial_condition> |
3932 | + <boundary_conditions name="Ends"> |
3933 | + <surface_ids> |
3934 | + <integer_value shape="2" rank="1">1 2</integer_value> |
3935 | + </surface_ids> |
3936 | + <type name="dirichlet"> |
3937 | + <align_bc_with_cartesian> |
3938 | + <x_component> |
3939 | + <constant> |
3940 | + <real_value rank="0">0.0</real_value> |
3941 | + </constant> |
3942 | + </x_component> |
3943 | + </align_bc_with_cartesian> |
3944 | + </type> |
3945 | + </boundary_conditions> |
3946 | + <tensor_field name="Viscosity" rank="2"> |
3947 | + <prescribed> |
3948 | + <value name="WholeMesh"> |
3949 | + <isotropic> |
3950 | + <constant> |
3951 | + <real_value rank="0">1.83e-5</real_value> |
3952 | + </constant> |
3953 | + </isotropic> |
3954 | + </value> |
3955 | + <output/> |
3956 | + </prescribed> |
3957 | + </tensor_field> |
3958 | + <output/> |
3959 | + <stat> |
3960 | + <include_in_stat/> |
3961 | + <previous_time_step> |
3962 | + <exclude_from_stat/> |
3963 | + </previous_time_step> |
3964 | + <nonlinear_field> |
3965 | + <exclude_from_stat/> |
3966 | + </nonlinear_field> |
3967 | + </stat> |
3968 | + <convergence> |
3969 | + <include_in_convergence/> |
3970 | + </convergence> |
3971 | + <detectors> |
3972 | + <include_in_detectors/> |
3973 | + </detectors> |
3974 | + <steady_state> |
3975 | + <include_in_steady_state/> |
3976 | + </steady_state> |
3977 | + <consistent_interpolation/> |
3978 | + </prognostic> |
3979 | + </vector_field> |
3980 | + <scalar_field name="PhaseVolumeFraction" rank="0"> |
3981 | + <diagnostic> |
3982 | + <mesh name="CoordinateMesh"/> |
3983 | + <algorithm name="Internal" material_phase_support="multiple"/> |
3984 | + <output/> |
3985 | + <stat/> |
3986 | + <detectors> |
3987 | + <include_in_detectors/> |
3988 | + </detectors> |
3989 | + </diagnostic> |
3990 | + </scalar_field> |
3991 | + <scalar_field name="CompressibleContinuityResidual" rank="0"> |
3992 | + <diagnostic> |
3993 | + <mesh name="PressureMesh"/> |
3994 | + <algorithm name="Internal" material_phase_support="multiple"/> |
3995 | + <output/> |
3996 | + <stat/> |
3997 | + <detectors> |
3998 | + <include_in_detectors/> |
3999 | + </detectors> |
4000 | + <solver> |
4001 | + <iterative_method name="gmres"> |
4002 | + <restart> |
4003 | + <integer_value rank="0">30</integer_value> |
4004 | + </restart> |
4005 | + </iterative_method> |
4006 | + <preconditioner name="sor"/> |
4007 | + <relative_error> |
4008 | + <real_value rank="0">1.0e-7</real_value> |
4009 | + </relative_error> |
4010 | + <max_iterations> |
4011 | + <integer_value rank="0">10000</integer_value> |
4012 | + </max_iterations> |
4013 | + <never_ignore_solver_failures/> |
4014 | + <diagnostics> |
4015 | + <monitors/> |
4016 | + </diagnostics> |
4017 | + </solver> |
4018 | + </diagnostic> |
4019 | + </scalar_field> |
4020 | + <scalar_field name="InternalEnergy" rank="0"> |
4021 | + <prognostic> |
4022 | + <mesh name="PressureMesh"/> |
4023 | + <equation name="InternalEnergy"> |
4024 | + <density name="Density"/> |
4025 | + </equation> |
4026 | + <spatial_discretisation> |
4027 | + <continuous_galerkin> |
4028 | + <stabilisation> |
4029 | + <streamline_upwind> |
4030 | + <nu_bar_optimal/> |
4031 | + <nu_scale name="0.5"> |
4032 | + <real_value shape="1" rank="0">0.5</real_value> |
4033 | + </nu_scale> |
4034 | + </streamline_upwind> |
4035 | + </stabilisation> |
4036 | + <advection_terms/> |
4037 | + <mass_terms/> |
4038 | + </continuous_galerkin> |
4039 | + <conservative_advection> |
4040 | + <real_value rank="0">0.0</real_value> |
4041 | + </conservative_advection> |
4042 | + </spatial_discretisation> |
4043 | + <temporal_discretisation> |
4044 | + <theta> |
4045 | + <real_value rank="0">1.0</real_value> |
4046 | + </theta> |
4047 | + </temporal_discretisation> |
4048 | + <solver> |
4049 | + <iterative_method name="preonly"/> |
4050 | + <preconditioner name="lu"/> |
4051 | + <relative_error> |
4052 | + <real_value rank="0">1.0e-7</real_value> |
4053 | + </relative_error> |
4054 | + <max_iterations> |
4055 | + <integer_value rank="0">10000</integer_value> |
4056 | + </max_iterations> |
4057 | + <never_ignore_solver_failures/> |
4058 | + <diagnostics> |
4059 | + <monitors/> |
4060 | + </diagnostics> |
4061 | + </solver> |
4062 | + <initial_condition name="WholeMesh"> |
4063 | + <constant> |
4064 | + <real_value rank="0">205894.308943</real_value> |
4065 | + <comment>InternalEnergy from the formula: (3/2)*R*Temperature = (3/2)*287.04*300 K = 203252 |
4066 | +InternalEnergy needed for a Pressure of 1.013e5 = 205894.308943</comment> |
4067 | + </constant> |
4068 | + </initial_condition> |
4069 | + <output/> |
4070 | + <stat/> |
4071 | + <convergence> |
4072 | + <include_in_convergence/> |
4073 | + </convergence> |
4074 | + <detectors> |
4075 | + <include_in_detectors/> |
4076 | + </detectors> |
4077 | + <steady_state> |
4078 | + <include_in_steady_state/> |
4079 | + </steady_state> |
4080 | + <consistent_interpolation/> |
4081 | + </prognostic> |
4082 | + </scalar_field> |
4083 | + <scalar_field name="CFLNumber" rank="0"> |
4084 | + <diagnostic> |
4085 | + <algorithm name="Internal" material_phase_support="multiple"/> |
4086 | + <mesh name="VelocityMesh"/> |
4087 | + <output/> |
4088 | + <stat/> |
4089 | + <convergence> |
4090 | + <include_in_convergence/> |
4091 | + </convergence> |
4092 | + <detectors> |
4093 | + <include_in_detectors/> |
4094 | + </detectors> |
4095 | + <steady_state> |
4096 | + <include_in_steady_state/> |
4097 | + </steady_state> |
4098 | + </diagnostic> |
4099 | + </scalar_field> |
4100 | + <multiphase_properties> |
4101 | + <effective_conductivity> |
4102 | + <real_value rank="0">0.026</real_value> |
4103 | + </effective_conductivity> |
4104 | + <specific_heat> |
4105 | + <real_value rank="0">758</real_value> |
4106 | + </specific_heat> |
4107 | + </multiphase_properties> |
4108 | + </material_phase> |
4109 | + <material_phase name="Dust"> |
4110 | + <equation_of_state> |
4111 | + <fluids> |
4112 | + <linear> |
4113 | + <reference_density> |
4114 | + <real_value rank="0">2500</real_value> |
4115 | + </reference_density> |
4116 | + </linear> |
4117 | + </fluids> |
4118 | + </equation_of_state> |
4119 | + <scalar_field name="Pressure" rank="0"> |
4120 | + <aliased material_phase_name="Air" field_name="Pressure"/> |
4121 | + </scalar_field> |
4122 | + <scalar_field name="Density" rank="0"> |
4123 | + <diagnostic> |
4124 | + <algorithm name="Internal" material_phase_support="multiple"/> |
4125 | + <mesh name="PressureMesh"/> |
4126 | + <output/> |
4127 | + <stat/> |
4128 | + <convergence> |
4129 | + <include_in_convergence/> |
4130 | + </convergence> |
4131 | + <detectors> |
4132 | + <include_in_detectors/> |
4133 | + </detectors> |
4134 | + <steady_state> |
4135 | + <include_in_steady_state/> |
4136 | + </steady_state> |
4137 | + </diagnostic> |
4138 | + </scalar_field> |
4139 | + <vector_field name="Velocity" rank="1"> |
4140 | + <prognostic> |
4141 | + <mesh name="VelocityMesh"/> |
4142 | + <equation name="LinearMomentum"/> |
4143 | + <spatial_discretisation> |
4144 | + <continuous_galerkin> |
4145 | + <stabilisation> |
4146 | + <no_stabilisation/> |
4147 | + </stabilisation> |
4148 | + <mass_terms> |
4149 | + <lump_mass_matrix/> |
4150 | + </mass_terms> |
4151 | + <advection_terms/> |
4152 | + <stress_terms> |
4153 | + <tensor_form/> |
4154 | + </stress_terms> |
4155 | + </continuous_galerkin> |
4156 | + <conservative_advection> |
4157 | + <real_value rank="0">0.0</real_value> |
4158 | + </conservative_advection> |
4159 | + </spatial_discretisation> |
4160 | + <temporal_discretisation> |
4161 | + <theta> |
4162 | + <real_value rank="0">1.0</real_value> |
4163 | + </theta> |
4164 | + <relaxation> |
4165 | + <real_value rank="0">0.5</real_value> |
4166 | + </relaxation> |
4167 | + </temporal_discretisation> |
4168 | + <solver> |
4169 | + <iterative_method name="preonly"/> |
4170 | + <preconditioner name="lu"/> |
4171 | + <relative_error> |
4172 | + <real_value rank="0">1e-7</real_value> |
4173 | + </relative_error> |
4174 | + <max_iterations> |
4175 | + <integer_value rank="0">10000</integer_value> |
4176 | + </max_iterations> |
4177 | + <never_ignore_solver_failures/> |
4178 | + <diagnostics> |
4179 | + <monitors/> |
4180 | + </diagnostics> |
4181 | + </solver> |
4182 | + <initial_condition name="WholeMesh"> |
4183 | + <constant> |
4184 | + <real_value shape="1" dim1="dim" rank="1">0.0</real_value> |
4185 | + </constant> |
4186 | + </initial_condition> |
4187 | + <boundary_conditions name="Ends"> |
4188 | + <surface_ids> |
4189 | + <integer_value shape="2" rank="1">1 2</integer_value> |
4190 | + </surface_ids> |
4191 | + <type name="dirichlet"> |
4192 | + <align_bc_with_cartesian> |
4193 | + <x_component> |
4194 | + <constant> |
4195 | + <real_value rank="0">0.0</real_value> |
4196 | + </constant> |
4197 | + </x_component> |
4198 | + </align_bc_with_cartesian> |
4199 | + </type> |
4200 | + </boundary_conditions> |
4201 | + <output/> |
4202 | + <stat> |
4203 | + <include_in_stat/> |
4204 | + <previous_time_step> |
4205 | + <exclude_from_stat/> |
4206 | + </previous_time_step> |
4207 | + <nonlinear_field> |
4208 | + <exclude_from_stat/> |
4209 | + </nonlinear_field> |
4210 | + </stat> |
4211 | + <convergence> |
4212 | + <include_in_convergence/> |
4213 | + </convergence> |
4214 | + <detectors> |
4215 | + <include_in_detectors/> |
4216 | + </detectors> |
4217 | + <steady_state> |
4218 | + <include_in_steady_state/> |
4219 | + </steady_state> |
4220 | + <consistent_interpolation/> |
4221 | + </prognostic> |
4222 | + </vector_field> |
4223 | + <scalar_field name="PhaseVolumeFraction" rank="0"> |
4224 | + <prognostic> |
4225 | + <mesh name="CoordinateMesh"/> |
4226 | + <equation name="AdvectionDiffusion"/> |
4227 | + <spatial_discretisation> |
4228 | + <control_volumes> |
4229 | + <face_value name="FiniteElement"> |
4230 | + <limit_face_value> |
4231 | + <limiter name="Sweby"/> |
4232 | + </limit_face_value> |
4233 | + </face_value> |
4234 | + <diffusion_scheme name="BassiRebay"/> |
4235 | + </control_volumes> |
4236 | + <conservative_advection> |
4237 | + <real_value rank="0">1.0</real_value> |
4238 | + </conservative_advection> |
4239 | + </spatial_discretisation> |
4240 | + <temporal_discretisation> |
4241 | + <theta> |
4242 | + <real_value rank="0">1.0</real_value> |
4243 | + </theta> |
4244 | + </temporal_discretisation> |
4245 | + <solver> |
4246 | + <iterative_method name="preonly"/> |
4247 | + <preconditioner name="lu"/> |
4248 | + <relative_error> |
4249 | + <real_value rank="0">1.0e-7</real_value> |
4250 | + </relative_error> |
4251 | + <max_iterations> |
4252 | + <integer_value rank="0">10000</integer_value> |
4253 | + </max_iterations> |
4254 | + <never_ignore_solver_failures/> |
4255 | + <diagnostics> |
4256 | + <monitors/> |
4257 | + </diagnostics> |
4258 | + </solver> |
4259 | + <initial_condition name="WholeMesh"> |
4260 | + <python> |
4261 | + <string_value lines="20" type="code" language="python">def val(X,t): |
4262 | + if(X[0] <= 0): |
4263 | + return 4.92e-8 |
4264 | + else: |
4265 | + return 4.92e-4</string_value> |
4266 | + <comment>alpha_p = 0.000492 in dusty gas region. |
4267 | + |
4268 | +We know that the ratio of gas density (rho_g) to particle mass concentration (sigma_p) is assumed to be 1.0, so rho_g = sigma_p = 1.23 kg/m^3. |
4269 | +The normalised mass concentration is 1.0 in the dusty gas region, so we divide 1.23 by the particle density (2500 kg/m^3) get obtain the particle phase volume fraction.</comment> |
4270 | + </python> |
4271 | + </initial_condition> |
4272 | + <boundary_conditions name="All"> |
4273 | + <surface_ids> |
4274 | + <integer_value shape="2" rank="1">1 2</integer_value> |
4275 | + </surface_ids> |
4276 | + <type name="zero_flux"/> |
4277 | + </boundary_conditions> |
4278 | + <output/> |
4279 | + <stat/> |
4280 | + <convergence> |
4281 | + <include_in_convergence/> |
4282 | + </convergence> |
4283 | + <detectors> |
4284 | + <include_in_detectors/> |
4285 | + </detectors> |
4286 | + <steady_state> |
4287 | + <include_in_steady_state/> |
4288 | + </steady_state> |
4289 | + <consistent_interpolation/> |
4290 | + </prognostic> |
4291 | + </scalar_field> |
4292 | + <scalar_field name="ParticleReynoldsNumber" rank="0"> |
4293 | + <diagnostic> |
4294 | + <algorithm name="particle_reynolds_number" material_phase_support="multiple"> |
4295 | + <depends> |
4296 | + <string_value lines="1">Dust::Velocity,Air::Velocity,Air::PhaseVolumeFraction,Air::Density</string_value> |
4297 | + </depends> |
4298 | + <particle_diameter> |
4299 | + <real_value rank="0">10e-6</real_value> |
4300 | + </particle_diameter> |
4301 | + <continuous_phase_name>Air</continuous_phase_name> |
4302 | + </algorithm> |
4303 | + <mesh name="VelocityMesh"/> |
4304 | + <output/> |
4305 | + <stat/> |
4306 | + <convergence> |
4307 | + <include_in_convergence/> |
4308 | + </convergence> |
4309 | + <detectors> |
4310 | + <include_in_detectors/> |
4311 | + </detectors> |
4312 | + <steady_state> |
4313 | + <include_in_steady_state/> |
4314 | + </steady_state> |
4315 | + </diagnostic> |
4316 | + </scalar_field> |
4317 | + <scalar_field name="CFLNumber" rank="0"> |
4318 | + <diagnostic> |
4319 | + <algorithm name="Internal" material_phase_support="multiple"/> |
4320 | + <mesh name="VelocityMesh"/> |
4321 | + <output/> |
4322 | + <stat/> |
4323 | + <convergence> |
4324 | + <include_in_convergence/> |
4325 | + </convergence> |
4326 | + <detectors> |
4327 | + <include_in_detectors/> |
4328 | + </detectors> |
4329 | + <steady_state> |
4330 | + <include_in_steady_state/> |
4331 | + </steady_state> |
4332 | + </diagnostic> |
4333 | + </scalar_field> |
4334 | + <scalar_field name="InternalEnergy" rank="0"> |
4335 | + <prognostic> |
4336 | + <mesh name="PressureMesh"/> |
4337 | + <equation name="InternalEnergy"> |
4338 | + <density name="Density"/> |
4339 | + </equation> |
4340 | + <spatial_discretisation> |
4341 | + <continuous_galerkin> |
4342 | + <stabilisation> |
4343 | + <streamline_upwind> |
4344 | + <nu_bar_optimal/> |
4345 | + <nu_scale name="0.5"> |
4346 | + <real_value shape="1" rank="0">0.5</real_value> |
4347 | + </nu_scale> |
4348 | + </streamline_upwind> |
4349 | + </stabilisation> |
4350 | + <advection_terms/> |
4351 | + <mass_terms/> |
4352 | + </continuous_galerkin> |
4353 | + <conservative_advection> |
4354 | + <real_value rank="0">0.0</real_value> |
4355 | + </conservative_advection> |
4356 | + </spatial_discretisation> |
4357 | + <temporal_discretisation> |
4358 | + <theta> |
4359 | + <real_value rank="0">1.0</real_value> |
4360 | + </theta> |
4361 | + </temporal_discretisation> |
4362 | + <solver> |
4363 | + <iterative_method name="preonly"/> |
4364 | + <preconditioner name="lu"/> |
4365 | + <relative_error> |
4366 | + <real_value rank="0">1.0e-7</real_value> |
4367 | + </relative_error> |
4368 | + <max_iterations> |
4369 | + <integer_value rank="0">10000</integer_value> |
4370 | + </max_iterations> |
4371 | + <never_ignore_solver_failures/> |
4372 | + <diagnostics> |
4373 | + <monitors/> |
4374 | + </diagnostics> |
4375 | + </solver> |
4376 | + <initial_condition name="WholeMesh"> |
4377 | + <constant> |
4378 | + <real_value rank="0">205894.308943</real_value> |
4379 | + <comment>InternalEnergy from the formula: (3/2)*R*Temperature = (3/2)*287.04*300 K = 203252 |
4380 | +InternalEnergy needed for a Pressure of 1.013e5 = 205894.308943</comment> |
4381 | + </constant> |
4382 | + </initial_condition> |
4383 | + <output/> |
4384 | + <stat/> |
4385 | + <convergence> |
4386 | + <include_in_convergence/> |
4387 | + </convergence> |
4388 | + <detectors> |
4389 | + <include_in_detectors/> |
4390 | + </detectors> |
4391 | + <steady_state> |
4392 | + <include_in_steady_state/> |
4393 | + </steady_state> |
4394 | + <consistent_interpolation/> |
4395 | + </prognostic> |
4396 | + </scalar_field> |
4397 | + <multiphase_properties> |
4398 | + <particle_diameter> |
4399 | + <real_value rank="0">10e-6</real_value> |
4400 | + </particle_diameter> |
4401 | + <specific_heat> |
4402 | + <real_value rank="0">758</real_value> |
4403 | + <comment>1406</comment> |
4404 | + </specific_heat> |
4405 | + </multiphase_properties> |
4406 | + </material_phase> |
4407 | + <multiphase_interaction> |
4408 | + <fluid_particle_drag> |
4409 | + <drag_correlation name="wen_yu"/> |
4410 | + </fluid_particle_drag> |
4411 | + <heat_transfer> |
4412 | + <heat_transfer_coefficient name="gunn"/> |
4413 | + </heat_transfer> |
4414 | + </multiphase_interaction> |
4415 | +</fluidity_options> |
4416 | |
4417 | === added file 'tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.xml' |
4418 | --- tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.xml 1970-01-01 00:00:00 +0000 |
4419 | +++ tests/mphase_dusty_gas_shock_tube/mphase_dusty_gas_shock_tube.xml 2012-08-09 00:46:20 +0000 |
4420 | @@ -0,0 +1,156 @@ |
4421 | +<?xml version="1.0" encoding="UTF-8" ?> |
4422 | + |
4423 | +<testproblem> |
4424 | + |
4425 | + <!-- A 1D gas-solid shock tube problem. The setup is the same as the one given in the paper by Miura & Glass (1982). --> |
4426 | + <!-- This test uses results at approximately t = 3.77e-4 seconds (normalised time of tau = 4 in the paper by Miura & Glass (1982)). --> |
4427 | + |
4428 | + <name>mphase_dusty_gas_shock_tube</name> |
4429 | + <owner userid="ctj10"/> |
4430 | + <tags>flml</tags> |
4431 | + |
4432 | + <problem_definition length="medium" nprocs="1"> |
4433 | + <command_line>make run</command_line> |
4434 | + </problem_definition> |
4435 | + |
4436 | + <variables> |
4437 | + <variable name="time" language="python"> |
4438 | +from fluidity_tools import stat_parser |
4439 | +s = stat_parser("mphase_dusty_gas_shock_tube.stat") |
4440 | +time=s["ElapsedTime"]["value"] |
4441 | + </variable> |
4442 | + |
4443 | + <variable name="positions" language="python"> |
4444 | +import vtktools |
4445 | + |
4446 | +file = vtktools.vtu('mphase_dusty_gas_shock_tube_378.vtu') |
4447 | +file.GetFieldNames() |
4448 | +p = file.GetLocations() |
4449 | +positions = p[:,1] |
4450 | + </variable> |
4451 | + |
4452 | + <variable name="air_velocity" language="python"> |
4453 | +import vtktools |
4454 | + |
4455 | +file = vtktools.vtu('mphase_dusty_gas_shock_tube_378.vtu') |
4456 | +file.GetFieldNames() |
4457 | +air_velocity = file.GetVectorField('Air::Velocity') |
4458 | + </variable> |
4459 | + |
4460 | + <variable name="dust_vfrac" language="python"> |
4461 | +import vtktools |
4462 | + |
4463 | +file = vtktools.vtu('mphase_dusty_gas_shock_tube_378.vtu') |
4464 | +file.GetFieldNames() |
4465 | +dust_vfrac = file.GetScalarField('Dust::PhaseVolumeFraction') |
4466 | + </variable> |
4467 | + |
4468 | + <variable name="pressure" language="python"> |
4469 | +import vtktools |
4470 | + |
4471 | +file = vtktools.vtu('mphase_dusty_gas_shock_tube_378.vtu') |
4472 | +file.GetFieldNames() |
4473 | +pressure = file.GetScalarField('Air::Pressure') |
4474 | + </variable> |
4475 | + |
4476 | + <variable name="max_dust_density" language="python"> |
4477 | +from fluidity_tools import stat_parser |
4478 | +s = stat_parser("mphase_dusty_gas_shock_tube.stat") |
4479 | +max_dust_density=s["Dust"]["Density"]["max"][-1] |
4480 | + </variable> |
4481 | + |
4482 | + <variable name="min_dust_density" language="python"> |
4483 | +from fluidity_tools import stat_parser |
4484 | +s = stat_parser("mphase_dusty_gas_shock_tube.stat") |
4485 | +min_dust_density=s["Dust"]["Density"]["min"][-1] |
4486 | + </variable> |
4487 | + </variables> |
4488 | + |
4489 | + <pass_tests> |
4490 | + <test name="Solvers have converged" language="python"> |
4491 | +import os |
4492 | +files = os.listdir("./") |
4493 | +solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files |
4494 | +assert solvers_converged |
4495 | + </test> |
4496 | + |
4497 | + <test name="Finish time is 3.77e-4 seconds" language="python"> |
4498 | +assert (time[-1] >= 3.77e-4) |
4499 | +assert (time[-1] < (3.77e-4 + 0.000002)) |
4500 | + </test> |
4501 | + |
4502 | + <test name="Normalised Gas::Velocity is 0 for normalised x in [-50, -10]" language="python"> |
4503 | +nnodes = len(air_velocity) |
4504 | + |
4505 | +vals = [] |
4506 | +for i in range(nnodes): |
4507 | + if(positions[i]/0.0271002710027 < -10.0): |
4508 | + assert abs(air_velocity[i,0]/286.980353992) < 1.0e-12 |
4509 | + </test> |
4510 | + |
4511 | + <test name="Normalised Gas::Velocity is in [0, 0.9] for normalised x in [-10, 10]" language="python"> |
4512 | +nnodes = len(air_velocity) |
4513 | + |
4514 | +vals = [] |
4515 | +for i in range(nnodes): |
4516 | + if( (positions[i]/0.0271002710027 > -10.0) and (positions[i]/0.0271002710027 < 10.0) ): |
4517 | + assert abs(air_velocity[i,0]/286.980353992) < 0.9 |
4518 | + <!-- -1.0e-9 is needed to deal with small under-shoots --> |
4519 | + assert abs(air_velocity[i,0]/286.980353992) > -1.0e-12 |
4520 | + </test> |
4521 | + |
4522 | + <test name="Normalised Gas::Velocity is 0 for normalised x in [10, 50]" language="python"> |
4523 | +nnodes = len(air_velocity) |
4524 | + |
4525 | +vals = [] |
4526 | +for i in range(nnodes): |
4527 | + if(positions[i]/0.0271002710027 > 10.0): |
4528 | + assert abs(air_velocity[i,0]/286.980353992) < 1.0e-12 |
4529 | + </test> |
4530 | + |
4531 | + <test name="Normalised Pressure is 10.0 for normalised x in [-50, -10]" language="python"> |
4532 | +nnodes = len(pressure) |
4533 | + |
4534 | +vals = [] |
4535 | +for i in range(nnodes): |
4536 | + if(positions[i]/0.0271002710027 < -10.0): |
4537 | + assert abs(pressure[i]/1.013e5 - 10.0) < 1.0e-6 |
4538 | + </test> |
4539 | + |
4540 | + <test name="Normalised Pressure is in [1.0, 10.0] for normalised x in [-10, 10]" language="python"> |
4541 | +nnodes = len(pressure) |
4542 | +<!-- Note: values of 0.9 and 10.1 have been used instead of 1.0 and 10.0 because of small over-shoots/under-shoots --> |
4543 | +vals = [] |
4544 | +for i in range(nnodes): |
4545 | + if( (positions[i]/0.0271002710027 > -10.0) and (positions[i]/0.0271002710027 < 10.0) ): |
4546 | + assert abs(pressure[i]/1.013e5) < 10.1 |
4547 | + assert abs(pressure[i]/1.013e5) > 0.9 |
4548 | + </test> |
4549 | + |
4550 | + <test name="Normalised Pressure is 1.0 for normalised x in [10, 50]" language="python"> |
4551 | +nnodes = len(air_velocity) |
4552 | + |
4553 | +vals = [] |
4554 | +for i in range(nnodes): |
4555 | + if(positions[i]/0.0271002710027 > 10.0): |
4556 | + assert abs(pressure[i]/1.013e5) >= 1.0 |
4557 | + </test> |
4558 | + |
4559 | + <test name="Mass concentration of Dust is between 0 and 2.5" language="python"> |
4560 | +nnodes = len(air_velocity) |
4561 | + |
4562 | +vals = [] |
4563 | +for i in range(nnodes): |
4564 | + assert abs(dust_vfrac[i]*2500.0/1.23) < 2.5 |
4565 | + assert abs(dust_vfrac[i]*2500.0/1.23) > 0 |
4566 | + </test> |
4567 | + |
4568 | + <test name="Dust density stays constant" language="python"> |
4569 | +assert abs(max_dust_density - min_dust_density) < 1.0e-16 |
4570 | + </test> |
4571 | + </pass_tests> |
4572 | + |
4573 | + <warn_tests> |
4574 | + </warn_tests> |
4575 | + |
4576 | +</testproblem> |
4577 | |
4578 | === added file 'tests/mphase_dusty_gas_shock_tube/plot.py' |
4579 | --- tests/mphase_dusty_gas_shock_tube/plot.py 1970-01-01 00:00:00 +0000 |
4580 | +++ tests/mphase_dusty_gas_shock_tube/plot.py 2012-08-09 00:46:20 +0000 |
4581 | @@ -0,0 +1,93 @@ |
4582 | +import pylab |
4583 | +import numpy |
4584 | +import shocktube |
4585 | +import vtktools |
4586 | +params = {'text.fontsize': 6, |
4587 | + 'legend.fontsize': 8, |
4588 | + 'xtick.labelsize': 8, |
4589 | + 'ytick.labelsize': 8, |
4590 | + 'lines.markersize': 6, |
4591 | + 'axes.titlesize': 'small'} |
4592 | +pylab.rcParams.update(params) |
4593 | + |
4594 | +filename='mphase_dusty_gas_shock_tube_378.vtu' # This corresponds to approximately t = 3.77e-4 seconds (normalised time of tau = 4 in the paper by Miura & Glass (1982)) |
4595 | +vt=vtktools.vtu(filename) |
4596 | + |
4597 | +# Time and coordinates |
4598 | +t=vt.GetScalarField('Air::Time')[0] |
4599 | +xyz=vt.GetLocations() |
4600 | +x=xyz[:,0] |
4601 | + |
4602 | +# Solution fields |
4603 | +p=vt.GetScalarField('Air::Pressure') |
4604 | +uvw_air=vt.GetVectorField('Air::Velocity') |
4605 | +u_air=uvw_air[:,0] |
4606 | +uvw_dust=vt.GetVectorField('Dust::Velocity') |
4607 | +u_dust=uvw_dust[:,0] |
4608 | +rho_air=vt.GetScalarField('Air::Density') |
4609 | +rho_dust=vt.GetScalarField('Dust::Density') |
4610 | +ie_air=vt.GetScalarField('Air::InternalEnergy') |
4611 | +ie_dust=vt.GetScalarField('Dust::InternalEnergy') |
4612 | +vfrac=vt.GetScalarField('Dust::PhaseVolumeFraction') |
4613 | + |
4614 | +# Plot the numerical results. |
4615 | +# Note that we have divided the results by the reference velocity/pressure/density/internal energy to normalise them. |
4616 | +# x is normalised by the reference length. |
4617 | +pylab.figure(1) |
4618 | +pylab.subplot(321) |
4619 | +pylab.plot( x/0.0271002710027, p/1.013e5,'b.', label="Numerical") |
4620 | + |
4621 | +pylab.subplot(322) |
4622 | +print len(x) |
4623 | +print len(u_air) |
4624 | +pylab.plot( x[:len(x)]/0.0271002710027, u_air/286.980353992,'b.', label="Numerical (Air)") |
4625 | +pylab.plot( x[:len(x)]/0.0271002710027, u_dust/286.980353992,'r.', label="Numerical (Particles)") |
4626 | + |
4627 | +pylab.subplot(323) |
4628 | +pylab.plot( x/0.0271002710027, rho_air/1.23,'b.', label="Numerical") |
4629 | + |
4630 | +pylab.subplot(324) |
4631 | +pylab.plot( x/0.0271002710027, rho_dust,'b.', label="Numerical") |
4632 | +pylab.title('Density of Particles') |
4633 | +pylab.legend(loc=2) |
4634 | + |
4635 | +pylab.subplot(325) |
4636 | +pylab.plot( x/0.0271002710027, ie_air/82357.7235772,'b.', label="Numerical (Air)") |
4637 | +pylab.plot( x/0.0271002710027, ie_dust/82357.7235772,'r.', label="Numerical (Dust)") |
4638 | + |
4639 | +# Multiply the PhaseVolumeFraction by the Dust Density to get the mass concentration, and then normalise. |
4640 | +pylab.subplot(326) |
4641 | +pylab.plot( x/0.0271002710027, vfrac*2500.0/1.23,'b.', label="Numerical") |
4642 | +pylab.title('Mass Concentration from PhaseVolumeFraction of Particles') |
4643 | +pylab.legend(loc=2) |
4644 | + |
4645 | +# Now work out the frozen flow ('analytical') solutions |
4646 | +sol=numpy.array([shocktube.solution(xi,4.0) for xi in x/0.0271002710027]) |
4647 | +p=sol[:,0] |
4648 | +u=sol[:,1] |
4649 | +rho=sol[:,2] |
4650 | +ie=p/rho/(shocktube.gamma-1.0) |
4651 | + |
4652 | +pylab.subplot(321) |
4653 | +pylab.plot( x/0.0271002710027, p,'-g', label='Frozen flow') |
4654 | +pylab.title('Normalised Pressure') |
4655 | +pylab.legend(loc=2) |
4656 | + |
4657 | +pylab.subplot(322) |
4658 | +pylab.plot( x/0.0271002710027, u,'g-', label='Frozen flow (Air)') |
4659 | +pylab.title('Normalised Velocity') |
4660 | +pylab.legend(loc=2) |
4661 | + |
4662 | +pylab.subplot(323) |
4663 | +pylab.plot( x/0.0271002710027, rho,'g-', label='Frozen flow') |
4664 | +pylab.title('Normalised Density of Air') |
4665 | +pylab.legend(loc=2) |
4666 | + |
4667 | +pylab.subplot(325) |
4668 | +pylab.plot( x/0.0271002710027, ie,'g-', label='Frozen flow') |
4669 | +pylab.title('Normalised Internal Energy of Air') |
4670 | +pylab.legend(loc=2) |
4671 | + |
4672 | +pylab.savefig('mphase_dusty_gas_shock_tube.png') |
4673 | + |
4674 | +pylab.show() |
4675 | |
4676 | === added file 'tests/mphase_dusty_gas_shock_tube/shocktube.py' |
4677 | --- tests/mphase_dusty_gas_shock_tube/shocktube.py 1970-01-01 00:00:00 +0000 |
4678 | +++ tests/mphase_dusty_gas_shock_tube/shocktube.py 2012-08-09 00:46:20 +0000 |
4679 | @@ -0,0 +1,157 @@ |
4680 | +# bit of python to compute analytical solution for the shock tube problem |
4681 | +# solution taken from "Mathematical and Computational Methods for Compressible Flow" |
4682 | +# by Miloslav Feistauer, Jiri Felcman and Ivan Straskraba, 2003 |
4683 | +# eqn numbers below refer to this books |
4684 | +from math import sqrt |
4685 | +from scipy.optimize import newton |
4686 | + |
4687 | +# ratio of specific heats |
4688 | +gamma=1.4 |
4689 | + |
4690 | +# left state: |
4691 | +ul=0. |
4692 | +rhol=10.0 |
4693 | +iel=2.5 |
4694 | + |
4695 | +# right state: |
4696 | +ur=0. |
4697 | +rhor=1.0 |
4698 | +ier=2.5 |
4699 | + |
4700 | +# this eos is assumed in the eqns below (so can't change this): |
4701 | +def p_eos(ie, rho): |
4702 | + return rho*ie*(gamma-1.0) |
4703 | + |
4704 | +def rho_eos(ie, p): |
4705 | + return p/ie/(gamma-1.0) |
4706 | + |
4707 | + |
4708 | +# derived quantities, left: |
4709 | +pl=p_eos(iel,rhol) |
4710 | +El=rhol*iel |
4711 | +al=sqrt(gamma*pl/rhol) |
4712 | +# derived quantities, right: |
4713 | +pr=p_eos(ier,rhor) |
4714 | +Er=rhor*ier |
4715 | +ar=sqrt(gamma*pr/rhor) |
4716 | + |
4717 | + |
4718 | +def F1l(p): |
4719 | + """Function F1l defines difference between us and ul: |
4720 | + us=ul+F1l(p) eqn. (3.1.159)""" |
4721 | + |
4722 | + # eqn. (3.1.160) |
4723 | + if p<=pl: |
4724 | + return 2.*al/(gamma-1.)*(1.-(p/pl)**((gamma-1.)/(2.*gamma))) |
4725 | + else: |
4726 | + return -(p-pl)*sqrt((2./(gamma+1.)/rhol)/(p+(gamma-1.)/(gamma+1.)*pl)) |
4727 | + |
4728 | +def F1lprime(p): |
4729 | + # derivative of the above |
4730 | + if p<=pl: |
4731 | + return -al/gamma*(p/pl)**((gamma-1.)/(2.*gamma)-1.) |
4732 | + else: |
4733 | + return -1.*sqrt((2./(gamma+1.)/rhol)/(p+(gamma-1.)/(gamma+1.)*pl)) + \ |
4734 | + -(p-pl)/2./sqrt((2./(gamma+1.)/rhol)/(p+(gamma-1.)/(gamma+1.)*pl)) |
4735 | + |
4736 | +def F3r(p): |
4737 | + """Function F3r defines difference between us and ur: |
4738 | + us=ur+F3r(p) eqn. (3.1.161)""" |
4739 | + |
4740 | + # eqn. (3.1.162) |
4741 | + if p<=pr: |
4742 | + return -2.*ar/(gamma-1.)*(1.-(p/pr)**((gamma-1.)/(2.*gamma))) |
4743 | + else: |
4744 | + return (p-pr)*sqrt((2./(gamma+1.)/rhor)/(p+(gamma-1.)/(gamma+1.)*pr)) |
4745 | + |
4746 | +def F3rprime(p): |
4747 | + # derivative of the above |
4748 | + if p<=pr: |
4749 | + return ar/gamma*(p/pr)**((gamma-1.)/(2.*gamma)-1.) |
4750 | + else: |
4751 | + return sqrt((2./(gamma+1.)/rhor)/(p+(gamma-1.)/(gamma+1.)*pr)) + \ |
4752 | + (p-pr)/2./sqrt((2./(gamma+1.)/rhor)/(p+(gamma-1.)/(gamma+1.)*pr)) |
4753 | + |
4754 | +def F(p): |
4755 | + # eqn (3.1.165) |
4756 | + return F3r(p)-F1l(p)+ur-ul |
4757 | + |
4758 | +def Fprime(p): |
4759 | + return F3rprime(p)-F1lprime(p) |
4760 | + |
4761 | +# inital guess: |
4762 | +p=(pl+pr)/2. |
4763 | + |
4764 | +# solve eqn (3.1.164): F(p*)=0, to find p* the pressure between the u-a and u+a waves |
4765 | +# (is constant over contact discontinuity at the u wave) |
4766 | +ps=newton(F, p, fprime=Fprime, maxiter=500) |
4767 | +print "should be zero: F(p*) = ", F(ps) |
4768 | +print "p* =", ps |
4769 | +us=ul+F1l(ps) |
4770 | +print "u* =", us |
4771 | +print "should be equal to:", ur+F3r(ps) |
4772 | + |
4773 | +# the star region is divide by the contact discontinuity which gives a discontinuity in |
4774 | +# density rhosl/rhosr and internal energy iesl/iesr |
4775 | + |
4776 | +# left: |
4777 | +if ps<pl: # rarefaction |
4778 | + rhosl=rhol*(ps/pl)**(1/gamma) # eqn (3.1.172) |
4779 | +else: # shock |
4780 | + rhosl=rhol*( (gamma-1.)/(gamma+1.)*pl/ps +1. )/(pl/ps+(gamma-1.)/(gamma+1.)) # eqn (3.1.176) |
4781 | + s1=ul-al*sqrt((gamma+1.)/2./gamma*ps/pl + (gamma-1.)/2./gamma) # shock speed, eqn (3.1.177) |
4782 | +iesl=ps/rhosl/(gamma-1.) |
4783 | +asl=sqrt(gamma*ps/rhosl) |
4784 | + |
4785 | +# right: |
4786 | +if ps<pr: # rarefaction |
4787 | + rhosr=rhor*(ps/pr)**(1/gamma) # eqn (3.1.178) |
4788 | +else: # shock |
4789 | + rhosr=rhor*( ps/pr + (gamma-1.)/(gamma+1.) )/( (gamma-1.)/(gamma+1.)*ps/pr +1. ) # eqn (3.1.182) |
4790 | + s3=ur+ar*sqrt((gamma+1.)/2./gamma*ps/pr + (gamma-1.)/2./gamma) # shock speed, eqn (3.1.183) |
4791 | +iesr=ps/rhosr/(gamma-1.) |
4792 | +asr=sqrt(gamma*pr/rhosr) |
4793 | + |
4794 | +def solution(x,t): |
4795 | + if x/t<us: |
4796 | + # before the contact discontinuity: |
4797 | + |
4798 | + if ps<pl: |
4799 | + # u-a is a rarefaction wave |
4800 | + if x/t<ul-al: # left |
4801 | + return (pl, ul, rhol) |
4802 | + elif x/t<us-asl: # within the rarefaction wave |
4803 | + p=pl*(2./(gamma+1.) + (gamma-1.)/(gamma+1.)/al*(ul-x/t))**(2*gamma/(gamma-1.)) # eqn (3.1.97) |
4804 | + u=2./(gamma+1.)*(al + (gamma-1.)/2.*ul + x/t) # eqn (3.1.97) |
4805 | + rho=rhol*(2./(gamma+1.) + (gamma-1.)/(gamma+1.)/al*(ul-x/t))**(2./(gamma-1.)) # eqn (3.1.97) |
4806 | + return (p, u, rho) |
4807 | + else: # after the rarefaction before the contact disc. |
4808 | + return (ps,us,rhosl) |
4809 | + else: |
4810 | + # u-a is a shock wave |
4811 | + if x/t<s1: # left |
4812 | + return (pl, ul, rhol) |
4813 | + else: # between u-a shock and contact disc. |
4814 | + return (psl, us, rhosl) |
4815 | + |
4816 | + else: |
4817 | + # after the contact discontinuity: |
4818 | + |
4819 | + if ps<pr: |
4820 | + # u+a is a rarefaction wave |
4821 | + if x/t>ur+ar: # right |
4822 | + return (pr, ur, rhor) |
4823 | + elif x/t>usr+asr: # within the rarefaction wave |
4824 | + p=pr*(2./(gamma+1.) + (gamma-1.)/(gamma+1.)/ar*(ur-x/t))**(2*gamma/(gamma-1.)) # eqn (3.1.98) |
4825 | + u=2./(gamma+1.)*(-ar + (gamma-1.)/2.*ur + x/t) # eqn (3.1.98) |
4826 | + rho=rhor*(2./(gamma+1.) - (gamma-1.)/(gamma+1.)/ar*(ur-x/t))**(2./(gamma-1.)) # eqn (3.1.98) |
4827 | + return (p, u, rho) |
4828 | + else: # after the contact disc. before the rarefaction |
4829 | + return (ps,us, rhosr) |
4830 | + else: |
4831 | + # u+a is a shock wave |
4832 | + if x/t>s3: # right |
4833 | + return (pr, ur, rhor) |
4834 | + else: # between contact disc. and u-a shock |
4835 | + return (ps, us, rhosr) |
4836 | + |
4837 | |
4838 | === added file 'tests/mphase_dusty_gas_shock_tube/single_phase_frozen_flow_test.flml' |
4839 | --- tests/mphase_dusty_gas_shock_tube/single_phase_frozen_flow_test.flml 1970-01-01 00:00:00 +0000 |
4840 | +++ tests/mphase_dusty_gas_shock_tube/single_phase_frozen_flow_test.flml 2012-08-09 00:46:20 +0000 |
4841 | @@ -0,0 +1,379 @@ |
4842 | +<?xml version='1.0' encoding='utf-8'?> |
4843 | +<fluidity_options> |
4844 | + <simulation_name> |
4845 | + <string_value lines="1">single_phase_frozen_flow_test</string_value> |
4846 | + </simulation_name> |
4847 | + <problem_type> |
4848 | + <string_value lines="1">fluids</string_value> |
4849 | + </problem_type> |
4850 | + <geometry> |
4851 | + <dimension> |
4852 | + <integer_value rank="0">1</integer_value> |
4853 | + </dimension> |
4854 | + <mesh name="CoordinateMesh"> |
4855 | + <from_file file_name="line"> |
4856 | + <format name="triangle"/> |
4857 | + <stat> |
4858 | + <include_in_stat/> |
4859 | + </stat> |
4860 | + </from_file> |
4861 | + </mesh> |
4862 | + <mesh name="VelocityMesh"> |
4863 | + <from_mesh> |
4864 | + <mesh name="CoordinateMesh"/> |
4865 | + <mesh_shape> |
4866 | + <polynomial_degree> |
4867 | + <integer_value rank="0">1</integer_value> |
4868 | + </polynomial_degree> |
4869 | + </mesh_shape> |
4870 | + <stat> |
4871 | + <exclude_from_stat/> |
4872 | + </stat> |
4873 | + </from_mesh> |
4874 | + </mesh> |
4875 | + <mesh name="PressureMesh"> |
4876 | + <from_mesh> |
4877 | + <mesh name="CoordinateMesh"/> |
4878 | + <mesh_shape> |
4879 | + <polynomial_degree> |
4880 | + <integer_value rank="0">1</integer_value> |
4881 | + </polynomial_degree> |
4882 | + </mesh_shape> |
4883 | + <stat> |
4884 | + <exclude_from_stat/> |
4885 | + </stat> |
4886 | + </from_mesh> |
4887 | + </mesh> |
4888 | + <mesh name="DensityMesh"> |
4889 | + <from_mesh> |
4890 | + <mesh name="CoordinateMesh"/> |
4891 | + <mesh_shape> |
4892 | + <polynomial_degree> |
4893 | + <integer_value rank="0">1</integer_value> |
4894 | + </polynomial_degree> |
4895 | + </mesh_shape> |
4896 | + <stat> |
4897 | + <exclude_from_stat/> |
4898 | + </stat> |
4899 | + </from_mesh> |
4900 | + </mesh> |
4901 | + <quadrature> |
4902 | + <degree> |
4903 | + <integer_value rank="0">5</integer_value> |
4904 | + </degree> |
4905 | + </quadrature> |
4906 | + </geometry> |
4907 | + <io> |
4908 | + <dump_format> |
4909 | + <string_value>vtk</string_value> |
4910 | + </dump_format> |
4911 | + <dump_period> |
4912 | + <constant> |
4913 | + <real_value rank="0">0.0</real_value> |
4914 | + </constant> |
4915 | + </dump_period> |
4916 | + <output_mesh name="VelocityMesh"/> |
4917 | + <stat> |
4918 | + <output_at_start/> |
4919 | + </stat> |
4920 | + </io> |
4921 | + <timestepping> |
4922 | + <current_time> |
4923 | + <real_value rank="0">0.0</real_value> |
4924 | + </current_time> |
4925 | + <timestep> |
4926 | + <real_value rank="0">0.000001</real_value> |
4927 | + </timestep> |
4928 | + <finish_time> |
4929 | + <real_value rank="0">0.00037772998222</real_value> |
4930 | + <comment>This corresponds to a non-dimensional time of tau = 4. Multiplying this by the reference time of l/u_ref gives 0.00037772998222 seconds. |
4931 | + |
4932 | +l = (4/3)*(rho_p/rho_g)*d |
4933 | +u_ref = sqrt(p/rho_g)</comment> |
4934 | + </finish_time> |
4935 | + <nonlinear_iterations> |
4936 | + <integer_value rank="0">3</integer_value> |
4937 | + <tolerance> |
4938 | + <real_value rank="0">1.0e-9</real_value> |
4939 | + <infinity_norm/> |
4940 | + </tolerance> |
4941 | + </nonlinear_iterations> |
4942 | + </timestepping> |
4943 | + <material_phase name="Air"> |
4944 | + <equation_of_state> |
4945 | + <compressible> |
4946 | + <stiffened_gas> |
4947 | + <ratio_specific_heats> |
4948 | + <real_value rank="0">1.4</real_value> |
4949 | + </ratio_specific_heats> |
4950 | + </stiffened_gas> |
4951 | + </compressible> |
4952 | + </equation_of_state> |
4953 | + <scalar_field name="Pressure" rank="0"> |
4954 | + <prognostic> |
4955 | + <mesh name="PressureMesh"/> |
4956 | + <spatial_discretisation> |
4957 | + <continuous_galerkin> |
4958 | + <remove_stabilisation_term/> |
4959 | + </continuous_galerkin> |
4960 | + </spatial_discretisation> |
4961 | + <scheme> |
4962 | + <poisson_pressure_solution> |
4963 | + <string_value lines="1">never</string_value> |
4964 | + </poisson_pressure_solution> |
4965 | + <use_compressible_projection_method/> |
4966 | + </scheme> |
4967 | + <solver> |
4968 | + <iterative_method name="preonly"/> |
4969 | + <preconditioner name="lu"/> |
4970 | + <relative_error> |
4971 | + <real_value rank="0">1e-7</real_value> |
4972 | + </relative_error> |
4973 | + <max_iterations> |
4974 | + <integer_value rank="0">10000</integer_value> |
4975 | + </max_iterations> |
4976 | + <never_ignore_solver_failures/> |
4977 | + <diagnostics> |
4978 | + <monitors/> |
4979 | + </diagnostics> |
4980 | + </solver> |
4981 | + <output/> |
4982 | + <stat/> |
4983 | + <convergence> |
4984 | + <include_in_convergence/> |
4985 | + </convergence> |
4986 | + <detectors> |
4987 | + <exclude_from_detectors/> |
4988 | + </detectors> |
4989 | + <steady_state> |
4990 | + <include_in_steady_state/> |
4991 | + </steady_state> |
4992 | + <no_interpolation/> |
4993 | + </prognostic> |
4994 | + </scalar_field> |
4995 | + <scalar_field name="Density" rank="0"> |
4996 | + <prognostic> |
4997 | + <mesh name="PressureMesh"/> |
4998 | + <spatial_discretisation> |
4999 | + <continuous_galerkin> |
5000 | + <stabilisation> |
Hey there. I've only looked at Advection_ Diffusion_ CG.F90 yet. Looking great so far, the maths of the changes in the i.e. equation seem to make sense to me. Here it goes:
* Currently the multiphase logical is only set within the type==FIELD_ EQUATION_ INTERNALENERGY case. Can we either: internal_ energy_ equation
equation_
- always set multiphase to be true or false, or
- rename it to mean multiphase_
As it is now suggested that multiphase can be used independently of
the value of equation_type.
* Just being picky here: can we get rid of the nullify(vfrac) in line 506. It seems
to suggest that vfrac is used later on, which it's not.
* Advection_ Diffusion_ CG.F90, line 517:
Why does the behaviour change if we're not multiphase? Do you maybe mean:
if (.not. multiphase .or. have_option( ....compressibl e) then
* Advection_ Diffusion_ CG.F90, line 794 having an allocate within an element
loop is not ideal. It's not efficient, and it's also not thread-safe. I can
see the problem with following the usual approach (nvfrac might
not be defined), but I think it's better to do some kludge with allocating
a dummy nvfrac or something rather.
* Advection_ Diffusion_ CG.F90, lines 1492 and 1502: I don't think you intended to
add (istate) there.
* Not too excited about having to pass in an array of states diffusion_ cg now tbh. One way to keep the damage limited would be equation_ cg. And keep everything else just limited to a
to advection_
to move the call to add_heat_transfer up one level inside
solve_field_
single state, i.e. not passing istate all the way down and changing state ->
state(istate) *everywhere*.