← Back to branch summary

~fluidity-core/fluidity/refactor-netcdf

« back to all changes in this revision

Viewing changes to femtools/tests/test_vector_lon_lat_height_2_cartesian.F90

  • Committer: Jon Hill
  • Date: 2013-02-16 09:01:40 UTC
  • mfrom: (3981.7.159 fluidity)
  • Revision ID: jon.hill@imperial.ac.uk-20130216090140-bplzxqzdk1eik4za
Megre from trunk, fixing several conflicts

Show diffs side-by-side

added added

removed removed

Lines of Context:
femtools/tests/test_vector_lon_lat_height_2_cartesian.F90
 
1
!    Copyright (C) 2012 Imperial College London and others.
 
2
!    
 
3
!    Please see the AUTHORS file in the main source directory for a full list
 
4
!    of copyright holders.
 
5
!
 
6
!    Prof. C Pain
 
7
!    Applied Modelling and Computation Group
 
8
!    Department of Earth Science and Engineering
 
9
!    Imperial College London
 
10
!
 
11
!    amcgsoftware@imperial.ac.uk
 
12
!    
 
13
!    This library is free software; you can redistribute it and/or
 
14
!    modify it under the terms of the GNU Lesser General Public
 
15
!    License as published by the Free Software Foundation,
 
16
!    version 2.1 of the License.
 
17
!
 
18
!    This library is distributed in the hope that it will be useful,
 
19
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
 
20
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
21
!    Lesser General Public License for more details.
 
22
!
 
23
!    You should have received a copy of the GNU Lesser General Public
 
24
!    License along with this library; if not, write to the Free Software
 
25
!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
 
26
!    USA
 
27
subroutine test_vector_lon_lat_height_2_cartesian
 
28
  !Subroutine/unit-test of correct vector basis change from a
 
29
  ! meridional-zonal-vertical system to a Cartesian system.
 
30
 
 
31
  use fields
 
32
  use vtk_interfaces
 
33
  use state_module
 
34
  use Coordinates
 
35
  use unittest_tools
 
36
  implicit none
 
37
 
 
38
  type(state_type) :: state
 
39
  type(mesh_type), pointer :: mesh
 
40
  type(vector_field), pointer :: CartesianCoordinate, LonLatHeightCoordinate
 
41
  type(vector_field), pointer :: UnitRadialVector_inCartesian
 
42
  type(vector_field), pointer :: UnitPolarVector_inCartesian
 
43
  type(vector_field), pointer :: UnitAzimuthalVector_inCartesian
 
44
  type(vector_field), pointer :: UnitRadialVector_inZonalMeridionalRadial
 
45
  type(vector_field), pointer :: UnitPolarVector_inZonalMeridionalRadial
 
46
  type(vector_field), pointer :: UnitAzimuthalVector_inZonalMeridionalRadial
 
47
  type(vector_field) :: radialVectorDifference, &
 
48
                        polarVectorDifference, &
 
49
                        azimuthalVectorDifference
 
50
  real, dimension(3) :: meridionalZonalVerticalVectorComponents, &
 
51
                        cartesianVectorComponents
 
52
  real, dimension(3) :: XYZ, LLH !Arrays containing a signle node's position vector
 
53
                                 ! components in Cartesian & lon-lat-height bases.
 
54
  integer :: node
 
55
  logical :: fail
 
56
 
 
57
  call vtk_read_state("data/on_sphere_rotations/spherical_shell_withFields.vtu", state)
 
58
  mesh => extract_mesh(state, "Mesh")
 
59
  CartesianCoordinate => extract_vector_field(state, "CartesianCoordinate")
 
60
  LonLatHeightCoordinate => extract_vector_field(state, "lonlatradius")
 
61
  UnitRadialVector_inCartesian => extract_vector_field(state, &
 
62
                                      "UnitRadialVector_inCartesian")
 
63
  UnitPolarVector_inCartesian => extract_vector_field(state, &
 
64
                                      "UnitPolarVector_inCartesian")
 
65
  UnitAzimuthalVector_inCartesian => extract_vector_field(state, &
 
66
                                      "UnitAzimuthalVector_inCartesian")
 
67
  UnitRadialVector_inZonalMeridionalRadial => extract_vector_field(state, &
 
68
                                      "UnitRadialVector_inZonalMeridionalRadial")
 
69
  UnitPolarVector_inZonalMeridionalRadial => extract_vector_field(state, &
 
70
                                      "UnitPolarVector_inZonalMeridionalRadial")
 
71
  UnitAzimuthalVector_inZonalMeridionalRadial => extract_vector_field(state, &
 
72
                                      "UnitAzimuthalVector_inZonalMeridionalRadial")
 
73
 
 
74
  !Test the change of basis from spherical-polar to Cartesian.
 
75
 
 
76
  call allocate(radialVectorDifference, 3 , mesh, 'radialVectorDifference')
 
77
  call allocate(polarVectorDifference, 3 , mesh, 'polarVectorDifference')
 
78
  call allocate(azimuthalVectorDifference, 3 , mesh, 'azimuthalVectorDifference')
 
