Merge lp:~fluidity-core/fluidity/shallow-water-dev into lp:fluidity
- shallow-water-dev
- Merge into dev-trunk
Status: | Work in progress |
---|---|
Proposed branch: | lp:~fluidity-core/fluidity/shallow-water-dev |
Merge into: | lp:fluidity |
Diff against target: |
68993 lines (has conflicts)
Text conflict in femtools/Makefile.dependencies Text conflict in femtools/Makefile.in Text conflict in femtools/doc/femtools_manual.pdf Text conflict in main/Shallow_Water.F90 Text conflict in tests/Helmholtz-anisotropic-vector-mms-p1p1cg/MMS_A.flml Text conflict in tests/Helmholtz-anisotropic-vector-mms-p1p1cg/MMS_B.flml Text conflict in tests/Helmholtz-anisotropic-vector-mms-p1p1cg/MMS_C.flml Text conflict in tests/Helmholtz-anisotropic-vector-mms-p1p1cg/MMS_D.flml Text conflict in tests/Helmholtz-anisotropic-vector-mms-p1p1cg/MMS_E.flml Text conflict in tests/data/cg_interpolation_A.flml Text conflict in tests/optimality/dummy.swml Contents conflict in tests/shallow_water_adjoint_default_controls_2d/adjoint_B.swml Contents conflict in tests/shallow_water_adjoint_default_controls_2d/adjoint_C.swml Contents conflict in tests/shallow_water_adjoint_default_controls_2d/adjoint_D.swml Contents conflict in tests/shallow_water_adjoint_default_controls_2d/adjoint_E.swml Text conflict in tests/shallow_water_adjoint_func_eval_2d/adjoint_template.swml Contents conflict in tests/shallow_water_optimisation/adjoint_B.swml Contents conflict in tests/shallow_water_optimisation/adjoint_C.swml Contents conflict in tests/shallow_water_optimisation/adjoint_D.swml Contents conflict in tests/shallow_water_optimisation/adjoint_E.swml Text conflict in tests/shallow_water_optimisation/adjoint_template.swml Text conflict in tests/shallow_water_optimisation_check_gradient_2d/adjoint_A.swml |
To merge this branch: | bzr merge lp:~fluidity-core/fluidity/shallow-water-dev |
Related bugs: |
Reviewer | Review Type | Date Requested | Status |
---|---|---|---|
Cian Wilson | Needs Information | ||
David Ham | Pending | ||
Review via email: mp+81751@code.launchpad.net |
Commit message
Description of the change
This shallow-water code includes the changes necessary to correctly interpolate vector fields from the sphere for BDFM1 elements, plus a suite of test cases.
Dr Colin J Cotter (colin-cotter) wrote : | # |
Hi Cian,
These options are read in populate_state and make things happen in
Elements.F90 and various other places.
--cjc
On 09/11/11 16:17, Cian Wilson wrote:
> Review: Needs Information
>
> Hi Colin,
>
> Only had a quick glance but can I check if these changes provide any new functionality to fluidity itself? I ask because you've changed mesh_options.rn[cg] and these options are now appearing in the fluidity schema too but I can't see how they're applicable to fluidity as I can only trace them back to the hybridized helmholtz solver. If they're not used in fluidity then I don't think they should appear in a piece of the schema that is visible within an flml but I may be missing something.
>
> Unrelated to this merge but the parent option (geometry/
>
> Cheers,
> Cian
- 3579. By Dr Colin J Cotter
-
A few fixes for adjoint, plus tweaking some swmls.
Cian Wilson (cwilson) wrote : | # |
Hi Colin,
Could you give me line numbers in fluidity code where this happens please?
When I open up diamond to set up a new flml I get option of setting:
/geometry/
and
/geometry/
However, when I grep for where these options are used I get:
cwilson@
assemble/
schemas/
schemas/
tests/bdfm1_
None of which appear to be fluidity related to me so shouldn't be available in an flml.
Cheers,
Cian
> Hi Cian,
> These options are read in populate_state and make things happen in
> Elements.F90 and various other places.
>
> --cjc
>
> On 09/11/11 16:17, Cian Wilson wrote:
> > Review: Needs Information
> >
> > Hi Colin,
> >
> > Only had a quick glance but can I check if these changes provide any new
> functionality to fluidity itself? I ask because you've changed
> mesh_options.rn[cg] and these options are now appearing in the fluidity schema
> too but I can't see how they're applicable to fluidity as I can only trace
> them back to the hybridized helmholtz solver. If they're not used in fluidity
> then I don't think they should appear in a piece of the schema that is visible
> within an flml but I may be missing something.
> >
> > Unrelated to this merge but the parent option
> (geometry/
> in the schema.
> >
> > Cheers,
> > Cian
Cian Wilson (cwilson) wrote : | # |
I can see where the already existing parent option /geometry/
The parent option does need some documentation in the schema however.
- 3580. By Patrick Farrell
-
Some XML fixes for Colin.
- 3581. By Patrick Farrell
-
Let's not put in a NonlinearVelocity in the adjoint run, for now.
- 3582. By Dr Colin J Cotter
-
Change place where option theta is looked for.
- 3583. By Dr Colin J Cotter
-
theta set.
- 3584. By Dr Colin J Cotter
-
Removing some dead tests and fixing some swmls.
Dr Colin J Cotter (colin-cotter) wrote : | # |
Oh I see, fair point. I'll fix that.
all the best
--cjc
On 09/11/11 19:01, Cian Wilson wrote:
> I can see where the already existing parent option /geometry/
>
> The parent option does need some documentation in the schema however.
- 3585. By Dr Colin J Cotter
-
Move continuity debugging options to shallow_
water_options. rnc as
requested by Cian. - 3586. By Dr Colin J Cotter
-
Merge from trunk.
- 3587. By Dr Colin J Cotter
-
Merged from the trunk.
- 3588. By Dr Colin J Cotter
-
Forgot that this file has been renamed.
David Ham (david-ham) wrote : | # |
I think this looks OK. Make sure it is current with trunk and check that all the major functionality is tested.
- 3589. By Dr Colin J Cotter
-
Merged from trunk.
- 3590. By Dr Colin J Cotter
-
Uncommented code that I previously forgot to uncomment.
- 3591. By Dr Colin J Cotter
-
Get theta from correct place in adjoint callbacks.
- 3592. By Dr Colin J Cotter
-
Get theta from correct place in adjoint control callbacks too.
- 3593. By Dr Colin J Cotter
-
Some dead tests.
- 3594. By Dr Colin J Cotter
-
Remove pointer status from scalar field.
- 3595. By Dr Colin J Cotter
-
Don't output to the statfile twice at the end of the simulation.
This fixes shallow_
water_adjoint_ func_eval. - 3596. By Dr Colin J Cotter
-
Merge from trunk.
- 3597. By Dr Colin J Cotter
-
Increased timesteps to 2 in test. This is not the solution, but I
can't work out why it doesn't work otherwise. - 3598. By Dr Colin J Cotter
-
Need to calculate functional values before writing diagnostics.
- 3599. By Dr Colin J Cotter
-
Made makefiles and got a merge from the trunk.
- 3600. By Dr Colin J Cotter
-
Merge from trunk.
- 3601. By Dr Colin J Cotter
-
Error in merging makefiles.
- 3602. By Jemma Shipton
-
merging changes from Colin's bdfm1-nonlinear-sw branch
- 3603. By Jemma Shipton
-
get global numbering fix for quads from trunk
Unmerged revisions
- 3603. By Jemma Shipton
-
get global numbering fix for quads from trunk
- 3602. By Jemma Shipton
-
merging changes from Colin's bdfm1-nonlinear-sw branch
- 3601. By Dr Colin J Cotter
-
Error in merging makefiles.
- 3600. By Dr Colin J Cotter
-
Merge from trunk.
- 3599. By Dr Colin J Cotter
-
Made makefiles and got a merge from the trunk.
- 3598. By Dr Colin J Cotter
-
Need to calculate functional values before writing diagnostics.
- 3597. By Dr Colin J Cotter
-
Increased timesteps to 2 in test. This is not the solution, but I
can't work out why it doesn't work otherwise. - 3596. By Dr Colin J Cotter
-
Merge from trunk.
- 3595. By Dr Colin J Cotter
-
Don't output to the statfile twice at the end of the simulation.
This fixes shallow_
water_adjoint_ func_eval. - 3594. By Dr Colin J Cotter
-
Remove pointer status from scalar field.
Preview Diff
1 | === modified file 'Makefile.in' | |||
2 | --- Makefile.in 2012-09-19 11:09:17 +0000 | |||
3 | +++ Makefile.in 2012-10-23 15:04:25 +0000 | |||
4 | @@ -186,7 +186,9 @@ | |||
5 | 186 | @echo " LD $(PROGRAM)" | 186 | @echo " LD $(PROGRAM)" |
6 | 187 | @$(EVAL) $(LINKER) -o $(PROGRAM) main.o $(LIBS) | 187 | @$(EVAL) $(LINKER) -o $(PROGRAM) main.o $(LIBS) |
7 | 188 | 188 | ||
9 | 189 | sw: bin/shallow_water | 189 | sw: bin/shallow_water |
10 | 190 | swm: bin/shallow_water_mimetic | ||
11 | 191 | test_sw: bin/test_sw_advection | ||
12 | 190 | be: bin/burgers_equation | 192 | be: bin/burgers_equation |
13 | 191 | hh: bin/hybridized_helmholtz_solver | 193 | hh: bin/hybridized_helmholtz_solver |
14 | 192 | 194 | ||
15 | @@ -199,6 +201,22 @@ | |||
16 | 199 | @$(EVAL) $(FLLINKER) -o bin/hybridized_helmholtz_solver main/Hybridized_Helmholtz_Solver.o $(LIBS) | 201 | @$(EVAL) $(FLLINKER) -o bin/hybridized_helmholtz_solver main/Hybridized_Helmholtz_Solver.o $(LIBS) |
17 | 200 | @$(EVAL) $(FLLINKER) -o bin/shallow_water main/Shallow_Water.o $(LIBS) | 202 | @$(EVAL) $(FLLINKER) -o bin/shallow_water main/Shallow_Water.o $(LIBS) |
18 | 201 | 203 | ||
19 | 204 | bin/shallow_water_mimetic: fluidity_library main/Shallow_Water_Mimetic.F90 | ||
20 | 205 | @cd main; $(MAKE) Shallow_Water_Mimetic.o | ||
21 | 206 | @echo "BUILD shallow_water_mimetic" | ||
22 | 207 | @echo " MKDIR bin" | ||
23 | 208 | @mkdir -p bin | ||
24 | 209 | @echo " LD shallow_water_mimetic" | ||
25 | 210 | @$(EVAL) $(FLLINKER) -o bin/shallow_water_mimetic main/Shallow_Water_Mimetic.o $(LIBS) | ||
26 | 211 | |||
27 | 212 | bin/test_sw_advection: fluidity_library main/test_sw_advection.F90 | ||
28 | 213 | @cd main; $(MAKE) test_sw_advection.o | ||
29 | 214 | @echo "BUILD test_sw_advection" | ||
30 | 215 | @echo " MKDIR bin" | ||
31 | 216 | @mkdir -p bin | ||
32 | 217 | @echo " LD test_sw_advection" | ||
33 | 218 | @$(EVAL) $(LINKER) -o bin/test_sw_advection main/test_sw_advection.o $(LIBS) | ||
34 | 219 | |||
35 | 202 | bin/hybridized_helmholtz_solver: fluidity_library main/Hybridized_Helmholtz_Solver.F90 | 220 | bin/hybridized_helmholtz_solver: fluidity_library main/Hybridized_Helmholtz_Solver.F90 |
36 | 203 | @cd main; $(MAKE) Hybridized_Helmholtz_Solver.o | 221 | @cd main; $(MAKE) Hybridized_Helmholtz_Solver.o |
37 | 204 | @echo "BUILD hybridized_helmholtz_solver" | 222 | @echo "BUILD hybridized_helmholtz_solver" |
38 | 205 | 223 | ||
39 | === modified file 'adjoint/Burgers_Control_Callbacks.F90' | |||
40 | --- adjoint/Burgers_Control_Callbacks.F90 2011-07-07 16:30:31 +0000 | |||
41 | +++ adjoint/Burgers_Control_Callbacks.F90 2012-10-23 15:04:25 +0000 | |||
42 | @@ -34,7 +34,7 @@ | |||
43 | 34 | use state_module | 34 | use state_module |
44 | 35 | use python_state | 35 | use python_state |
45 | 36 | use sparse_matrices_fields | 36 | use sparse_matrices_fields |
47 | 37 | use manifold_projections | 37 | use manifold_tools |
48 | 38 | use mangle_dirichlet_rows_module | 38 | use mangle_dirichlet_rows_module |
49 | 39 | 39 | ||
50 | 40 | implicit none | 40 | implicit none |
51 | 41 | 41 | ||
52 | === modified file 'adjoint/Makefile.dependencies' | |||
53 | --- adjoint/Makefile.dependencies 2011-07-19 16:55:03 +0000 | |||
54 | +++ adjoint/Makefile.dependencies 2012-10-23 15:04:25 +0000 | |||
55 | @@ -70,10 +70,9 @@ | |||
56 | 70 | Burgers_Control_Callbacks.o ../include/burgers_adjoint_controls.mod: \ | 70 | Burgers_Control_Callbacks.o ../include/burgers_adjoint_controls.mod: \ |
57 | 71 | Burgers_Control_Callbacks.F90 ../include/fdebug.h ../include/fields.mod \ | 71 | Burgers_Control_Callbacks.F90 ../include/fdebug.h ../include/fields.mod \ |
58 | 72 | ../include/global_parameters.mod \ | 72 | ../include/global_parameters.mod \ |
63 | 73 | ../include/mangle_dirichlet_rows_module.mod \ | 73 | ../include/mangle_dirichlet_rows_module.mod ../include/manifold_tools.mod \ |
64 | 74 | ../include/manifold_projections.mod ../include/python_state.mod \ | 74 | ../include/python_state.mod ../include/sparse_matrices_fields.mod \ |
65 | 75 | ../include/sparse_matrices_fields.mod ../include/spud.mod \ | 75 | ../include/spud.mod ../include/state_module.mod |
62 | 76 | ../include/state_module.mod | ||
66 | 77 | 76 | ||
67 | 78 | ../include/forward_main_loop.mod: Forward_Main_Loop.o | 77 | ../include/forward_main_loop.mod: Forward_Main_Loop.o |
68 | 79 | @true | 78 | @true |
69 | @@ -128,7 +127,7 @@ | |||
70 | 128 | Shallow_Water_Adjoint_Callbacks.F90 ../include/adjoint_global_variables.mod \ | 127 | Shallow_Water_Adjoint_Callbacks.F90 ../include/adjoint_global_variables.mod \ |
71 | 129 | ../include/fdebug.h ../include/fields.mod ../include/global_parameters.mod \ | 128 | ../include/fdebug.h ../include/fields.mod ../include/global_parameters.mod \ |
72 | 130 | ../include/libadjoint_data_callbacks.mod ../include/mangle_options_tree.mod \ | 129 | ../include/libadjoint_data_callbacks.mod ../include/mangle_options_tree.mod \ |
74 | 131 | ../include/manifold_projections.mod ../include/populate_state_module.mod \ | 130 | ../include/manifold_tools.mod ../include/populate_state_module.mod \ |
75 | 132 | ../include/sparse_matrices_fields.mod ../include/spud.mod \ | 131 | ../include/sparse_matrices_fields.mod ../include/spud.mod \ |
76 | 133 | ../include/state_module.mod | 132 | ../include/state_module.mod |
77 | 134 | 133 | ||
78 | @@ -140,7 +139,7 @@ | |||
79 | 140 | ../include/shallow_water_adjoint_controls.mod: \ | 139 | ../include/shallow_water_adjoint_controls.mod: \ |
80 | 141 | Shallow_Water_Control_Callbacks.F90 ../include/fdebug.h \ | 140 | Shallow_Water_Control_Callbacks.F90 ../include/fdebug.h \ |
81 | 142 | ../include/fields.mod ../include/global_parameters.mod \ | 141 | ../include/fields.mod ../include/global_parameters.mod \ |
83 | 143 | ../include/manifold_projections.mod ../include/python_state.mod \ | 142 | ../include/manifold_tools.mod ../include/python_state.mod \ |
84 | 144 | ../include/sparse_matrices_fields.mod ../include/spud.mod \ | 143 | ../include/sparse_matrices_fields.mod ../include/spud.mod \ |
85 | 145 | ../include/state_module.mod | 144 | ../include/state_module.mod |
86 | 146 | 145 | ||
87 | 147 | 146 | ||
88 | === modified file 'adjoint/Shallow_Water_Adjoint_Callbacks.F90' | |||
89 | --- adjoint/Shallow_Water_Adjoint_Callbacks.F90 2011-06-09 09:43:33 +0000 | |||
90 | +++ adjoint/Shallow_Water_Adjoint_Callbacks.F90 2012-10-23 15:04:25 +0000 | |||
91 | @@ -39,7 +39,7 @@ | |||
92 | 39 | use spud, only: get_option | 39 | use spud, only: get_option |
93 | 40 | use adjoint_global_variables, only: adj_path_lookup | 40 | use adjoint_global_variables, only: adj_path_lookup |
94 | 41 | use mangle_options_tree, only: adjoint_field_path | 41 | use mangle_options_tree, only: adjoint_field_path |
96 | 42 | use manifold_projections | 42 | use manifold_tools |
97 | 43 | use global_parameters, only: running_adjoint | 43 | use global_parameters, only: running_adjoint |
98 | 44 | use populate_state_module | 44 | use populate_state_module |
99 | 45 | implicit none | 45 | implicit none |
100 | @@ -488,7 +488,7 @@ | |||
101 | 488 | path = adjoint_field_path(path) | 488 | path = adjoint_field_path(path) |
102 | 489 | end if | 489 | end if |
103 | 490 | 490 | ||
105 | 491 | call get_option(trim(path) // "/prognostic/temporal_discretisation/theta", theta) | 491 | call get_option("/timestepping/theta", theta) |
106 | 492 | call get_option(trim(path) // "/prognostic/mean_layer_thickness", d0) | 492 | call get_option(trim(path) // "/prognostic/mean_layer_thickness", d0) |
107 | 493 | 493 | ||
108 | 494 | if (hermitian==ADJ_FALSE) then | 494 | if (hermitian==ADJ_FALSE) then |
109 | @@ -858,7 +858,7 @@ | |||
110 | 858 | 858 | ||
111 | 859 | ierr = adj_dict_find(adj_path_lookup, "Fluid::LayerThickness", path) | 859 | ierr = adj_dict_find(adj_path_lookup, "Fluid::LayerThickness", path) |
112 | 860 | call adj_chkierr(ierr) | 860 | call adj_chkierr(ierr) |
114 | 861 | call get_option(trim(path) // "/prognostic/temporal_discretisation/theta", theta) | 861 | call get_option("/timestepping/theta",theta) |
115 | 862 | call get_option(trim(path) // "/prognostic/mean_layer_thickness", d0) | 862 | call get_option(trim(path) // "/prognostic/mean_layer_thickness", d0) |
116 | 863 | 863 | ||
117 | 864 | eta_mesh => extract_mesh(matrices, "LayerThicknessMesh") | 864 | eta_mesh => extract_mesh(matrices, "LayerThicknessMesh") |
118 | 865 | 865 | ||
119 | === modified file 'adjoint/Shallow_Water_Control_Callbacks.F90' | |||
120 | --- adjoint/Shallow_Water_Control_Callbacks.F90 2011-06-11 13:26:55 +0000 | |||
121 | +++ adjoint/Shallow_Water_Control_Callbacks.F90 2012-10-23 15:04:25 +0000 | |||
122 | @@ -34,7 +34,7 @@ | |||
123 | 34 | use state_module | 34 | use state_module |
124 | 35 | use python_state | 35 | use python_state |
125 | 36 | use sparse_matrices_fields | 36 | use sparse_matrices_fields |
127 | 37 | use manifold_projections | 37 | use manifold_tools |
128 | 38 | 38 | ||
129 | 39 | implicit none | 39 | implicit none |
130 | 40 | 40 | ||
131 | @@ -154,7 +154,7 @@ | |||
132 | 154 | adj_vfield => extract_vector_field(states(state_id), "AdjointLocalVelocityDelta") | 154 | adj_vfield => extract_vector_field(states(state_id), "AdjointLocalVelocityDelta") |
133 | 155 | ! We need the AdjointVelocity for the option path only | 155 | ! We need the AdjointVelocity for the option path only |
134 | 156 | adj_sfield2 => extract_scalar_field(states(state_id), "AdjointLayerThickness") | 156 | adj_sfield2 => extract_scalar_field(states(state_id), "AdjointLayerThickness") |
136 | 157 | call get_option(trim(adj_sfield2%option_path) // "/prognostic/temporal_discretisation/theta", theta) | 157 | call get_option("/timestepping/theta", theta) |
137 | 158 | call get_option(trim(adj_sfield2%option_path) // "/prognostic/mean_layer_thickness", d0) | 158 | call get_option(trim(adj_sfield2%option_path) // "/prognostic/mean_layer_thickness", d0) |
138 | 159 | 159 | ||
139 | 160 | big_mat => extract_block_csr_matrix(matrices, "InverseBigMatrix") | 160 | big_mat => extract_block_csr_matrix(matrices, "InverseBigMatrix") |
140 | 161 | 161 | ||
141 | === modified file 'assemble/Advection_Diffusion_DG.F90' | |||
142 | --- assemble/Advection_Diffusion_DG.F90 2012-10-11 17:52:01 +0000 | |||
143 | +++ assemble/Advection_Diffusion_DG.F90 2012-10-23 15:04:25 +0000 | |||
144 | @@ -1547,7 +1547,7 @@ | |||
145 | 1547 | end if | 1547 | end if |
146 | 1548 | end if | 1548 | end if |
147 | 1549 | end if | 1549 | end if |
149 | 1550 | 1550 | ||
150 | 1551 | ! Add stabilisation to the advection term if requested by the user. | 1551 | ! Add stabilisation to the advection term if requested by the user. |
151 | 1552 | if (stabilisation_scheme==UPWIND) then | 1552 | if (stabilisation_scheme==UPWIND) then |
152 | 1553 | ! NOTE: U_nl_q may (or may not) have been modified by the grid velocity | 1553 | ! NOTE: U_nl_q may (or may not) have been modified by the grid velocity |
153 | @@ -2603,7 +2603,6 @@ | |||
154 | 2603 | nnAdvection_in=shape_shape(T_shape, T_shape_2, & | 2603 | nnAdvection_in=shape_shape(T_shape, T_shape_2, & |
155 | 2604 | & outer_advection_integral * detwei) | 2604 | & outer_advection_integral * detwei) |
156 | 2605 | 2605 | ||
157 | 2606 | |||
158 | 2607 | end if | 2606 | end if |
159 | 2608 | if (include_diffusion.or.have_buoyancy_adjustment_by_vertical_diffusion) then | 2607 | if (include_diffusion.or.have_buoyancy_adjustment_by_vertical_diffusion) then |
160 | 2609 | 2608 | ||
161 | 2610 | 2609 | ||
162 | === added file 'assemble/Advection_Local_DG.F90' | |||
163 | --- assemble/Advection_Local_DG.F90 1970-01-01 00:00:00 +0000 | |||
164 | +++ assemble/Advection_Local_DG.F90 2012-10-23 15:04:25 +0000 | |||
165 | @@ -0,0 +1,2538 @@ | |||
166 | 1 | ! Copyright (C) 2006 Imperial College London and others. | ||
167 | 2 | ! | ||
168 | 3 | ! Please see the AUTHORS file in the main source directory for a full list | ||
169 | 4 | ! of copyright holders. | ||
170 | 5 | ! | ||
171 | 6 | ! Prof. C Pain | ||
172 | 7 | ! Applied Modelling and Computation Group | ||
173 | 8 | ! Department of Earth Science and Engineering | ||
174 | 9 | ! Imperial College London | ||
175 | 10 | ! | ||
176 | 11 | ! amcgsoftware@imperial.ac.uk | ||
177 | 12 | ! | ||
178 | 13 | ! This library is free software; you can redistribute it and/or | ||
179 | 14 | ! modify it under the terms of the GNU Lesser General Public | ||
180 | 15 | ! License as published by the Free Software Foundation, | ||
181 | 16 | ! version 2.1 of the License. | ||
182 | 17 | ! | ||
183 | 18 | ! This library is distributed in the hope that it will be useful, | ||
184 | 19 | ! but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
185 | 20 | ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | ||
186 | 21 | ! Lesser General Public License for more details. | ||
187 | 22 | ! | ||
188 | 23 | ! You should have received a copy of the GNU Lesser General Public | ||
189 | 24 | ! License along with this library; if not, write to the Free Software | ||
190 | 25 | ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 | ||
191 | 26 | ! USA | ||
192 | 27 | #include "fdebug.h" | ||
193 | 28 | |||
194 | 29 | module advection_local_DG | ||
195 | 30 | !!< This module contains the Discontinuous Galerkin form of the advection | ||
196 | 31 | !!< equation when vector fields are represented in the local coordinates | ||
197 | 32 | !!< from the Piola transformation. | ||
198 | 33 | use elements | ||
199 | 34 | use global_parameters, only:current_debug_level, OPTION_PATH_LEN | ||
200 | 35 | use sparse_tools | ||
201 | 36 | use fetools | ||
202 | 37 | use dgtools | ||
203 | 38 | use fields | ||
204 | 39 | use fefields | ||
205 | 40 | use state_module | ||
206 | 41 | use shape_functions | ||
207 | 42 | use global_numbering | ||
208 | 43 | use transform_elements | ||
209 | 44 | use vector_tools | ||
210 | 45 | use fldebug | ||
211 | 46 | use vtk_interfaces | ||
212 | 47 | use Coordinates | ||
213 | 48 | use Tensors | ||
214 | 49 | use petsc_solve_state_module | ||
215 | 50 | use spud | ||
216 | 51 | use slope_limiters_dg | ||
217 | 52 | use sparsity_patterns | ||
218 | 53 | use sparse_matrices_fields | ||
219 | 54 | use sparsity_patterns_meshes | ||
220 | 55 | use Vector_Tools | ||
221 | 56 | use manifold_tools | ||
222 | 57 | use bubble_tools | ||
223 | 58 | use diagnostic_fields, only: calculate_diagnostic_variable | ||
224 | 59 | use global_parameters, only : FIELD_NAME_LEN | ||
225 | 60 | |||
226 | 61 | implicit none | ||
227 | 62 | |||
228 | 63 | ! Buffer for output messages. | ||
229 | 64 | character(len=255), private :: message | ||
230 | 65 | |||
231 | 66 | private | ||
232 | 67 | public solve_advection_dg_subcycle, solve_vector_advection_dg_subcycle,& | ||
233 | 68 | & solve_advection_cg_tracer, compute_U_residual, get_pv | ||
234 | 69 | |||
235 | 70 | !parameters specifying vector upwind options | ||
236 | 71 | integer, parameter :: VECTOR_UPWIND_EDGE=1, VECTOR_UPWIND_SPHERE=2 | ||
237 | 72 | |||
238 | 73 | ! Local private control parameters. These are module-global parameters | ||
239 | 74 | ! because it would be expensive and/or inconvenient to re-evaluate them | ||
240 | 75 | ! on a per-element or per-face basis | ||
241 | 76 | real :: dt | ||
242 | 77 | |||
243 | 78 | ! Whether to include various terms | ||
244 | 79 | logical :: include_advection | ||
245 | 80 | contains | ||
246 | 81 | |||
247 | 82 | subroutine solve_advection_dg_subcycle(field_name, state, velocity_name,& | ||
248 | 83 | &continuity, flux) | ||
249 | 84 | !!< Construct and solve the advection equation for the given | ||
250 | 85 | !!< field using discontinuous elements. | ||
251 | 86 | |||
252 | 87 | !! Name of the field to be solved for. | ||
253 | 88 | character(len=*), intent(in) :: field_name | ||
254 | 89 | !! Collection of fields defining system state. | ||
255 | 90 | type(state_type), intent(inout) :: state | ||
256 | 91 | !! Optional velocity name | ||
257 | 92 | character(len = *), intent(in) :: velocity_name | ||
258 | 93 | !! Solve continuity equation as opposed to advection equation | ||
259 | 94 | logical, intent(in), optional :: continuity | ||
260 | 95 | !! Return a flux variable such that T-T_old=div(Flux) | ||
261 | 96 | type(vector_field), intent(inout), optional :: Flux | ||
262 | 97 | |||
263 | 98 | !! Tracer to be solved for. | ||
264 | 99 | type(scalar_field), pointer :: T, T_old, s_field | ||
265 | 100 | |||
266 | 101 | !! Coordinate and velocity fields | ||
267 | 102 | type(vector_field), pointer :: X, U_nl | ||
268 | 103 | |||
269 | 104 | !! Velocity field in Cartesian coordinates for slope limiter | ||
270 | 105 | type(vector_field) :: U_nl_cartesian | ||
271 | 106 | |||
272 | 107 | !! Change in T over one timestep. | ||
273 | 108 | type(scalar_field) :: delta_T, delta_T_total | ||
274 | 109 | |||
275 | 110 | !! Sparsity of advection matrix. | ||
276 | 111 | type(csr_sparsity), pointer :: sparsity | ||
277 | 112 | |||
278 | 113 | !! System matrix. | ||
279 | 114 | type(csr_matrix) :: matrix, mass, inv_mass | ||
280 | 115 | |||
281 | 116 | !! Sparsity of mass matrix. | ||
282 | 117 | type(csr_sparsity) :: mass_sparsity | ||
283 | 118 | |||
284 | 119 | !! Right hand side vector. | ||
285 | 120 | type(scalar_field) :: rhs | ||
286 | 121 | |||
287 | 122 | !! Whether to invoke the slope limiter | ||
288 | 123 | logical :: limit_slope | ||
289 | 124 | !! Which limiter to use | ||
290 | 125 | integer :: limiter | ||
291 | 126 | |||
292 | 127 | !! Number of advection subcycles. | ||
293 | 128 | integer :: subcycles | ||
294 | 129 | real :: max_courant_number | ||
295 | 130 | |||
296 | 131 | !! Flux reconstruction stuff: | ||
297 | 132 | !! Upwind fluxes trace variable | ||
298 | 133 | type(scalar_field) :: UpwindFlux | ||
299 | 134 | !! Mesh where UpwindFlux lives | ||
300 | 135 | type(mesh_type), pointer :: lambda_mesh | ||
301 | 136 | !! Upwind Flux matrix | ||
302 | 137 | type(csr_matrix) :: UpwindFluxMatrix | ||
303 | 138 | |||
304 | 139 | character(len=FIELD_NAME_LEN) :: limiter_name | ||
305 | 140 | integer :: i, ele | ||
306 | 141 | real :: meancheck | ||
307 | 142 | ewrite(1,*) 'subroutine solve_advection_dg_subcycle' | ||
308 | 143 | |||
309 | 144 | |||
310 | 145 | T=>extract_scalar_field(state, field_name) | ||
311 | 146 | T_old=>extract_scalar_field(state, "Old"//field_name) | ||
312 | 147 | X=>extract_vector_field(state, "Coordinate") | ||
313 | 148 | U_nl=>extract_vector_field(state, velocity_name) | ||
314 | 149 | |||
315 | 150 | ewrite(2,*) 'UVALS in dg', maxval(abs(U_nl%val)) | ||
316 | 151 | |||
317 | 152 | ! Reset T to value at the beginning of the timestep. | ||
318 | 153 | call set(T, T_old) | ||
319 | 154 | if(present(Flux)) then | ||
320 | 155 | call zero(Flux) | ||
321 | 156 | end if | ||
322 | 157 | |||
323 | 158 | sparsity => get_csr_sparsity_firstorder(state, T%mesh, T%mesh) | ||
324 | 159 | call allocate(matrix, sparsity) ! Add data space to the sparsity | ||
325 | 160 | ! pattern. | ||
326 | 161 | call zero(matrix) | ||
327 | 162 | |||
328 | 163 | !Flux reconstruction allocations | ||
329 | 164 | if(present(Flux)) then | ||
330 | 165 | !Allocate trace variable for upwind fluxes | ||
331 | 166 | lambda_mesh=>extract_mesh(state, "VelocityMeshTrace") | ||
332 | 167 | call allocate(UpwindFlux,lambda_mesh, "UpwindFlux") | ||
333 | 168 | call zero(UpwindFlux) | ||
334 | 169 | !Allocate matrix that maps T values to upwind fluxes | ||
335 | 170 | sparsity => get_csr_sparsity_firstorder(state,lambda_mesh,T%mesh) | ||
336 | 171 | call allocate(UpwindFluxMatrix,sparsity) | ||
337 | 172 | call zero(UpwindFluxMatrix) | ||
338 | 173 | end if | ||
339 | 174 | |||
340 | 175 | mass_sparsity=make_sparsity_dg_mass(T%mesh) | ||
341 | 176 | call allocate(mass, mass_sparsity) | ||
342 | 177 | call zero(mass) | ||
343 | 178 | call allocate(inv_mass, mass_sparsity) | ||
344 | 179 | call zero(inv_mass) | ||
345 | 180 | |||
346 | 181 | ! Ensure delta_T inherits options from T. | ||
347 | 182 | call allocate(delta_T, T%mesh, "delta_T") | ||
348 | 183 | call zero(delta_T) | ||
349 | 184 | delta_T%option_path = T%option_path | ||
350 | 185 | if(present(Flux)) then | ||
351 | 186 | call allocate(delta_T_total, T%mesh, "deltaT_total") | ||
352 | 187 | call zero(delta_T_total) | ||
353 | 188 | end if | ||
354 | 189 | call allocate(rhs, T%mesh, trim(field_name)//" RHS") | ||
355 | 190 | call zero(rhs) | ||
356 | 191 | |||
357 | 192 | if(present(Flux)) then | ||
358 | 193 | call construct_advection_dg(matrix, rhs, field_name, state, & | ||
359 | 194 | mass, velocity_name=velocity_name,continuity=continuity, & | ||
360 | 195 | UpwindFlux=upwindflux, UpwindFluxMatrix=UpwindFluxMatrix) | ||
361 | 196 | else | ||
362 | 197 | call construct_advection_dg(matrix, rhs, field_name, state, & | ||
363 | 198 | mass, velocity_name=velocity_name,continuity=continuity) | ||
364 | 199 | end if | ||
365 | 200 | |||
366 | 201 | call get_dg_inverse_mass_matrix(inv_mass, mass) | ||
367 | 202 | |||
368 | 203 | ! Note that since dt is module global, these lines have to | ||
369 | 204 | ! come after construct_advection_diffusion_dg. | ||
370 | 205 | call get_option("/timestepping/timestep", dt) | ||
371 | 206 | |||
372 | 207 | if(have_option(trim(T%option_path)//& | ||
373 | 208 | &"/prognostic/temporal_discretisation/discontinuous_galerkin/& | ||
374 | 209 | &/number_advection_subcycles")) then | ||
375 | 210 | call get_option(trim(T%option_path)//& | ||
376 | 211 | &"/prognostic/temporal_discretisation/discontinuous_galerkin/& | ||
377 | 212 | &/number_advection_subcycles", subcycles) | ||
378 | 213 | else | ||
379 | 214 | call get_option(trim(T%option_path)//& | ||
380 | 215 | &"/prognostic/temporal_discretisation/discontinuous_galerkin/& | ||
381 | 216 | &maximum_courant_number_per_subcycle", Max_Courant_number) | ||
382 | 217 | |||
383 | 218 | s_field => extract_scalar_field(state, "DG_CourantNumber") | ||
384 | 219 | call calculate_diagnostic_variable(state, "DG_CourantNumber_Local", & | ||
385 | 220 | & s_field) | ||
386 | 221 | ewrite_minmax(s_field) | ||
387 | 222 | subcycles = ceiling( maxval(s_field%val)/Max_Courant_number) | ||
388 | 223 | call allmax(subcycles) | ||
389 | 224 | ewrite(2,*) 'Number of subcycles for tracer eqn: ', subcycles | ||
390 | 225 | end if | ||
391 | 226 | |||
392 | 227 | limit_slope=.false. | ||
393 | 228 | if (have_option(trim(T%option_path)//"/prognostic/spatial_discretisation/& | ||
394 | 229 | &discontinuous_galerkin/slope_limiter")) then | ||
395 | 230 | limit_slope=.true. | ||
396 | 231 | |||
397 | 232 | call get_option(trim(T%option_path)//"/prognostic/spatial_discretisation/discontinuous_galerkin/slope_limiter/name",limiter_name) | ||
398 | 233 | |||
399 | 234 | select case(trim(limiter_name)) | ||
400 | 235 | case("Vertex_Based") | ||
401 | 236 | limiter=LIMITER_VB | ||
402 | 237 | case default | ||
403 | 238 | FLAbort('No such limiter') | ||
404 | 239 | end select | ||
405 | 240 | |||
406 | 241 | ! Note unsafe for mixed element meshes | ||
407 | 242 | if (element_degree(T,1)==0) then | ||
408 | 243 | FLExit("Slope limiters make no sense for degree 0 fields") | ||
409 | 244 | end if | ||
410 | 245 | |||
411 | 246 | end if | ||
412 | 247 | |||
413 | 248 | if (limit_slope) then | ||
414 | 249 | ! Filter wiggles from T | ||
415 | 250 | call limit_vb_manifold(state,T,delta_T) | ||
416 | 251 | if(present(flux)) then | ||
417 | 252 | call zero(upwindflux) | ||
418 | 253 | call mult(delta_T_total, mass, delta_T) | ||
419 | 254 | call update_flux(Flux, Delta_T_total, UpwindFlux, Mass) | ||
420 | 255 | end if | ||
421 | 256 | end if | ||
422 | 257 | |||
423 | 258 | do i=1, subcycles | ||
424 | 259 | ewrite(1,*) 'SUBCYCLE', i | ||
425 | 260 | |||
426 | 261 | ! dT = Advection * T | ||
427 | 262 | call mult(delta_T, matrix, T) | ||
428 | 263 | ! dT = dT + RHS | ||
429 | 264 | call addto(delta_T, RHS, -1.0) | ||
430 | 265 | if(present(Flux)) then | ||
431 | 266 | if(maxval(abs(RHS%val))>1.0e-8) then | ||
432 | 267 | FLAbort('Flux reconstruction doesn''t work if diffusion/bcs present at the moment') | ||
433 | 268 | end if | ||
434 | 269 | end if | ||
435 | 270 | call scale(delta_T, -dt/subcycles) | ||
436 | 271 | ! dT = M^(-1) dT | ||
437 | 272 | if(present(flux)) then | ||
438 | 273 | call mult(UpwindFlux, upwindfluxmatrix, T) | ||
439 | 274 | call scale(UpwindFlux,-dt/subcycles) | ||
440 | 275 | call update_flux(Flux, Delta_T, UpwindFlux, Mass) | ||
441 | 276 | end if | ||
442 | 277 | ! T = T + dt/s * dT | ||
443 | 278 | call dg_apply_mass(inv_mass, delta_T) | ||
444 | 279 | ewrite(2,*) 'Delta_T', maxval(abs(delta_T%val)) | ||
445 | 280 | call addto(T, delta_T) | ||
446 | 281 | !Probably need to halo_update(Flux) as well | ||
447 | 282 | call halo_update(T) | ||
448 | 283 | |||
449 | 284 | if (limit_slope) then | ||
450 | 285 | ! Filter wiggles from T | ||
451 | 286 | call limit_vb_manifold(state,T,delta_T) | ||
452 | 287 | if(present(flux)) then | ||
453 | 288 | call mult(delta_t_total, mass, delta_t) | ||
454 | 289 | call zero(UpwindFlux) | ||
455 | 290 | call update_flux(Flux, Delta_T_total, UpwindFlux, Mass) | ||
456 | 291 | end if | ||
457 | 292 | end if | ||
458 | 293 | end do | ||
459 | 294 | |||
460 | 295 | if(present(Flux)) then | ||
461 | 296 | if(have_option('/material_phase::Fluid/scalar_field::LayerThickness/p& | ||
462 | 297 | &rognostic/spatial_discretisation/debug')) then | ||
463 | 298 | do ele = 1, ele_count(Flux) | ||
464 | 299 | call check_flux(Flux,T,T_old,X,mass,ele) | ||
465 | 300 | end do | ||
466 | 301 | end if | ||
467 | 302 | call deallocate(UpwindFlux) | ||
468 | 303 | call deallocate(UpwindFluxMatrix) | ||
469 | 304 | call deallocate(delta_T_total) | ||
470 | 305 | end if | ||
471 | 306 | |||
472 | 307 | call deallocate(delta_T) | ||
473 | 308 | call deallocate(matrix) | ||
474 | 309 | call deallocate(mass) | ||
475 | 310 | call deallocate(inv_mass) | ||
476 | 311 | call deallocate(mass_sparsity) | ||
477 | 312 | call deallocate(rhs) | ||
478 | 313 | |||
479 | 314 | end subroutine solve_advection_dg_subcycle | ||
480 | 315 | |||
481 | 316 | subroutine check_flux(Flux,T,T_old,X,mass,ele) | ||
482 | 317 | type(vector_field), intent(in) :: Flux, X | ||
483 | 318 | type(scalar_field), intent(in) :: T,T_old | ||
484 | 319 | integer, intent(in) :: ele | ||
485 | 320 | type(csr_matrix), intent(in) :: mass | ||
486 | 321 | ! | ||
487 | 322 | real, dimension(Flux%dim,ele_loc(Flux,ele)) :: Flux_vals | ||
488 | 323 | real, dimension(ele_ngi(T,ele)) :: Div_Flux_gi | ||
489 | 324 | real, dimension(ele_ngi(T,ele)) :: Delta_T_gi, T_gi | ||
490 | 325 | real, dimension(ele_loc(T,ele)) :: Delta_T_rhs, Div_Flux_rhs,T_rhs | ||
491 | 326 | integer :: dim1, loc | ||
492 | 327 | real, dimension(ele_ngi(X,ele)) :: detwei | ||
493 | 328 | real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J_mat | ||
494 | 329 | real :: residual | ||
495 | 330 | Delta_T_gi = ele_val_at_quad(T,ele)-ele_val_at_quad(T_old,ele) | ||
496 | 331 | T_gi = ele_val_at_quad(T,ele) | ||
497 | 332 | Flux_vals = ele_val(Flux,ele) | ||
498 | 333 | |||
499 | 334 | Div_Flux_gi = 0. | ||
500 | 335 | !dn is loc x ngi x dim | ||
501 | 336 | do dim1 = 1, Flux%dim | ||
502 | 337 | do loc = 1, ele_loc(Flux,ele) | ||
503 | 338 | Div_Flux_gi = Div_Flux_gi + Flux%mesh%shape%dn(loc,:,dim1)*& | ||
504 | 339 | &Flux_vals(dim1,loc) | ||
505 | 340 | end do | ||
506 | 341 | end do | ||
507 | 342 | |||
508 | 343 | call compute_jacobian(ele_val(X,ele), ele_shape(X,ele),& | ||
509 | 344 | detwei=detwei,J=J_mat) | ||
510 | 345 | |||
511 | 346 | Delta_T_rhs = shape_rhs(ele_shape(T,ele),Delta_T_gi*detwei) | ||
512 | 347 | T_rhs = shape_rhs(ele_shape(T,ele),T_gi*detwei) | ||
513 | 348 | Div_Flux_rhs = shape_rhs(ele_shape(T,ele),Div_Flux_gi*& | ||
514 | 349 | Flux%mesh%shape%quadrature%weight) | ||
515 | 350 | |||
516 | 351 | residual = & | ||
517 | 352 | maxval(abs(Delta_T_rhs-Div_Flux_rhs))/max(1.0& | ||
518 | 353 | &,maxval(abs(T_rhs))) | ||
519 | 354 | if(residual>1.0e-6) then | ||
520 | 355 | ewrite(2,*) residual, maxval(abs(T_rhs)) | ||
521 | 356 | FLExit('Flux residual error.') | ||
522 | 357 | end if | ||
523 | 358 | |||
524 | 359 | end subroutine check_flux | ||
525 | 360 | |||
526 | 361 | subroutine update_flux(Flux, Delta_T, UpwindFlux, Mass) | ||
527 | 362 | !! Subroutine to compute div-conforming Flux such that | ||
528 | 363 | !! div Flux = Delta_T | ||
529 | 364 | !! with Flux.n = UpwindFlux on element boundaries | ||
530 | 365 | !! and then to add to Flux | ||
531 | 366 | type(vector_field), intent(inout) :: Flux | ||
532 | 367 | type(scalar_field), intent(in) :: Delta_T | ||
533 | 368 | type(scalar_field), intent(in) :: UpwindFlux | ||
534 | 369 | type(csr_matrix), intent(in) :: Mass | ||
535 | 370 | ! | ||
536 | 371 | integer :: ele, i | ||
537 | 372 | real :: Area | ||
538 | 373 | integer, dimension(:), pointer :: T_ele | ||
539 | 374 | |||
540 | 375 | !Checking the Flux is in local representation | ||
541 | 376 | assert(Flux%dim==mesh_dim(Flux)) | ||
542 | 377 | |||
543 | 378 | do ele = 1, ele_count(Flux) | ||
544 | 379 | T_ele => ele_nodes(delta_T,ele) | ||
545 | 380 | area = 0. | ||
546 | 381 | do i = 1, size(T_ele) | ||
547 | 382 | area = area + & | ||
548 | 383 | sum(row_val(Mass,t_ele(i))) | ||
549 | 384 | end do | ||
550 | 385 | call update_flux_ele(Flux, Delta_T, UpwindFlux,ele,area) | ||
551 | 386 | end do | ||
552 | 387 | |||
553 | 388 | end subroutine update_flux | ||
554 | 389 | |||
555 | 390 | subroutine update_flux_ele(Flux, Delta_T, UpwindFlux,ele,area) | ||
556 | 391 | !! Subroutine to compute div-conforming Flux such that | ||
557 | 392 | !! div Flux = Delta_T | ||
558 | 393 | !! with Flux.n = UpwindFlux on element boundaries | ||
559 | 394 | type(vector_field), intent(inout) :: Flux | ||
560 | 395 | type(scalar_field), intent(in) :: Delta_T | ||
561 | 396 | type(scalar_field), intent(in) :: UpwindFlux | ||
562 | 397 | integer, intent(in) :: ele | ||
563 | 398 | real, intent(in) :: area | ||
564 | 399 | ! | ||
565 | 400 | real, dimension(Flux%dim,ele_loc(Flux,ele)) :: Flux_vals | ||
566 | 401 | real, dimension(ele_loc(Delta_T,ele)) :: Delta_T_vals | ||
567 | 402 | real, dimension(Flux%dim*ele_loc(Flux,ele),& | ||
568 | 403 | Flux%dim*ele_loc(Flux,ele)) :: Flux_mat | ||
569 | 404 | real, dimension(Flux%dim*ele_loc(Flux,ele)) :: Flux_rhs | ||
570 | 405 | integer :: row, ni, ele2, face, face2,floc, i, dim1, dim2 | ||
571 | 406 | integer, dimension(:), pointer :: neigh | ||
572 | 407 | type(constraints_type), pointer :: flux_constraint | ||
573 | 408 | type(element_type), pointer :: flux_shape, T_shape | ||
574 | 409 | real, dimension(Flux%dim,ele_loc(Delta_T,ele),& | ||
575 | 410 | &ele_loc(Flux,ele)) :: grad_mat | ||
576 | 411 | real, dimension(ele_loc(Delta_t,ele)) :: delta_t_val | ||
577 | 412 | real, dimension(ele_loc(flux,ele),ele_loc(flux,ele)) & | ||
578 | 413 | &:: flux_mass_mat | ||
579 | 414 | real, dimension(ele_loc(Delta_t,ele)) :: div_flux_val | ||
580 | 415 | real, dimension(ele_ngi(Delta_t,ele)) :: div_flux_gi | ||
581 | 416 | real :: total_flux, residual | ||
582 | 417 | !! Need to set up and solve equation for flux DOFs | ||
583 | 418 | !! Rows are as follows: | ||
584 | 419 | !! All of the normal fluxes through edges, numbered according to | ||
585 | 420 | !! UpwindFlux numbering | ||
586 | 421 | !! Inner product of Flux with grad basis equaling gradient of delta_T | ||
587 | 422 | !! plus boundary conditions | ||
588 | 423 | !! Inner product of solution with curl basis equals zero | ||
589 | 424 | !! inner product of solution with constraint basis equals zero | ||
590 | 425 | |||
591 | 426 | !! Debugging check: | ||
592 | 427 | !! \int \Delta T dV = \int div Flux dV = \int Flux.n dS | ||
593 | 428 | !! = \int \hat{Flux}.\hat{n dS} | ||
594 | 429 | neigh => ele_neigh(Flux,ele) | ||
595 | 430 | total_flux = 0. | ||
596 | 431 | do ni = 1, size(neigh) | ||
597 | 432 | ele2 = neigh(ni) | ||
598 | 433 | face = ele_face(Flux,ele,ele2) | ||
599 | 434 | if(ele2>0) then | ||
600 | 435 | face2 = ele_face(Flux,ele2,ele) | ||
601 | 436 | else | ||
602 | 437 | face2 = face | ||
603 | 438 | end if | ||
604 | 439 | if(face>face2) then | ||
605 | 440 | total_flux=total_flux+sum(face_val(UpwindFlux,face)) | ||
606 | 441 | else | ||
607 | 442 | total_flux=total_flux-sum(face_val(UpwindFlux,face)) | ||
608 | 443 | end if | ||
609 | 444 | end do | ||
610 | 445 | residual = abs(total_flux-sum(ele_val(Delta_T,ele)))& | ||
611 | 446 | &/max(1.0,maxval(ele_val(Delta_T,ele)))/area | ||
612 | 447 | if(residual>1.0e-10) then | ||
613 | 448 | ewrite(0,*) 'Residual = ', residual | ||
614 | 449 | FLExit('Bad residual') | ||
615 | 450 | end if | ||
616 | 451 | |||
617 | 452 | Flux_rhs = 0. | ||
618 | 453 | Flux_mat = 0. | ||
619 | 454 | flux_shape => ele_shape(flux,ele) | ||
620 | 455 | T_shape => ele_shape(delta_t,ele) | ||
621 | 456 | flux_constraint => flux%mesh%shape%constraints | ||
622 | 457 | |||
623 | 458 | row = 1 | ||
624 | 459 | !first do the face DOFs | ||
625 | 460 | neigh => ele_neigh(Flux,ele) | ||
626 | 461 | do ni = 1, size(neigh) | ||
627 | 462 | ele2 = neigh(ni) | ||
628 | 463 | face = ele_face(Flux,ele,ele2) | ||
629 | 464 | if(ele2>0) then | ||
630 | 465 | face2 = ele_face(Flux,ele2,ele) | ||
631 | 466 | else | ||
632 | 467 | face2 = face | ||
633 | 468 | end if | ||
634 | 469 | floc = face_loc(upwindflux,face) | ||
635 | 470 | |||
636 | 471 | call update_flux_face(& | ||
637 | 472 | Flux,UpwindFlux,face,face2,& | ||
638 | 473 | ele,ele2,Flux_mat(row:row+floc-1,:),& | ||
639 | 474 | Flux_rhs(row:row+floc-1)) | ||
640 | 475 | row = row + floc | ||
641 | 476 | end do | ||
642 | 477 | |||
643 | 478 | !Next the gradient basis | ||
644 | 479 | !Start with \nabla\cdot F = \Delta T | ||
645 | 480 | !Integrate with test function over element | ||
646 | 481 | ! <\phi,Div F>_E = <\phi, \Delta T>_E | ||
647 | 482 | |||
648 | 483 | grad_mat = shape_dshape(T_shape, flux_shape%dn,& | ||
649 | 484 | & T_shape%quadrature%weight) | ||
650 | 485 | delta_t_val = ele_val(delta_t,ele) | ||
651 | 486 | do i = 1, ele_loc(Delta_T,ele)-1 | ||
652 | 487 | do dim1 = 1, flux%dim | ||
653 | 488 | Flux_mat(row,(dim1-1)*ele_loc(flux,ele)+1:dim1*ele_loc(flux,ele))& | ||
654 | 489 | &=grad_mat(dim1,i,:) | ||
655 | 490 | end do | ||
656 | 491 | flux_rhs(row) = delta_t_val(i) | ||
657 | 492 | row = row + 1 | ||
658 | 493 | end do | ||
659 | 494 | |||
660 | 495 | !Then the curl basis | ||
661 | 496 | !Adding in the curl terms | ||
662 | 497 | flux_mass_mat = shape_shape(flux_shape,flux_shape,T_shape%quadrature%weight) | ||
663 | 498 | do i = 1, flux_constraint%n_curl_basis | ||
664 | 499 | flux_mat(row,:) = 0. | ||
665 | 500 | flux_rhs(row) = 0. | ||
666 | 501 | do dim1 = 1, mesh_dim(flux) | ||
667 | 502 | do dim2 = 1, mesh_dim(flux) | ||
668 | 503 | flux_mat(& | ||
669 | 504 | row,(dim2-1)*ele_loc(flux,ele)+1:dim2*ele_loc(flux,ele)) =& | ||
670 | 505 | flux_mat(& | ||
671 | 506 | row,(dim2-1)*ele_loc(flux,ele)+1:dim2*ele_loc(flux,ele)) +& | ||
672 | 507 | matmul(flux_constraint%curl_basis(i,:,dim1),& | ||
673 | 508 | flux_mass_mat) | ||
674 | 509 | end do | ||
675 | 510 | end do | ||
676 | 511 | row = row + 1 | ||
677 | 512 | end do | ||
678 | 513 | |||
679 | 514 | !Then the constraints | ||
680 | 515 | do i = 1, flux_constraint%n_constraints | ||
681 | 516 | do dim1 = 1, mesh_dim(flux) | ||
682 | 517 | flux_mat(& | ||
683 | 518 | row,(dim1-1)*ele_loc(flux,ele)+1:& | ||
684 | 519 | dim1*ele_loc(flux,ele)) =& | ||
685 | 520 | flux_constraint%orthogonal(i,:,dim1) | ||
686 | 521 | end do | ||
687 | 522 | row = row + 1 | ||
688 | 523 | end do | ||
689 | 524 | |||
690 | 525 | call solve(flux_mat,flux_rhs) | ||
691 | 526 | |||
692 | 527 | do dim1 = 1, mesh_dim(flux) | ||
693 | 528 | flux_vals(dim1,:) = & | ||
694 | 529 | flux_rhs((dim1-1)*ele_loc(flux,ele)+1:dim1*ele_loc(flux,ele)) | ||
695 | 530 | end do | ||
696 | 531 | |||
697 | 532 | !A debugging check. | ||
698 | 533 | !Computing < \phi, div F> in local coordinates | ||
699 | 534 | !and comparing with <\phi, delta_T> in global coordinates | ||
700 | 535 | !(works because of pullback formula) | ||
701 | 536 | |||
702 | 537 | !dn is loc x ngi x dim | ||
703 | 538 | !Flux is dim x loc | ||
704 | 539 | div_flux_gi = 0. | ||
705 | 540 | do dim1 =1, size(flux_vals,1) | ||
706 | 541 | do i=1,size(flux_vals,2) | ||
707 | 542 | div_flux_gi = div_flux_gi + flux_vals(dim1,i)*& | ||
708 | 543 | &flux_shape%dn(i,:,dim1) | ||
709 | 544 | end do | ||
710 | 545 | end do | ||
711 | 546 | div_flux_val = shape_rhs(T_shape,div_flux_gi*T_shape%quadrature%weight) | ||
712 | 547 | residual = maxval(abs(div_flux_val-delta_T_val))/max(1.0& | ||
713 | 548 | &,maxval(abs(delta_T_val)))/area | ||
714 | 549 | assert(residual<1.0e-10) | ||
715 | 550 | |||
716 | 551 | !Add the flux in | ||
717 | 552 | call addto(flux,ele_nodes(flux,ele),flux_vals) | ||
718 | 553 | end subroutine update_flux_ele | ||
719 | 554 | |||
720 | 555 | subroutine update_flux_face(& | ||
721 | 556 | Flux,UpwindFlux,face,face2,& | ||
722 | 557 | ele,ele2,Flux_mat_rows,Flux_rhs_rows) | ||
723 | 558 | type(vector_field), intent(in) :: Flux | ||
724 | 559 | type(scalar_Field), intent(in) :: UpwindFlux | ||
725 | 560 | integer, intent(in) :: face,face2,ele,ele2 | ||
726 | 561 | real, dimension(face_loc(upwindflux,ele),& | ||
727 | 562 | &mesh_dim(flux)*ele_loc(flux,ele)), & | ||
728 | 563 | &intent(inout) :: flux_mat_rows | ||
729 | 564 | real, dimension(flux%mesh%shape%constraints%n_face_basis), intent(inout)& | ||
730 | 565 | &:: flux_rhs_rows | ||
731 | 566 | ! | ||
732 | 567 | integer :: dim1, row, flux_dim1 | ||
733 | 568 | real, dimension(mesh_dim(flux), face_ngi(flux, face)) :: n_local | ||
734 | 569 | real, dimension(face_ngi(flux,face)) :: detwei_f | ||
735 | 570 | real :: weight | ||
736 | 571 | type(element_type), pointer :: flux_face_shape, upwindflux_face_shape | ||
737 | 572 | integer, dimension(face_loc(flux,face)) :: flux_face_nodes | ||
738 | 573 | real, dimension(mesh_dim(flux),face_loc(upwindflux,face),& | ||
739 | 574 | face_loc(flux,face)) :: face_mat | ||
740 | 575 | ! | ||
741 | 576 | flux_face_shape=>face_shape(flux, face) | ||
742 | 577 | upwindflux_face_shape=>face_shape(upwindflux, face) | ||
743 | 578 | flux_face_nodes=face_local_nodes(flux, face) | ||
744 | 579 | |||
745 | 580 | !Get normal in local coordinates | ||
746 | 581 | call get_local_normal(n_local,weight,flux,& | ||
747 | 582 | &local_face_number(flux%mesh,face)) | ||
748 | 583 | |||
749 | 584 | detwei_f = weight*flux_face_shape%quadrature%weight | ||
750 | 585 | face_mat = shape_shape_vector(& | ||
751 | 586 | upwindflux_face_shape,flux_face_shape,detwei_f,n_local) | ||
752 | 587 | !Equation is: | ||
753 | 588 | ! <\phi,Flux.n> = <\phi,UpwindFlux> for all trace test functions \phi. | ||
754 | 589 | |||
755 | 590 | do row = 1, size(flux_mat_rows,1) | ||
756 | 591 | flux_mat_rows(row,:) = 0. | ||
757 | 592 | do dim1 = 1,mesh_dim(flux) | ||
758 | 593 | flux_dim1 = (dim1-1)*ele_loc(flux,ele) | ||
759 | 594 | flux_mat_rows(row,flux_dim1+flux_face_nodes)& | ||
760 | 595 | &=face_mat(dim1,row,:) | ||
761 | 596 | end do | ||
762 | 597 | end do | ||
763 | 598 | !Need to adopt a sign convention: | ||
764 | 599 | !If face>face_2 | ||
765 | 600 | !We store flux with normal pointing from face into face_2 | ||
766 | 601 | !Otherwise | ||
767 | 602 | !We store flux with normal pointing from face_2 into face | ||
768 | 603 | ! Outflow boundary integral. | ||
769 | 604 | if(face>face2) then | ||
770 | 605 | flux_rhs_rows = face_val(UpwindFlux,face) | ||
771 | 606 | else | ||
772 | 607 | flux_rhs_rows = -face_val(UpwindFlux,face) | ||
773 | 608 | end if | ||
774 | 609 | end subroutine update_flux_face | ||
775 | 610 | |||
776 | 611 | subroutine construct_advection_dg(big_m, rhs, field_name,& | ||
777 | 612 | & state, mass, velocity_name, continuity,& | ||
778 | 613 | & UpwindFlux, UpwindFluxMatrix) | ||
779 | 614 | !!< Construct the advection equation for discontinuous elements in | ||
780 | 615 | !!< acceleration form. | ||
781 | 616 | !!< | ||
782 | 617 | !!< If mass is provided then the mass matrix is not added into big_m or | ||
783 | 618 | !!< rhs. It is instead returned as mass. This may be useful for testing | ||
784 | 619 | !!< or for solving equations otherwise than in acceleration form. | ||
785 | 620 | |||
786 | 621 | !! Main advection matrix. | ||
787 | 622 | type(csr_matrix), intent(inout) :: big_m | ||
788 | 623 | !! Right hand side vector. | ||
789 | 624 | type(scalar_field), intent(inout) :: rhs | ||
790 | 625 | |||
791 | 626 | !! Name of the field to be advected. | ||
792 | 627 | character(len=*), intent(in) :: field_name | ||
793 | 628 | !! Collection of fields defining system state. | ||
794 | 629 | type(state_type), intent(inout) :: state | ||
795 | 630 | !! Optional separate mass matrix. | ||
796 | 631 | type(csr_matrix), intent(inout), optional :: mass | ||
797 | 632 | !! Optional velocity name | ||
798 | 633 | character(len = *), intent(in), optional :: velocity_name | ||
799 | 634 | !! solve the continuity equation | ||
800 | 635 | logical, intent(in), optional :: continuity | ||
801 | 636 | !! Upwind flux field | ||
802 | 637 | type(scalar_field), intent(in), optional :: UpwindFlux | ||
803 | 638 | !! Upwind flux matrix | ||
804 | 639 | type(csr_matrix), intent(inout), optional :: UpwindFluxMatrix | ||
805 | 640 | |||
806 | 641 | !! Position, and velocity fields. | ||
807 | 642 | type(vector_field) :: X, U, U_nl | ||
808 | 643 | !! Tracer to be solved for. | ||
809 | 644 | type(scalar_field) :: T | ||
810 | 645 | |||
811 | 646 | !! Local velocity name | ||
812 | 647 | character(len = FIELD_NAME_LEN) :: lvelocity_name | ||
813 | 648 | |||
814 | 649 | !! Element index | ||
815 | 650 | integer :: ele | ||
816 | 651 | |||
817 | 652 | !! Status variable for field extraction. | ||
818 | 653 | integer :: stat | ||
819 | 654 | |||
820 | 655 | ewrite(1,*) "Writing advection equation for "& | ||
821 | 656 | &//trim(field_name) | ||
822 | 657 | |||
823 | 658 | ! These names are based on the CGNS SIDS. | ||
824 | 659 | T=extract_scalar_field(state, field_name) | ||
825 | 660 | X=extract_vector_field(state, "Coordinate") | ||
826 | 661 | |||
827 | 662 | if(present(velocity_name)) then | ||
828 | 663 | lvelocity_name = velocity_name | ||
829 | 664 | else | ||
830 | 665 | lvelocity_name = "NonlinearVelocity" | ||
831 | 666 | end if | ||
832 | 667 | |||
833 | 668 | if (.not.have_option(trim(T%option_path)//"/prognostic& | ||
834 | 669 | &/spatial_discretisation/discontinuous_galerkin& | ||
835 | 670 | &/advection_scheme/none")) then | ||
836 | 671 | U_nl=extract_vector_field(state, lvelocity_name) | ||
837 | 672 | call incref(U_nl) | ||
838 | 673 | include_advection=.true. | ||
839 | 674 | else | ||
840 | 675 | ! Forcing a zero NonlinearVelocity will disable advection. | ||
841 | 676 | U=extract_vector_field(state, "Velocity", stat) | ||
842 | 677 | if (stat/=0) then | ||
843 | 678 | FLExit("Oh dear, no velocity field. A velocity field is required for advection!") | ||
844 | 679 | end if | ||
845 | 680 | call allocate(U_nl, U%dim, U%mesh, "LocalNonlinearVelocity") | ||
846 | 681 | call zero(U_nl) | ||
847 | 682 | include_advection=.false. | ||
848 | 683 | end if | ||
849 | 684 | |||
850 | 685 | assert(has_faces(X%mesh)) | ||
851 | 686 | assert(has_faces(T%mesh)) | ||
852 | 687 | |||
853 | 688 | call zero(big_m) | ||
854 | 689 | call zero(RHS) | ||
855 | 690 | call zero(mass) | ||
856 | 691 | |||
857 | 692 | element_loop: do ele=1,element_count(T) | ||
858 | 693 | |||
859 | 694 | call construct_adv_element_dg(ele, big_m, rhs,& | ||
860 | 695 | & X, T, U_nl, mass, continuity=continuity,& | ||
861 | 696 | & upwindflux=upwindflux, upwindfluxmatrix=upwindfluxmatrix) | ||
862 | 697 | |||
863 | 698 | end do element_loop | ||
864 | 699 | |||
865 | 700 | ! Drop any extra field references. | ||
866 | 701 | |||
867 | 702 | call deallocate(U_nl) | ||
868 | 703 | |||
869 | 704 | end subroutine construct_advection_dg | ||
870 | 705 | |||
871 | 706 | subroutine construct_adv_element_dg(ele, big_m, rhs,& | ||
872 | 707 | & X, T, U_nl, mass, continuity, upwindflux, upwindfluxmatrix) | ||
873 | 708 | !!< Construct the advection_diffusion equation for discontinuous elements in | ||
874 | 709 | !!< acceleration form. | ||
875 | 710 | implicit none | ||
876 | 711 | !! Index of current element | ||
877 | 712 | integer :: ele | ||
878 | 713 | !! Main advection matrix. | ||
879 | 714 | type(csr_matrix), intent(inout) :: big_m | ||
880 | 715 | !! Right hand side vector. | ||
881 | 716 | type(scalar_field), intent(inout) :: rhs | ||
882 | 717 | !! Optional separate mass matrix. | ||
883 | 718 | type(csr_matrix), intent(inout) :: mass | ||
884 | 719 | |||
885 | 720 | !! Position and velocity. | ||
886 | 721 | type(vector_field), intent(in) :: X, U_nl | ||
887 | 722 | |||
888 | 723 | type(scalar_field), intent(in) :: T | ||
889 | 724 | |||
890 | 725 | !! Upwind flux field | ||
891 | 726 | type(scalar_field), intent(in), optional :: UpwindFlux | ||
892 | 727 | !! Upwind flux matrix | ||
893 | 728 | type(csr_matrix), intent(inout), optional :: UpwindFluxMatrix | ||
894 | 729 | |||
895 | 730 | |||
896 | 731 | !! Solve the continuity equation rather than advection | ||
897 | 732 | logical, intent(in), optional :: continuity | ||
898 | 733 | |||
899 | 734 | ! Bilinear forms. | ||
900 | 735 | real, dimension(ele_loc(T,ele), ele_loc(T,ele)) :: & | ||
901 | 736 | mass_mat | ||
902 | 737 | real, dimension(ele_loc(T,ele), ele_loc(T,ele)) :: & | ||
903 | 738 | Advection_mat, Ad_mat1, Ad_mat2 | ||
904 | 739 | |||
905 | 740 | ! Local assembly matrices. | ||
906 | 741 | real, dimension(ele_loc(T,ele), ele_loc(T,ele)) :: l_T_mat | ||
907 | 742 | real, dimension(ele_loc(T,ele)) :: l_T_rhs | ||
908 | 743 | |||
909 | 744 | ! Local variables. | ||
910 | 745 | |||
911 | 746 | ! Neighbour element, face and neighbour face. | ||
912 | 747 | integer :: ele_2, face, face_2 | ||
913 | 748 | ! Loops over faces. | ||
914 | 749 | integer :: ni | ||
915 | 750 | |||
916 | 751 | ! Transform from local to physical coordinates. | ||
917 | 752 | real, dimension(ele_ngi(X,ele)) :: detwei | ||
918 | 753 | real, dimension(mesh_dim(U_nl), X%dim, ele_ngi(X,ele)) :: J_mat | ||
919 | 754 | |||
920 | 755 | ! Different velocities at quad points. | ||
921 | 756 | real, dimension(U_nl%dim, ele_ngi(U_nl, ele)) :: U_quad | ||
922 | 757 | real, dimension(U_nl%dim, ele_ngi(U_nl, ele)) :: U_nl_q | ||
923 | 758 | real, dimension(ele_ngi(U_nl, ele)) :: U_nl_div_q | ||
924 | 759 | |||
925 | 760 | ! Node and shape pointers. | ||
926 | 761 | integer, dimension(:), pointer :: T_ele | ||
927 | 762 | type(element_type), pointer :: T_shape, U_shape | ||
928 | 763 | ! Neighbours of this element. | ||
929 | 764 | integer, dimension(:), pointer :: neigh | ||
930 | 765 | |||
931 | 766 | integer :: gi,i,j | ||
932 | 767 | |||
933 | 768 | !---------------------------------------------------------------------- | ||
934 | 769 | ! Establish local node lists | ||
935 | 770 | !---------------------------------------------------------------------- | ||
936 | 771 | |||
937 | 772 | T_ele=>ele_nodes(T,ele) ! Tracer node numbers | ||
938 | 773 | |||
939 | 774 | !---------------------------------------------------------------------- | ||
940 | 775 | ! Establish local shape functions | ||
941 | 776 | !---------------------------------------------------------------------- | ||
942 | 777 | |||
943 | 778 | T_shape=>ele_shape(T, ele) | ||
944 | 779 | U_shape=>ele_shape(U_nl, ele) | ||
945 | 780 | |||
946 | 781 | !========================== | ||
947 | 782 | ! Coordinates | ||
948 | 783 | !========================== | ||
949 | 784 | |||
950 | 785 | ! Get J_mat | ||
951 | 786 | call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), detwei=detwei, J=J_mat) | ||
952 | 787 | |||
953 | 788 | !---------------------------------------------------------------------- | ||
954 | 789 | ! Construct bilinear forms. | ||
955 | 790 | !---------------------------------------------------------------------- | ||
956 | 791 | |||
957 | 792 | ! Element density matrix. | ||
958 | 793 | ! / | ||
959 | 794 | ! | T T dV | ||
960 | 795 | ! / | ||
961 | 796 | mass_mat = shape_shape(T_shape, T_shape, detwei) | ||
962 | 797 | |||
963 | 798 | if (include_advection) then | ||
964 | 799 | |||
965 | 800 | ! Advecting velocity at quadrature points. | ||
966 | 801 | U_quad=ele_val_at_quad(U_nl,ele) | ||
967 | 802 | U_nl_q=U_quad | ||
968 | 803 | |||
969 | 804 | U_nl_div_q=ele_div_at_quad(U_nl, ele, U_shape%dn) | ||
970 | 805 | |||
971 | 806 | ! Element advection matrix | ||
972 | 807 | ! / / | ||
973 | 808 | ! - | (grad phi dot U_nl) T dV -| phi ( div U_nl ) T dV | ||
974 | 809 | ! / / | ||
975 | 810 | Advection_mat = -dshape_dot_vector_shape(T_shape%dn, U_nl_q, T_shape,& | ||
976 | 811 | & T_shape%quadrature%weight) | ||
977 | 812 | if(.not.present_and_true(continuity)) then | ||
978 | 813 | Advection_mat = Advection_mat& | ||
979 | 814 | &-shape_shape(T_shape, T_shape, U_nl_div_q * T_shape& | ||
980 | 815 | &%quadrature%weight) | ||
981 | 816 | end if | ||
982 | 817 | else | ||
983 | 818 | Advection_mat=0.0 | ||
984 | 819 | end if | ||
985 | 820 | |||
986 | 821 | !---------------------------------------------------------------------- | ||
987 | 822 | ! Perform global assembly. | ||
988 | 823 | !---------------------------------------------------------------------- | ||
989 | 824 | |||
990 | 825 | l_T_rhs=0.0 | ||
991 | 826 | |||
992 | 827 | ! Assemble matrix. | ||
993 | 828 | |||
994 | 829 | ! Advection. | ||
995 | 830 | l_T_mat= Advection_mat | ||
996 | 831 | |||
997 | 832 | call addto(mass, t_ele, t_ele, mass_mat) | ||
998 | 833 | |||
999 | 834 | call addto(big_m, t_ele, t_ele, l_T_mat) | ||
1000 | 835 | |||
1001 | 836 | !------------------------------------------------------------------- | ||
1002 | 837 | ! Interface integrals | ||
1003 | 838 | !------------------------------------------------------------------- | ||
1004 | 839 | |||
1005 | 840 | neigh=>ele_neigh(T, ele) | ||
1006 | 841 | |||
1007 | 842 | neighbourloop: do ni=1,size(neigh) | ||
1008 | 843 | |||
1009 | 844 | !---------------------------------------------------------------------- | ||
1010 | 845 | ! Find the relevant faces. | ||
1011 | 846 | !---------------------------------------------------------------------- | ||
1012 | 847 | |||
1013 | 848 | ! These finding routines are outside the inner loop so as to allow | ||
1014 | 849 | ! for local stack variables of the right size in | ||
1015 | 850 | ! construct_add_diff_interface_dg. | ||
1016 | 851 | |||
1017 | 852 | ele_2=neigh(ni) | ||
1018 | 853 | |||
1019 | 854 | ! Note that although face is calculated on field U, it is in fact | ||
1020 | 855 | ! applicable to any field which shares the same mesh topology. | ||
1021 | 856 | face=ele_face(T, ele, ele_2) | ||
1022 | 857 | |||
1023 | 858 | if (ele_2>0) then | ||
1024 | 859 | ! Internal faces. | ||
1025 | 860 | face_2=ele_face(T, ele_2, ele) | ||
1026 | 861 | else | ||
1027 | 862 | ! External face. | ||
1028 | 863 | face_2=face | ||
1029 | 864 | end if | ||
1030 | 865 | |||
1031 | 866 | call construct_adv_interface_dg(ele, face, face_2,& | ||
1032 | 867 | & big_m, rhs, X, T, U_nl, upwindflux, upwindfluxmatrix) | ||
1033 | 868 | |||
1034 | 869 | end do neighbourloop | ||
1035 | 870 | |||
1036 | 871 | end subroutine construct_adv_element_dg | ||
1037 | 872 | |||
1038 | 873 | subroutine construct_adv_interface_dg(ele, face, face_2, & | ||
1039 | 874 | big_m, rhs, X, T, U_nl,upwindflux, upwindfluxmatrix) | ||
1040 | 875 | |||
1041 | 876 | !!< Construct the DG element boundary integrals on the ni-th face of | ||
1042 | 877 | !!< element ele. | ||
1043 | 878 | implicit none | ||
1044 | 879 | |||
1045 | 880 | integer, intent(in) :: ele, face, face_2 | ||
1046 | 881 | type(csr_matrix), intent(inout) :: big_m | ||
1047 | 882 | type(scalar_field), intent(inout) :: rhs | ||
1048 | 883 | ! We pass these additional fields to save on state lookups. | ||
1049 | 884 | type(vector_field), intent(in) :: X, U_nl | ||
1050 | 885 | type(scalar_field), intent(in) :: T | ||
1051 | 886 | !! Upwind flux field | ||
1052 | 887 | type(scalar_field), intent(in), optional :: UpwindFlux | ||
1053 | 888 | !! Upwind flux matrix | ||
1054 | 889 | type(csr_matrix), intent(inout), optional :: UpwindFluxMatrix | ||
1055 | 890 | |||
1056 | 891 | ! Face objects and numberings. | ||
1057 | 892 | type(element_type), pointer :: T_shape, T_shape_2 | ||
1058 | 893 | integer, dimension(face_loc(T,face)) :: T_face | ||
1059 | 894 | integer, dimension(face_loc(T,face_2)) :: T_face_2 | ||
1060 | 895 | integer, dimension(face_loc(U_nl,face)) :: U_face | ||
1061 | 896 | integer, dimension(face_loc(U_nl,face_2)) :: U_face_2 | ||
1062 | 897 | |||
1063 | 898 | ! Note that both sides of the face can be assumed to have the same | ||
1064 | 899 | ! number of quadrature points. | ||
1065 | 900 | real, dimension(U_nl%dim, face_ngi(U_nl, face)) :: n1, n2, normal, U_nl_q,& | ||
1066 | 901 | & u_f_q, u_f2_q, div_u_f_q | ||
1067 | 902 | logical, dimension(face_ngi(U_nl, face)) :: inflow | ||
1068 | 903 | real, dimension(face_ngi(U_nl, face)) :: U_nl_q_dotn1, U_nl_q_dotn2, U_nl_q_dotn, income | ||
1069 | 904 | ! Variable transform times quadrature weights. | ||
1070 | 905 | real, dimension(face_ngi(T,face)) :: detwei | ||
1071 | 906 | real, dimension(U_nl%dim, X%dim, face_ngi(T,face)) :: J | ||
1072 | 907 | real, dimension(face_ngi(T,face)) :: inner_advection_integral, outer_advection_integral | ||
1073 | 908 | |||
1074 | 909 | ! Bilinear forms | ||
1075 | 910 | real, dimension(face_loc(T,face),face_loc(T,face)) :: nnAdvection_out | ||
1076 | 911 | real, dimension(face_loc(T,face),face_loc(T,face_2)) :: nnAdvection_in | ||
1077 | 912 | |||
1078 | 913 | ! Normal weights | ||
1079 | 914 | real :: w1, w2 | ||
1080 | 915 | integer :: i | ||
1081 | 916 | |||
1082 | 917 | T_face=face_global_nodes(T, face) | ||
1083 | 918 | T_shape=>face_shape(T, face) | ||
1084 | 919 | |||
1085 | 920 | T_face_2=face_global_nodes(T, face_2) | ||
1086 | 921 | T_shape_2=>face_shape(T, face_2) | ||
1087 | 922 | |||
1088 | 923 | !Get normals | ||
1089 | 924 | call get_local_normal(n1, w1, U_nl, local_face_number(U_nl%mesh,face)) | ||
1090 | 925 | call get_local_normal(n2, w2, U_nl, local_face_number(U_nl%mesh,face_2)) | ||
1091 | 926 | |||
1092 | 927 | if (include_advection) then | ||
1093 | 928 | |||
1094 | 929 | ! Advecting velocity at quadrature points. | ||
1095 | 930 | U_f_q = face_val_at_quad(U_nl, face) | ||
1096 | 931 | U_f2_q = face_val_at_quad(U_nl, face_2) | ||
1097 | 932 | |||
1098 | 933 | u_nl_q_dotn1 = sum(U_f_q*w1*n1,1) | ||
1099 | 934 | u_nl_q_dotn2 = -sum(U_f2_q*w2*n2,1) | ||
1100 | 935 | U_nl_q_dotn=0.5*(u_nl_q_dotn1+u_nl_q_dotn2) | ||
1101 | 936 | |||
1102 | 937 | ! Inflow is true if the flow at this gauss point is directed | ||
1103 | 938 | ! into this element. | ||
1104 | 939 | inflow = u_nl_q_dotn<0.0 | ||
1105 | 940 | income = merge(1.0,0.0,inflow) | ||
1106 | 941 | |||
1107 | 942 | !---------------------------------------------------------------------- | ||
1108 | 943 | ! Construct bilinear forms. | ||
1109 | 944 | !---------------------------------------------------------------------- | ||
1110 | 945 | |||
1111 | 946 | ! Calculate outflow boundary integral. | ||
1112 | 947 | ! can anyone think of a way of optimising this more to avoid | ||
1113 | 948 | ! superfluous operations (i.e. multiplying things by 0 or 1)? | ||
1114 | 949 | |||
1115 | 950 | ! first the integral around the inside of the element | ||
1116 | 951 | ! (this is the flux *out* of the element) | ||
1117 | 952 | inner_advection_integral = (1.-income)*u_nl_q_dotn | ||
1118 | 953 | nnAdvection_out=shape_shape(T_shape, T_shape, & | ||
1119 | 954 | & inner_advection_integral*T_shape%quadrature%weight) | ||
1120 | 955 | |||
1121 | 956 | ! now the integral around the outside of the element | ||
1122 | 957 | ! (this is the flux *in* to the element) | ||
1123 | 958 | outer_advection_integral = income*u_nl_q_dotn | ||
1124 | 959 | nnAdvection_in=shape_shape(T_shape, T_shape_2, & | ||
1125 | 960 | & outer_advection_integral*T_shape%quadrature%weight) | ||
1126 | 961 | |||
1127 | 962 | end if | ||
1128 | 963 | |||
1129 | 964 | !---------------------------------------------------------------------- | ||
1130 | 965 | ! Perform global assembly. | ||
1131 | 966 | !---------------------------------------------------------------------- | ||
1132 | 967 | |||
1133 | 968 | ! Insert advection in matrix. | ||
1134 | 969 | if (include_advection) then | ||
1135 | 970 | |||
1136 | 971 | ! Outflow boundary integral. | ||
1137 | 972 | call addto(big_M, T_face, T_face,& | ||
1138 | 973 | nnAdvection_out) | ||
1139 | 974 | |||
1140 | 975 | ! Inflow boundary integral. | ||
1141 | 976 | call addto(big_M, T_face, T_face_2,& | ||
1142 | 977 | nnAdvection_in) | ||
1143 | 978 | |||
1144 | 979 | if(present(UpwindFlux).and.present(UpwindFluxMatrix)) then | ||
1145 | 980 | call construct_upwindflux_interface(upwindflux,upwindfluxmatrix) | ||
1146 | 981 | end if | ||
1147 | 982 | end if | ||
1148 | 983 | |||
1149 | 984 | contains | ||
1150 | 985 | subroutine construct_upwindflux_interface(upwindflux,upwindfluxmatrix) | ||
1151 | 986 | type(scalar_field), intent(in) :: upwindflux | ||
1152 | 987 | type(csr_matrix), intent(inout) :: upwindfluxmatrix | ||
1153 | 988 | ! | ||
1154 | 989 | ! Bilinear forms | ||
1155 | 990 | real, dimension(face_loc(upwindflux,face),& | ||
1156 | 991 | face_loc(T,face)) :: upwindflux_mat_in,upwindflux_mat_out | ||
1157 | 992 | type(element_type), pointer :: flux_shape | ||
1158 | 993 | integer, dimension(face_loc(upwindflux,face)) :: flux_face | ||
1159 | 994 | integer :: sign | ||
1160 | 995 | |||
1161 | 996 | flux_shape=>face_shape(upwindflux, face) | ||
1162 | 997 | flux_face=face_global_nodes(upwindflux, face) | ||
1163 | 998 | |||
1164 | 999 | upwindflux_mat_in=shape_shape(flux_shape, T_shape, & | ||
1165 | 1000 | income*u_nl_q_dotn*T_shape%quadrature%weight) | ||
1166 | 1001 | upwindflux_mat_out=shape_shape(flux_shape, T_shape, & | ||
1167 | 1002 | (1.0-income)*u_nl_q_dotn*T_shape%quadrature%weight) | ||
1168 | 1003 | |||
1169 | 1004 | !Need to adopt a sign convention: | ||
1170 | 1005 | !If face>face_2 | ||
1171 | 1006 | !We store flux with normal pointing from face into face_2 | ||
1172 | 1007 | !Otherwise | ||
1173 | 1008 | !We store flux with normal pointing from face_2 into face | ||
1174 | 1009 | ! Outflow boundary integral. | ||
1175 | 1010 | |||
1176 | 1011 | !inflow = u_nl_q_dotn<0.0 | ||
1177 | 1012 | !income = 1 if(inflow) and 0 otherwise | ||
1178 | 1013 | if(face>face_2) then | ||
1179 | 1014 | sign = 1 | ||
1180 | 1015 | else | ||
1181 | 1016 | sign = -1 | ||
1182 | 1017 | end if | ||
1183 | 1018 | |||
1184 | 1019 | !Factor of 0.5 is because we visit from both sides | ||
1185 | 1020 | if(face.ne.face_2) then | ||
1186 | 1021 | call addto(upwindfluxmatrix, flux_face, T_face_2,& | ||
1187 | 1022 | 0.5*sign*upwindflux_mat_in) | ||
1188 | 1023 | end if | ||
1189 | 1024 | call addto(upwindfluxmatrix, flux_face, T_face,& | ||
1190 | 1025 | 0.5*sign*upwindflux_mat_out) | ||
1191 | 1026 | |||
1192 | 1027 | end subroutine construct_upwindflux_interface | ||
1193 | 1028 | |||
1194 | 1029 | end subroutine construct_adv_interface_dg | ||
1195 | 1030 | |||
1196 | 1031 | subroutine solve_advection_cg_tracer(Q,D,D_old,Flux,U_nl,QF,state) | ||
1197 | 1032 | !!< Solve the continuity equation for D*Q with Flux | ||
1198 | 1033 | !!< d/dt (D*Q) + div(Flux*Q) = diffusion terms. | ||
1199 | 1034 | !!< Done on the vorticity mesh | ||
1200 | 1035 | !!< Return a PV flux Q defined on the velocity mesh | ||
1201 | 1036 | !!< which is the projection of Flux*Q into the velocity space | ||
1202 | 1037 | !!< Note that Flux and Q contain a factor of dt for convenience. | ||
1203 | 1038 | type(state_type), intent(inout) :: state | ||
1204 | 1039 | type(scalar_field), intent(in),target :: D,D_old | ||
1205 | 1040 | type(scalar_field), intent(inout),target :: Q | ||
1206 | 1041 | type(vector_field), intent(in) :: Flux,U_nl | ||
1207 | 1042 | type(vector_field), intent(inout) :: QF | ||
1208 | 1043 | ! | ||
1209 | 1044 | type(csr_sparsity), pointer :: Q_sparsity | ||
1210 | 1045 | type(csr_matrix) :: Adv_mat | ||
1211 | 1046 | type(scalar_field) :: Q_rhs, Qtest1,Qtest2 | ||
1212 | 1047 | type(vector_field), pointer :: X, down | ||
1213 | 1048 | type(scalar_field), pointer :: Q_old | ||
1214 | 1049 | integer :: ele | ||
1215 | 1050 | real :: dt, t_theta, residual | ||
1216 | 1051 | |||
1217 | 1052 | ewrite(1,*) ' subroutine solve_advection_cg_tracer(' | ||
1218 | 1053 | |||
1219 | 1054 | X=>extract_vector_field(state, "Coordinate") | ||
1220 | 1055 | down=>extract_vector_field(state, "GravityDirection") | ||
1221 | 1056 | Q_old=>extract_scalar_field(state, "Old"//trim(Q%name)) | ||
1222 | 1057 | call get_option('/timestepping/timestep',dt) | ||
1223 | 1058 | call get_option('/timestepping/theta',t_theta) | ||
1224 | 1059 | |||
1225 | 1060 | !set up matrix and rhs | ||
1226 | 1061 | Q_sparsity => get_csr_sparsity_firstorder(state, Q%mesh, Q%mesh) | ||
1227 | 1062 | call allocate(adv_mat, Q_sparsity) ! Add data space to the sparsity | ||
1228 | 1063 | ! pattern. | ||
1229 | 1064 | call zero(adv_mat) | ||
1230 | 1065 | call allocate(Q_rhs,Q%mesh,trim(Q%name)//"RHS") | ||
1231 | 1066 | call zero(Q_rhs) | ||
1232 | 1067 | call zero(QF) | ||
1233 | 1068 | |||
1234 | 1069 | do ele = 1, ele_count(Q) | ||
1235 | 1070 | call construct_advection_cg_tracer_ele(Q_rhs,adv_mat,Q,D,D_old,Flux& | ||
1236 | 1071 | &,X,down,U_nl,dt,t_theta,ele) | ||
1237 | 1072 | end do | ||
1238 | 1073 | |||
1239 | 1074 | ewrite(2,*) 'Q_RHS', maxval(abs(Q_rhs%val)) | ||
1240 | 1075 | |||
1241 | 1076 | call petsc_solve(Q,adv_mat,Q_rhs) | ||
1242 | 1077 | |||
1243 | 1078 | !! Compute the PV flux to pass to velocity equation | ||
1244 | 1079 | do ele = 1, ele_count(Q) | ||
1245 | 1080 | call construct_pv_flux_ele(QF,Q,Q_old,D,D_old,Flux,& | ||
1246 | 1081 | X,down,U_nl,t_theta,ele) | ||
1247 | 1082 | end do | ||
1248 | 1083 | |||
1249 | 1084 | if(have_option('/material_phase::Fluid/scalar_field::PotentialVorticity/& | ||
1250 | 1085 | &prognostic/debug')) then | ||
1251 | 1086 | call allocate(Qtest1,Q%mesh,trim(Q%name)//"test1") | ||
1252 | 1087 | call zero(Qtest1) | ||
1253 | 1088 | call allocate(Qtest2,Q%mesh,trim(Q%name)//"test2") | ||
1254 | 1089 | call zero(Qtest2) | ||
1255 | 1090 | do ele = 1, ele_count(Q) | ||
1256 | 1091 | call test_pv_flux_ele(Qtest1,Qtest2,QF,Q,Q_old,D,D_old,& | ||
1257 | 1092 | Flux,X,down,t_theta,ele) | ||
1258 | 1093 | end do | ||
1259 | 1094 | ewrite(2,*) 'Error = ', maxval(abs(Qtest2%val)), maxval(abs(Qtest1%val)) | ||
1260 | 1095 | ewrite(2,*) 'Error = ', maxval(abs(Qtest1%val-Qtest2%val)), maxval(Qtest1%val) | ||
1261 | 1096 | residual = maxval(abs(Qtest1%val-Qtest2%val))/max(1.0& | ||
1262 | 1097 | &,maxval(abs(Qtest1%val))) | ||
1263 | 1098 | ewrite(2,*) 'residual from qtest', residual | ||
1264 | 1099 | assert(residual<1.0e-6) | ||
1265 | 1100 | ewrite(2,*) 'test passed' | ||
1266 | 1101 | call deallocate(Qtest1) | ||
1267 | 1102 | call deallocate(Qtest2) | ||
1268 | 1103 | end if | ||
1269 | 1104 | |||
1270 | 1105 | call deallocate(adv_mat) | ||
1271 | 1106 | call deallocate(Q_rhs) | ||
1272 | 1107 | |||
1273 | 1108 | ewrite(1,*) 'END subroutine solve_advection_cg_tracer(' | ||
1274 | 1109 | |||
1275 | 1110 | end subroutine solve_advection_cg_tracer | ||
1276 | 1111 | |||
1277 | 1112 | subroutine construct_advection_cg_tracer_ele(Q_rhs,adv_mat,Q,D,D_old,Flux,& | ||
1278 | 1113 | & X,down,U_nl,dt,t_theta,ele) | ||
1279 | 1114 | type(scalar_field), intent(in) :: D,D_old,Q | ||
1280 | 1115 | type(scalar_field), intent(inout) :: Q_rhs | ||
1281 | 1116 | type(csr_matrix), intent(inout) :: Adv_mat | ||
1282 | 1117 | type(vector_field), intent(in) :: X, Flux, down,U_nl | ||
1283 | 1118 | integer, intent(in) :: ele | ||
1284 | 1119 | real, intent(in) :: dt, t_theta | ||
1285 | 1120 | ! | ||
1286 | 1121 | real, dimension(ele_loc(Q,ele),ele_loc(Q,ele)) :: l_adv_mat,tmp_mat | ||
1287 | 1122 | real, dimension(ele_loc(Q,ele)) :: l_rhs | ||
1288 | 1123 | real, dimension(U_nl%dim,ele_ngi(U_nl,ele)) :: U_nl_gi | ||
1289 | 1124 | real, dimension(X%dim, ele_ngi(U_nl,ele)) :: U_cart_gi | ||
1290 | 1125 | real, dimension(Flux%dim,ele_ngi(Flux,ele)) :: Flux_gi | ||
1291 | 1126 | real, dimension(Flux%dim,FLux%dim,ele_ngi(Flux,ele)) :: Grad_Flux_gi,& | ||
1292 | 1127 | & tensor_gi | ||
1293 | 1128 | real, dimension(ele_ngi(X,ele)) :: detwei, Q_gi, D_gi, & | ||
1294 | 1129 | & D_old_gi,div_flux_gi, detwei_l, detJ | ||
1295 | 1130 | real, dimension(ele_loc(D,ele)) :: D_val | ||
1296 | 1131 | real, dimension(ele_loc(Q,ele)) :: Q_val | ||
1297 | 1132 | real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J | ||
1298 | 1133 | type(element_type), pointer :: Q_shape, Flux_shape | ||
1299 | 1134 | integer :: loc,dim1,dim2,gi | ||
1300 | 1135 | real :: tau, alpha, tol, area, h,c_sc,disc_indicator | ||
1301 | 1136 | real, dimension(ele_ngi(X,ele),ele_ngi(X,ele)) :: Metric, MetricT | ||
1302 | 1137 | real, dimension(mesh_dim(X),ele_ngi(X,ele)) :: gradQ | ||
1303 | 1138 | type(element_type), pointer :: D_shape | ||
1304 | 1139 | ! | ||
1305 | 1140 | |||
1306 | 1141 | ! Get J and detwei | ||
1307 | 1142 | call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), detwei& | ||
1308 | 1143 | &=detwei, J=J, detJ=detJ) | ||
1309 | 1144 | |||
1310 | 1145 | D_shape => ele_shape(D,ele) | ||
1311 | 1146 | Q_shape => ele_shape(Q,ele) | ||
1312 | 1147 | Flux_shape => ele_shape(Flux,ele) | ||
1313 | 1148 | D_val = invert_pi_ele(ele_val(D,ele),D_shape,detwei) | ||
1314 | 1149 | D_gi = matmul(transpose(D_shape%n),D_val) | ||
1315 | 1150 | D_val = invert_pi_ele(ele_val(D_old,ele),D_shape,detwei) | ||
1316 | 1151 | D_old_gi = matmul(transpose(D_shape%n),D_val) | ||
1317 | 1152 | Q_gi = ele_val_at_quad(Q,ele) | ||
1318 | 1153 | Flux_gi = ele_val_at_quad(Flux,ele) | ||
1319 | 1154 | U_nl_gi = ele_val_at_quad(U_nl,ele) | ||
1320 | 1155 | Q_val = ele_val(Q,ele) | ||
1321 | 1156 | |||
1322 | 1157 | !Equations | ||
1323 | 1158 | ! D^{n+1} = D^n + \pi div Flux | ||
1324 | 1159 | ! If define \int_K\phi\hat{D}/J J dS = \int_K\phi D J dS, then | ||
1325 | 1160 | ! <\phi,\pi(\hat{D}/J)> = <\phi, D> for all \phi, | ||
1326 | 1161 | ! => \pi \hat{D} = D | ||
1327 | 1162 | ! \hat{D}^{n+1}/J = \hat{D}^n/J + div Flux | ||
1328 | 1163 | |||
1329 | 1164 | ! So (integrating by parts) | ||
1330 | 1165 | ! <\gamma, \hat{D}^{n+1}/J> = | ||
1331 | 1166 | ! <\gamma, D^n/J> - <\nabla\gamma, F> | ||
1332 | 1167 | |||
1333 | 1168 | ! and a consistent theta-method for Q is | ||
1334 | 1169 | ! <\gamma,\hat{D}^{n+1}Q^{n+1}/J> = | ||
1335 | 1170 | ! <\gamma, \hat{D}^nQ^n/J> - | ||
1336 | 1171 | ! \theta<\nabla\gamma,FQ^{n+1}> | ||
1337 | 1172 | ! -(1-\theta)<\nabla\gamma,FQ^n> | ||
1338 | 1173 | |||
1339 | 1174 | detwei_l = Q_shape%quadrature%weight | ||
1340 | 1175 | !Advection terms | ||
1341 | 1176 | !! <-grad gamma, FQ > | ||
1342 | 1177 | !! Can be evaluated locally using pullback | ||
1343 | 1178 | tmp_mat = dshape_dot_vector_shape(& | ||
1344 | 1179 | Q_shape%dn,Flux_gi,Q_shape,detwei_l) | ||
1345 | 1180 | l_adv_mat = t_theta*tmp_mat | ||
1346 | 1181 | l_rhs = -(1-t_theta)*matmul(tmp_mat,Q_val) | ||
1347 | 1182 | |||
1348 | 1183 | !Mass terms | ||
1349 | 1184 | !! <gamma, QD> | ||
1350 | 1185 | !! Requires transformed integral | ||
1351 | 1186 | tmp_mat = shape_shape(ele_shape(Q,ele),ele_shape(Q,ele),D_gi*detwei_l) | ||
1352 | 1187 | l_adv_mat = l_adv_mat + tmp_mat | ||
1353 | 1188 | tmp_mat = shape_shape(ele_shape(Q,ele),ele_shape(Q,ele),D_old_gi*detwei_l) | ||
1354 | 1189 | l_rhs = l_rhs + matmul(tmp_mat,Q_val) | ||
1355 | 1190 | |||
1356 | 1191 | !Laplacian filter options | ||
1357 | 1192 | if(have_option(trim(Q%option_path)//'/prognostic/spatial_discretisation/& | ||
1358 | 1193 | &continuous_galerkin/laplacian_filter')) then | ||
1359 | 1194 | FLExit('Laplacian filter not coded yet.') | ||
1360 | 1195 | end if | ||
1361 | 1196 | !Streamline upwinding options | ||
1362 | 1197 | if(have_option(trim(Q%option_path)//'/prognostic/spatial_discretisation/& | ||
1363 | 1198 | &continuous_galerkin/streamline_upwinding')) then | ||
1364 | 1199 | !! Replace test function \gamma with \gamma + \tau u\cdot\nabla\gamma | ||
1365 | 1200 | call get_option(trim(Q%option_path)//'/prognostic/spatial_discretisat& | ||
1366 | 1201 | &ion/continuous_galerkin/streamline_upwinding/alpha',alpha) | ||
1367 | 1202 | call get_option(trim(Q%option_path)//'/prognostic/spatial_discretisat& | ||
1368 | 1203 | &ion/continuous_galerkin/streamline_upwinding/tol',tol) | ||
1369 | 1204 | |||
1370 | 1205 | !Gradients of fluxes | ||
1371 | 1206 | Grad_flux_gi = ele_grad_at_quad(Flux,ele,Flux_shape%dn) | ||
1372 | 1207 | div_flux_gi = 0. | ||
1373 | 1208 | tensor_gi = 0. | ||
1374 | 1209 | do dim1 = 1, Flux%dim | ||
1375 | 1210 | div_flux_gi = div_flux_gi + grad_flux_gi(dim1,dim1,:) | ||
1376 | 1211 | do dim2 = 1, Flux%dim | ||
1377 | 1212 | tensor_gi(dim1,dim2,:) = U_nl_gi(dim1,:)*Flux_gi(dim2,:) | ||
1378 | 1213 | end do | ||
1379 | 1214 | end do | ||
1380 | 1215 | !real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J | ||
1381 | 1216 | do gi = 1, ele_ngi(X,ele) | ||
1382 | 1217 | U_cart_gi(:,gi) = matmul(transpose(J(:,:,gi)),U_nl_gi(:,gi))/detJ(gi) | ||
1383 | 1218 | end do | ||
1384 | 1219 | |||
1385 | 1220 | if(maxval(sqrt(sum(U_cart_gi**2,1)))<tol) then | ||
1386 | 1221 | tau = 0. | ||
1387 | 1222 | else | ||
1388 | 1223 | area = sum(detwei) | ||
1389 | 1224 | h = sqrt(4*area/sqrt(3.)) | ||
1390 | 1225 | tau = h*alpha/maxval(sqrt(sum(U_cart_gi**2,1))) | ||
1391 | 1226 | end if | ||
1392 | 1227 | !! Mass terms | ||
1393 | 1228 | !! tau < u . grad gamma, QD> | ||
1394 | 1229 | !! = tau < grad gamma, u QD> | ||
1395 | 1230 | !! can do locally because of pullback formula | ||
1396 | 1231 | tmp_mat = dshape_dot_vector_shape(& | ||
1397 | 1232 | Q_shape%dn,U_nl_gi,Q_shape,D_gi*detwei_l/detJ) | ||
1398 | 1233 | l_adv_mat = l_adv_mat + tau*tmp_mat | ||
1399 | 1234 | tmp_mat = dshape_dot_vector_shape(& | ||
1400 | 1235 | Q_shape%dn,U_nl_gi,Q_shape,D_old_gi*detwei_l/detJ) | ||
1401 | 1236 | l_rhs = l_rhs + tau*matmul(tmp_mat,Q_val) | ||
1402 | 1237 | !! Advection terms | ||
1403 | 1238 | !! tau < u. grad gamma, div (F Q)> | ||
1404 | 1239 | !! = tau <grad gamma, u (F.grad Q + Q div F)> | ||
1405 | 1240 | !! can do locally because of pullback formula | ||
1406 | 1241 | !! tensor = u outer F (F already contains factor of dt) | ||
1407 | 1242 | tmp_mat = dshape_dot_vector_shape(& | ||
1408 | 1243 | Q_shape%dn,U_nl_gi,Q_shape,div_flux_gi*detwei_l/detJ) & | ||
1409 | 1244 | + dshape_tensor_dshape(& | ||
1410 | 1245 | Q_shape%dn,tensor_gi,Q_shape%dn,detwei_l/detJ) | ||
1411 | 1246 | l_adv_mat = l_adv_mat - tau*t_theta*tmp_mat | ||
1412 | 1247 | l_rhs = l_rhs + tau*(1-t_theta)*matmul(tmp_mat,Q_val) | ||
1413 | 1248 | end if | ||
1414 | 1249 | !Discontinuity capturing options | ||
1415 | 1250 | if(have_option(trim(Q%option_path)//'/prognostic/spatial_discretisation/& | ||
1416 | 1251 | &continuous_galerkin/discontinuity_capturing')) then | ||
1417 | 1252 | ! call have_option(& | ||
1418 | 1253 | ! &trim(Q%option_path)//'/prognostic/spatial_discretisation/& | ||
1419 | 1254 | ! &continuous_galerkin/discontinuity_capturing/& | ||
1420 | 1255 | ! &scaling_coefficient',c_sc) | ||
1421 | 1256 | |||
1422 | 1257 | ! do gi=1,ele_ngi(Q,ele) | ||
1423 | 1258 | ! Metric(:,:,gi)=matmul(J(:,:,gi), transpose(J(:,:,gi)))/detJ(gi) | ||
1424 | 1259 | ! end do | ||
1425 | 1260 | ! MetricT(1,1,:) = MetricT(2,2,:) | ||
1426 | 1261 | ! MetricT(2,2,:) = MetricT(1,1,:) | ||
1427 | 1262 | ! MetricT(1,2,:) =-MetricT(2,1,:) | ||
1428 | 1263 | ! MetricT(2,1,:) =-MetricT(1,2,:) | ||
1429 | 1264 | |||
1430 | 1265 | ! area = sum(detwei) | ||
1431 | 1266 | ! D_gi = ele_val_at_quad(D_old,ele) | ||
1432 | 1267 | ! tmp_mat = dshape_tensor_dshape(Q_shape%dn, & | ||
1433 | 1268 | ! area*D_gi*MetricT, Q_shape%dn,& | ||
1434 | 1269 | ! Q_shape%quadrature%weight) | ||
1435 | 1270 | ! l_rhs = l_rhs - dt*c_sc*matmul(tmp_mat,Q_val) | ||
1436 | 1271 | FLAbort('not ready yet') | ||
1437 | 1272 | end if | ||
1438 | 1273 | |||
1439 | 1274 | call addto(Q_rhs,ele_nodes(Q_rhs,ele),l_rhs) | ||
1440 | 1275 | call addto(adv_mat,ele_nodes(Q_rhs,ele),ele_nodes(Q_rhs,ele),& | ||
1441 | 1276 | l_adv_mat) | ||
1442 | 1277 | |||
1443 | 1278 | end subroutine construct_advection_cg_tracer_ele | ||
1444 | 1279 | |||
1445 | 1280 | subroutine construct_pv_flux_ele(QFlux,Q,Q_old,D,D_old,& | ||
1446 | 1281 | Flux,X,down,U_nl,t_theta,ele) | ||
1447 | 1282 | type(scalar_field), intent(in) :: Q,Q_old,D,D_old | ||
1448 | 1283 | type(vector_field), intent(inout) :: QFlux | ||
1449 | 1284 | type(vector_field), intent(in) :: X, Flux, Down, U_nl | ||
1450 | 1285 | integer, intent(in) :: ele | ||
1451 | 1286 | real, intent(in) :: t_theta | ||
1452 | 1287 | ! | ||
1453 | 1288 | real, dimension(mesh_dim(X),ele_loc(QFlux,ele)) :: QFlux_perp_rhs | ||
1454 | 1289 | real, dimension(Flux%dim,ele_ngi(Flux,ele)) :: Flux_gi, Flux_perp_gi,& | ||
1455 | 1290 | QFlux_gi | ||
1456 | 1291 | real, dimension(ele_ngi(X,ele)) :: Q_gi,Q_old_gi,D_gi,D_old_gi,Qbar_gi | ||
1457 | 1292 | real, dimension(ele_loc(D,ele)) :: D_val | ||
1458 | 1293 | real, dimension(X%dim, ele_ngi(Q, ele)) :: up_gi | ||
1459 | 1294 | real, dimension(mesh_dim(QFlux),mesh_dim(QFlux),ele_loc(QFlux,ele)& | ||
1460 | 1295 | &,ele_loc(QFlux,ele)) :: l_u_mat | ||
1461 | 1296 | real, dimension(mesh_dim(QFlux)*ele_loc(QFlux,ele),mesh_dim(QFlux)& | ||
1462 | 1297 | &*ele_loc(QFlux,ele)) :: solve_mat | ||
1463 | 1298 | real, dimension(ele_loc(Q,ele)) :: Q_test1_rhs, Q_test2_rhs | ||
1464 | 1299 | real, dimension(mesh_dim(Qflux)*ele_loc(QFlux,ele)) :: solve_rhs | ||
1465 | 1300 | real, dimension(mesh_dim(Qflux),ele_ngi(Qflux,ele)) :: & | ||
1466 | 1301 | contravariant_velocity_gi, velocity_perp_gi | ||
1467 | 1302 | type(element_type), pointer :: Q_shape, QFlux_shape, Flux_shape | ||
1468 | 1303 | integer :: loc, orientation,gi, dim1, dim2 | ||
1469 | 1304 | real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J | ||
1470 | 1305 | real, dimension(ele_ngi(X,ele)) :: detwei, detJ, detwei_l, div_flux_gi | ||
1471 | 1306 | real, dimension(mesh_dim(Q),ele_ngi(Q,ele)) :: gradQ_gi | ||
1472 | 1307 | real, dimension(Flux%dim,FLux%dim,ele_ngi(Flux,ele)) :: Grad_Flux_gi | ||
1473 | 1308 | real, dimension(U_nl%dim,ele_ngi(U_nl,ele)) :: U_nl_gi | ||
1474 | 1309 | real, dimension(X%dim,ele_ngi(U_nl,ele)) :: U_cart_gi | ||
1475 | 1310 | real :: residual, alpha, tol, tau, area, h | ||
1476 | 1311 | type(element_type), pointer :: D_shape | ||
1477 | 1312 | ! | ||
1478 | 1313 | |||
1479 | 1314 | D_shape => ele_shape(D,ele) | ||
1480 | 1315 | |||
1481 | 1316 | call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), J, detwei, detJ) | ||
1482 | 1317 | |||
1483 | 1318 | Q_shape => ele_shape(Q,ele) | ||
1484 | 1319 | Flux_shape => ele_shape(Flux,ele) | ||
1485 | 1320 | QFlux_shape => ele_shape(QFlux,ele) | ||
1486 | 1321 | Q_gi = ele_val_at_quad(Q,ele) | ||
1487 | 1322 | Q_old_gi = ele_val_at_quad(Q_old,ele) | ||
1488 | 1323 | Qbar_gi = t_theta*Q_gi + (1-t_theta)*Q_old_gi | ||
1489 | 1324 | D_val = invert_pi_ele(ele_val(D,ele),D_shape,detwei) | ||
1490 | 1325 | D_gi = matmul(transpose(D_shape%n),D_val) | ||
1491 | 1326 | D_val = invert_pi_ele(ele_val(D_old,ele),D_shape,detwei) | ||
1492 | 1327 | D_old_gi = matmul(transpose(D_shape%n),D_val) | ||
1493 | 1328 | Flux_gi = ele_val_at_quad(Flux,ele) | ||
1494 | 1329 | U_nl_gi = ele_val_at_quad(U_nl,ele) | ||
1495 | 1330 | |||
1496 | 1331 | detwei_l = Q_shape%quadrature%weight | ||
1497 | 1332 | up_gi = -ele_val_at_quad(down,ele) | ||
1498 | 1333 | |||
1499 | 1334 | call get_up_gi(X,ele,up_gi,orientation) | ||
1500 | 1335 | |||
1501 | 1336 | !Equations | ||
1502 | 1337 | ! <\nabla \gamma, (-v,u)> = <(\gamma_x, \gamma_y),(-v,u)> | ||
1503 | 1338 | ! = <(\gamma_y,-\gamma_x),( u,v)> | ||
1504 | 1339 | ! = <-\nabla^\perp\gamma,(u,v)> | ||
1505 | 1340 | |||
1506 | 1341 | ! and a consistent theta-method for Q is | ||
1507 | 1342 | ! <\gamma, D^{n+1}Q^{n+1}> = <\gamma, D^nQ^n> - | ||
1508 | 1343 | ! \theta<\nabla\gamma,FQ^{n+1}> | ||
1509 | 1344 | ! -(1-\theta)<\nabla\gamma,FQ^n> | ||
1510 | 1345 | ! So vorticity update is | ||
1511 | 1346 | ! <\gamma, \zeta^{n+1}> = <-\nabla^\perp\gamma,u^{n+1}> | ||
1512 | 1347 | != <-\nabla^\perp\gamma,u^n> | ||
1513 | 1348 | ! +\theta<-\nabla^\perp\gamma,(FQ^{n+1})^\perp> | ||
1514 | 1349 | ! +(1-\theta)<-\nabla^\perp\gamma,(FQ^n)^\perp> | ||
1515 | 1350 | ! taking w = -\nabla^\perp\gamma, | ||
1516 | 1351 | ! <w,u^{n+1}>=<w,u^n> + <w,F(\theta Q^{n+1}+(1-\theta)Q^n)^\perp> | ||
1517 | 1352 | ! <w,u^{n+1}>=<w,u^n> + <w,\bar{Q}^\perp> | ||
1518 | 1353 | ! We actually store the inner product with test function (FQ_rhs) | ||
1519 | 1354 | ! (avoids having to do a solve) | ||
1520 | 1355 | |||
1521 | 1356 | do gi = 1, ele_ngi(QFlux,ele) | ||
1522 | 1357 | QFlux_gi(:,gi) = Flux_gi(:,gi)*Qbar_gi(gi) | ||
1523 | 1358 | end do | ||
1524 | 1359 | |||
1525 | 1360 | !! Additional dissipative terms. | ||
1526 | 1361 | !! They all take the form | ||
1527 | 1362 | !! <\nabla gamma, G> | ||
1528 | 1363 | !! Which can be evaluated in local coordinates due to pullback magic | ||
1529 | 1364 | |||
1530 | 1365 | !Laplacian filter options | ||
1531 | 1366 | if(have_option(trim(Q%option_path)//'/prognostic/spatial_discretisation/& | ||
1532 | 1367 | &continuous_galerkin/laplacian_filter')) then | ||
1533 | 1368 | FLExit('Laplacian filter not coded yet.') | ||
1534 | 1369 | end if | ||
1535 | 1370 | !Streamline upwinding options | ||
1536 | 1371 | if(have_option(trim(Q%option_path)//'/prognostic/spatial_discretisation/& | ||
1537 | 1372 | &continuous_galerkin/streamline_upwinding')) then | ||
1538 | 1373 | !! Replace test function \gamma with \gamma + \tau u\cdot\nabla\gamma | ||
1539 | 1374 | call get_option(trim(Q%option_path)//'/prognostic/spatial_discretisat& | ||
1540 | 1375 | &ion/continuous_galerkin/streamline_upwinding/alpha',alpha) | ||
1541 | 1376 | call get_option(trim(Q%option_path)//'/prognostic/spatial_discretisat& | ||
1542 | 1377 | &ion/continuous_galerkin/streamline_upwinding/tol',tol) | ||
1543 | 1378 | |||
1544 | 1379 | !Gradients of fluxes | ||
1545 | 1380 | Grad_flux_gi = ele_grad_at_quad(Flux,ele,Flux_shape%dn) | ||
1546 | 1381 | div_flux_gi = 0. | ||
1547 | 1382 | do dim1 = 1, Flux%dim | ||
1548 | 1383 | div_flux_gi = div_flux_gi + grad_flux_gi(dim1,dim1,:) | ||
1549 | 1384 | end do | ||
1550 | 1385 | !real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J | ||
1551 | 1386 | do gi = 1, ele_ngi(X,ele) | ||
1552 | 1387 | U_cart_gi(:,gi) = matmul(transpose(J(:,:,gi)),U_nl_gi(:,gi))/detJ(gi) | ||
1553 | 1388 | end do | ||
1554 | 1389 | |||
1555 | 1390 | if(maxval(sqrt(sum(U_cart_gi**2,1)))<tol) then | ||
1556 | 1391 | tau = 0. | ||
1557 | 1392 | else | ||
1558 | 1393 | area = sum(detwei) | ||
1559 | 1394 | h = sqrt(4*area/sqrt(3.)) | ||
1560 | 1395 | tau = h*alpha/maxval(sqrt(sum(U_cart_gi**2,1))) | ||
1561 | 1396 | end if | ||
1562 | 1397 | !! Mass terms | ||
1563 | 1398 | !! tau < u . grad gamma, -(QD)^{n+1}+(QD)^n> | ||
1564 | 1399 | !! = tau < grad gamma, u(-(QD)^{n+1}+(QD)^n)> | ||
1565 | 1400 | !! can do locally because of pullback formula | ||
1566 | 1401 | !! Extra flux is tau u\Delta(QD) | ||
1567 | 1402 | !! minus sign from definition of vorticity | ||
1568 | 1403 | !! tau | ||
1569 | 1404 | do gi = 1, ele_ngi(Q,ele) | ||
1570 | 1405 | QFLux_gi(:,gi) = QFlux_gi(:,gi) + & | ||
1571 | 1406 | & tau*U_nl_gi(:,gi)*(Q_gi(gi)*D_gi(gi)-Q_old_gi(gi)& | ||
1572 | 1407 | &*D_old_gi(gi))/detJ(gi) | ||
1573 | 1408 | |||
1574 | 1409 | end do | ||
1575 | 1410 | !! Advection terms | ||
1576 | 1411 | !! tau < u. grad gamma, div (F\bar{Q})> | ||
1577 | 1412 | !! can do locally because of pullback formula | ||
1578 | 1413 | !! = <grad gamma, tau u (F.grad\bar{Q} + \bar{Q} div F)/det J>_{Ehat} | ||
1579 | 1414 | !! minus sign because of definition of vorticity | ||
1580 | 1415 | |||
1581 | 1416 | GradQ_gi = t_theta*ele_grad_at_quad(Q,ele,Q_shape%dn) + & | ||
1582 | 1417 | (1-t_theta)*ele_grad_at_quad(Q_old,ele,Q_shape%dn) | ||
1583 | 1418 | |||
1584 | 1419 | do gi = 1, ele_ngi(Q,ele) | ||
1585 | 1420 | QFlux_gi(:,gi) = QFlux_gi(:,gi) - tau*u_nl_gi(:,gi)*(& | ||
1586 | 1421 | & sum(Flux_gi(:,gi)*gradQ_gi(:,gi)) + & | ||
1587 | 1422 | & div_flux_gi(gi)*Qbar_gi(gi))/detJ(gi) | ||
1588 | 1423 | end do | ||
1589 | 1424 | end if | ||
1590 | 1425 | !Discontinuity capturing options | ||
1591 | 1426 | if(have_option(trim(Q%option_path)//'/prognostic/spatial_discretisation/& | ||
1592 | 1427 | &continuous_galerkin/discontinuity_capturing')) then | ||
1593 | 1428 | FLExit('Discontinuity capturing not coded yet.') | ||
1594 | 1429 | end if | ||
1595 | 1430 | |||
1596 | 1431 | !! Evaluate | ||
1597 | 1432 | !! < w, F^\perp > in local coordinates | ||
1598 | 1433 | |||
1599 | 1434 | Flux_perp_gi(1,:) = -orientation*QFlux_gi(2,:) | ||
1600 | 1435 | Flux_perp_gi(2,:) = orientation*QFlux_gi(1,:) | ||
1601 | 1436 | |||
1602 | 1437 | QFlux_perp_rhs = shape_vector_rhs(QFlux_shape,Flux_perp_gi,& | ||
1603 | 1438 | & QFlux_shape%quadrature%weight) | ||
1604 | 1439 | |||
1605 | 1440 | call set(QFlux,ele_nodes(QFlux,ele),QFlux_perp_rhs) | ||
1606 | 1441 | |||
1607 | 1442 | end subroutine construct_pv_flux_ele | ||
1608 | 1443 | |||
1609 | 1444 | subroutine test_pv_flux_ele(Qtest1,Qtest2,QFlux,Q,Q_old,D,D_old,& | ||
1610 | 1445 | Flux,X,down,t_theta,ele) | ||
1611 | 1446 | type(scalar_field), intent(in) :: Q,Q_old,D,D_old | ||
1612 | 1447 | type(scalar_field), intent(inout) :: Qtest1, Qtest2 | ||
1613 | 1448 | type(vector_field), intent(in) :: X, Flux, Down,QFlux | ||
1614 | 1449 | integer, intent(in) :: ele | ||
1615 | 1450 | real, intent(in) :: t_theta | ||
1616 | 1451 | ! | ||
1617 | 1452 | real, dimension(mesh_dim(X),ele_loc(QFlux,ele)) :: QFlux_perp_rhs | ||
1618 | 1453 | real, dimension(ele_ngi(X,ele)) :: Q_gi,Q_old_gi,D_gi,D_old_gi | ||
1619 | 1454 | real, dimension(Flux%dim,ele_ngi(Flux,ele)) :: Flux_gi, Flux_perp_gi | ||
1620 | 1455 | real, dimension(X%dim, ele_ngi(Q, ele)) :: up_gi | ||
1621 | 1456 | real, dimension(mesh_dim(QFlux), mesh_dim(QFlux), ele_ngi(QFlux,ele)) & | ||
1622 | 1457 | :: Metric, Metricf | ||
1623 | 1458 | real, dimension(mesh_dim(QFlux),mesh_dim(QFlux),ele_loc(QFlux,ele)& | ||
1624 | 1459 | &,ele_loc(QFlux,ele)) :: l_u_mat | ||
1625 | 1460 | real, dimension(mesh_dim(QFlux)*ele_loc(QFlux,ele),mesh_dim(QFlux)& | ||
1626 | 1461 | &*ele_loc(QFlux,ele)) :: solve_mat | ||
1627 | 1462 | real, dimension(ele_loc(Q,ele)) :: Q_test1_rhs, Q_test2_rhs | ||
1628 | 1463 | real, dimension(mesh_dim(Qflux)*ele_loc(QFlux,ele)) :: solve_rhs | ||
1629 | 1464 | real, dimension(mesh_dim(Qflux),ele_ngi(Qflux,ele)) :: & | ||
1630 | 1465 | contravariant_velocity_gi, velocity_perp_gi | ||
1631 | 1466 | type(element_type), pointer :: Q_shape, QFlux_shape,D_shape | ||
1632 | 1467 | integer :: loc, orientation,gi, dim1, dim2 | ||
1633 | 1468 | real, dimension(mesh_dim(X), X%dim, ele_ngi(X,ele)) :: J | ||
1634 | 1469 | real, dimension(ele_ngi(X,ele)) :: detwei, detJ | ||
1635 | 1470 | real, dimension(X%dim, X%dim, ele_ngi(X,ele)) :: rot | ||
1636 | 1471 | real, dimension(ele_loc(D,ele)) :: D_val | ||
1637 | 1472 | |||
1638 | 1473 | D_shape => ele_shape(D,ele) | ||
1639 | 1474 | Q_shape => ele_shape(Q,ele) | ||
1640 | 1475 | QFlux_shape => ele_shape(QFlux,ele) | ||
1641 | 1476 | Q_gi = ele_val_at_quad(Q,ele) | ||
1642 | 1477 | Q_old_gi = ele_val_at_quad(Q_old,ele) | ||
1643 | 1478 | |||
1644 | 1479 | call compute_jacobian(ele_val(X,ele), ele_shape(X,ele), J, detwei, detJ) | ||
1645 | 1480 | D_val = invert_pi_ele(ele_val(D,ele),D_shape,detwei) | ||
1646 | 1481 | D_gi = matmul(transpose(D_shape%n),D_val) | ||
1647 | 1482 | D_val = invert_pi_ele(ele_val(D_old,ele),D_shape,detwei) | ||
1648 | 1483 | D_old_gi = matmul(transpose(D_shape%n),D_val) | ||
1649 | 1484 | |||
1650 | 1485 | Flux_gi = ele_val_at_quad(Flux,ele) | ||
1651 | 1486 | Qflux_perp_rhs = ele_val(Qflux,ele) | ||
1652 | 1487 | |||
1653 | 1488 | do gi = 1, ele_ngi(QFlux,ele) | ||
1654 | 1489 | Flux_gi(:,gi) = Flux_gi(:,gi)*(t_theta*Q_gi(gi) + & | ||
1655 | 1490 | (1-t_theta)*Q_old_gi(gi)) | ||
1656 | 1491 | end do | ||
1657 | 1492 | |||
1658 | 1493 | up_gi = -ele_val_at_quad(down,ele) | ||
1659 | 1494 | call get_up_gi(X,ele,up_gi,orientation) | ||
1660 | 1495 | |||
1661 | 1496 | do gi=1, ele_ngi(QFlux,ele) | ||
1662 | 1497 | rot(1,:,gi)=(/0.,-up_gi(3,gi),up_gi(2,gi)/) | ||
1663 | 1498 | rot(2,:,gi)=(/up_gi(3,gi),0.,-up_gi(1,gi)/) | ||
1664 | 1499 | rot(3,:,gi)=(/-up_gi(2,gi),up_gi(1,gi),0./) | ||
1665 | 1500 | end do | ||
1666 | 1501 | do gi=1,ele_ngi(QFlux,ele) | ||
1667 | 1502 | Metric(:,:,gi)=matmul(J(:,:,gi), transpose(J(:,:,gi)))/detJ(gi) | ||
1668 | 1503 | Metricf(:,:,gi)=matmul(J(:,:,gi), & | ||
1669 | 1504 | matmul(rot(:,:,gi), transpose(J(:,:,gi))))/detJ(gi) | ||
1670 | 1505 | end do | ||
1671 | 1506 | |||
1672 | 1507 | l_u_mat = shape_shape_tensor(qflux_shape, qflux_shape, & | ||
1673 | 1508 | qflux_shape%quadrature%weight,Metric) | ||
1674 | 1509 | solve_mat = 0. | ||
1675 | 1510 | solve_rhs = 0. | ||
1676 | 1511 | do dim1 = 1, mesh_dim(qflux) | ||
1677 | 1512 | do dim2 = 1, mesh_dim(qflux) | ||
1678 | 1513 | solve_mat( (dim1-1)*ele_loc(qflux,ele)+1:dim1*ele_loc(qflux,ele),& | ||
1679 | 1514 | (dim2-1)*ele_loc(qflux,ele)+1:dim2*ele_loc(qflux,ele)) = & | ||
1680 | 1515 | l_u_mat(dim1,dim2,:,:) | ||
1681 | 1516 | end do | ||
1682 | 1517 | solve_rhs( (dim1-1)*ele_loc(qflux,ele)+1:dim1*ele_loc(qflux,ele))& | ||
1683 | 1518 | = QFlux_perp_rhs(dim1,:) | ||
1684 | 1519 | end do | ||
1685 | 1520 | call solve(solve_mat,solve_rhs) | ||
1686 | 1521 | !Don't need to include constraints since u eqn is | ||
1687 | 1522 | ! <w,Q^\perp> + <<[w],\lambda>> + \sum_i C_i\cdot w \Gamma_i = RHS | ||
1688 | 1523 | ! but if w=-\nabla^\perp\gamma then [w]=0 and C_i\cdot w=0. | ||
1689 | 1524 | |||
1690 | 1525 | do dim1 = 1, mesh_dim(qflux) | ||
1691 | 1526 | flux_perp_gi(dim1,:) = matmul(transpose(QFlux_shape%n),& | ||
1692 | 1527 | solve_rhs((dim1-1)*ele_loc(qflux,ele)+1:dim1*ele_loc(qflux& | ||
1693 | 1528 | &,ele))) | ||
1694 | 1529 | end do | ||
1695 | 1530 | !Now compute vorticity | ||
1696 | 1531 | do gi=1,ele_ngi(X,ele) | ||
1697 | 1532 | Metric(:,:,gi)=matmul(J(:,:,gi), transpose(J(:,:,gi)))/detJ(gi) | ||
1698 | 1533 | contravariant_velocity_gi(:,gi) = & | ||
1699 | 1534 | matmul(Metric(:,:,gi),Flux_perp_gi(:,gi)) | ||
1700 | 1535 | end do | ||
1701 | 1536 | |||
1702 | 1537 | !!Q_test1_rhs contains -<\nabla\gamma,Q> | ||
1703 | 1538 | |||
1704 | 1539 | select case(mesh_dim(X)) | ||
1705 | 1540 | case (2) | ||
1706 | 1541 | velocity_perp_gi(1,:) = -contravariant_velocity_gi(2,:) | ||
1707 | 1542 | velocity_perp_gi(2,:) = contravariant_velocity_gi(1,:) | ||
1708 | 1543 | Q_test1_rhs = dshape_dot_vector_rhs(Q%mesh%shape%dn, & | ||
1709 | 1544 | velocity_perp_gi,Q%mesh%shape%quadrature%weight) | ||
1710 | 1545 | Q_test1_rhs = q_test1_rhs*orientation | ||
1711 | 1546 | case default | ||
1712 | 1547 | FLAbort('Exterior derivative not implemented for given mesh dimension' |
Hi Colin,
Only had a quick glance but can I check if these changes provide any new functionality to fluidity itself? I ask because you've changed mesh_options.rn[cg] and these options are now appearing in the fluidity schema too but I can't see how they're applicable to fluidity as I can only trace them back to the hybridized helmholtz solver. If they're not used in fluidity then I don't think they should appear in a piece of the schema that is visible within an flml but I may be missing something.
Unrelated to this merge but the parent option (geometry/ mesh/from_ mesh/constrain_ type) for these new options is undocumented in the schema.
Cheers,
Cian