79
 
 
80
  !Convert unit-radial vector components into to Cartesian basis. Then compare
 
81
  ! with vector already in Cartesian basis, obtained from vtu.
 
82
  do node=1,node_count(LonLatHeightCoordinate)
 
83
    LLH = node_val(LonLatHeightCoordinate, node)
 
84
    meridionalZonalVerticalVectorComponents = &
 
85
                   node_val(UnitRadialVector_inZonalMeridionalRadial, node)
 
86
    call vector_lon_lat_height_2_cartesian(meridionalZonalVerticalVectorComponents(1), &
 
87
                                           meridionalZonalVerticalVectorComponents(2), &
 
88
                                           meridionalZonalVerticalVectorComponents(3), &
 
89
                                           LLH(1), &
 
90
                                           LLH(2), &
 
91
                                           LLH(3), &
 
92
                                           cartesianVectorComponents(1), &
 
93
                                           cartesianVectorComponents(2), &
 
94
                                           cartesianVectorComponents(3), &
 
95
                                           XYZ(1), &
 
96
                                           XYZ(2), &
 
97
                                           XYZ(3))
 
98
    call set(radialVectorDifference, node, cartesianVectorComponents)
 
99
  enddo
 
100
  call addto(radialVectorDifference, UnitRadialVector_inCartesian, -1.0)
 
101
  fail = any(radialVectorDifference%val > 1e-12)
 
102
  call report_test( &
 
103
   "[Vector basis change: Zonal-Meridional-Vertical to Cartesian of unit-radial vector.]", &
 
104
   fail, .false., "Radial unit vector components not transformed correctly.")
 
105
 
 
106
  !Convert unit-polar vector components into Cartesian basis. Then compare
 
107
  ! with vector already in Cartesian basis, obtained from vtu.
 
108
  do node=1,node_count(LonLatHeightCoordinate)
 
109
    LLH = node_val(LonLatHeightCoordinate, node)
 
110
    meridionalZonalVerticalVectorComponents = &
 
111
                   node_val(UnitPolarVector_inZonalMeridionalRadial, node)
 
112
    call vector_lon_lat_height_2_cartesian(meridionalZonalVerticalVectorComponents(1), &
 
113
                                           meridionalZonalVerticalVectorComponents(2), &
 
114
                                           meridionalZonalVerticalVectorComponents(3), &
 
115
                                           LLH(1), &
 
116
                                           LLH(2), &
 
117
                                           LLH(3), &
 
118
                                           cartesianVectorComponents(1), &
 
119
                                           cartesianVectorComponents(2), &
 
120
                                           cartesianVectorComponents(3), &
 
121
                                           XYZ(1), &
 
122
                                           XYZ(2), &
 
123
                                           XYZ(3))
 
124
    call set(polarVectorDifference, node, cartesianVectorComponents)
 
125
  enddo
 
126
  call addto(polarVectorDifference, UnitPolarVector_inCartesian, -1.0)
 
127
  fail = any(polarVectorDifference%val > 1e-12)
 
128
  call report_test( &
 
129
   "[Vector basis change: Zonal-Meridional-Vertical to Cartesian of unit-polar vector.]", &
 
130
   fail, .false., "Polar unit vector components not transformed correctly.")
 
131
 
 
132
  !Convert unit-azimuthal vector components into Cartesian basis. Then compare
 
133
  ! with vector already in Cartesian basis, obtained from vtu.
 
134
  do node=1,node_count(LonLatHeightCoordinate)
 
135
    LLH = node_val(LonLatHeightCoordinate, node)
 
136
    meridionalZonalVerticalVectorComponents = &
 
137
                   node_val(UnitAzimuthalVector_inZonalMeridionalRadial, node)
 
138
    call vector_lon_lat_height_2_cartesian(meridionalZonalVerticalVectorComponents(1), &
 
139
                                           meridionalZonalVerticalVectorComponents(2), &
 
140
                                           meridionalZonalVerticalVectorComponents(3), &
 
141
                                           LLH(1), &
 
142
                                           LLH(2), &
 
143
                                           LLH(3), &
 
144
                                           cartesianVectorComponents(1), &
 
145
                                           cartesianVectorComponents(2), &
 
146
                                           cartesianVectorComponents(3), &
 
147
                                           XYZ(1), &
 
148
                                           XYZ(2), &
 
149
                                           XYZ(3))
 
150
    call set(azimuthalVectorDifference, node, cartesianVectorComponents)
 
151
  enddo
 
152
  call addto(azimuthalVectorDifference, UnitAzimuthalVector_inCartesian, -1.0)
 
153
  fail = any(azimuthalVectorDifference%val > 1e-12)
 
154
  call report_test( &
 
155
   "[Vector basis change: Zonal-Meridional-Vertical to Cartesian of unit-azimuthal vector.]", &
 
156
   fail, .false., "Azimuthal unit vector components not transformed correctly.")
 
157
 
 
158
  call deallocate(radialVectorDifference)
 
159
  call deallocate(polarVectorDifference)
 
160
  call deallocate(azimuthalVectorDifference)
 
161
 
 
162
end subroutine