Octopus
curv_briggs.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
25
28 use debug_oct_m
29 use global_oct_m
32 use parser_oct_m
34
35 implicit none
36
37 private
38 public :: &
41
42 type, extends(coordinate_system_t) :: curv_briggs_t
43 private
44 real(real64), allocatable :: lsize(:)
45 real(real64) :: beta
46 contains
47 procedure :: to_cartesian => curv_briggs_to_cartesian
48 procedure :: from_cartesian => curv_briggs_from_cartesian
49 procedure :: write_info => curv_briggs_write_info
50 procedure :: surface_element => curv_briggs_surface_element
51 procedure :: jacobian => curv_briggs_jacobian
52 procedure :: jacobian_inverse => curv_briggs_jacobian_inverse
54 end type curv_briggs_t
55
56 interface curv_briggs_t
57 procedure curv_briggs_constructor
58 end interface curv_briggs_t
59
60contains
61
62 ! ---------------------------------------------------------
63 function curv_briggs_constructor(namespace, dim, lsize, spacing) result(briggs)
64 type(namespace_t), intent(in) :: namespace
65 integer, intent(in) :: dim
66 real(real64), intent(in) :: lsize(1:dim)
67 real(real64), intent(in) :: spacing(1:dim)
68 class(curv_briggs_t), pointer :: briggs
69
70 integer :: idim
71
73
74 safe_allocate(briggs)
75
76 briggs%local_basis = .true.
77 briggs%orthogonal = .true. ! This needs to be checked
78 briggs%dim = dim
79 safe_allocate(briggs%lsize(1:dim))
80 briggs%lsize(1:dim) = lsize(1:dim)
81
82 !%Variable CurvBriggsBeta
83 !%Type float
84 !%Default 0.5
85 !%Section Mesh::Curvilinear::Briggs
86 !%Description
87 !% Controls the Briggs curvilinear coordinates
88 !%End
89 call parse_variable(namespace, 'CurvBriggsBeta', m_half, briggs%beta)
90
91 if (briggs%beta < m_zero .or. briggs%beta > m_one) then
92 message(1) = 'The parameter "CurvBriggsBeta" must lie between 0 and 1.'
93 call messages_fatal(1, namespace=namespace)
94 end if
95
96 briggs%min_mesh_scaling_product = m_one
97 do idim = 1, briggs%dim
98 ! corresponds to the distance of grid points at [+spacing/2,-spacing/2]
99 briggs%min_mesh_scaling_product = briggs%min_mesh_scaling_product * (m_one / &
100 (m_one - briggs%lsize(idim) * briggs%beta / (m_pi * spacing(idim)) * sin(m_pi * spacing(idim) / briggs%lsize(idim))))
101 end do
102
104 end function curv_briggs_constructor
105
106 ! ---------------------------------------------------------
107 subroutine curv_briggs_finalize(this)
108 type(curv_briggs_t), intent(inout) :: this
109
110 push_sub(curv_briggs_finalize)
111
112 safe_deallocate_a(this%lsize)
113
114 pop_sub(curv_briggs_finalize)
115 end subroutine curv_briggs_finalize
116
117 ! ---------------------------------------------------------
118 subroutine curv_briggs_copy(this_out, this_in)
119 type(curv_briggs_t), intent(inout) :: this_out
120 type(curv_briggs_t), intent(in) :: this_in
121
122 push_sub(curv_briggs_copy)
123
124 safe_allocate_source_a(this_out%lsize, this_in%lsize)
125 this_out%beta = this_in%beta
126
127 pop_sub(curv_briggs_copy)
128 end subroutine curv_briggs_copy
129
130 ! ---------------------------------------------------------
131 pure function curv_briggs_to_cartesian(this, chi) result(xx)
132 class(curv_briggs_t), target, intent(in) :: this
133 real(real64), intent(in) :: chi(:)
134 real(real64) :: xx(1:this%dim)
136 ! no PUSH_SUB, called too often
138 xx = chi - this%lsize*this%beta/(m_two*m_pi)*sin(m_two*m_pi*chi/this%lsize)
139
142 ! ---------------------------------------------------------
143 function curv_briggs_from_cartesian(this, xx) result(chi)
144 class(curv_briggs_t), target, intent(in) :: this
145 real(real64), intent(in) :: xx(:)
146 real(real64) :: chi(1:this%dim)
147
148 ! no PUSH_SUB, called too often
149
150 chi = m_zero
151 message(1) = "Internal error in curv_briggs_from_cartesian"
152 call messages_fatal(1)
153
154 end function curv_briggs_from_cartesian
155
156 ! ---------------------------------------------------------
157 subroutine curv_briggs_write_info(this, iunit, namespace)
158 class(curv_briggs_t), intent(in) :: this
159 integer, optional, intent(in) :: iunit
160 type(namespace_t), optional, intent(in) :: namespace
161
162 push_sub(curv_briggs_write_info)
163
164 write(message(1), '(a)') ' Curvilinear Method = briggs'
165 call messages_info(1, iunit=iunit, namespace=namespace)
166
168 end subroutine curv_briggs_write_info
169
170 ! ---------------------------------------------------------
171 real(real64) function curv_briggs_surface_element(this, idir) result(ds)
172 class(curv_briggs_t), intent(in) :: this
173 integer, intent(in) :: idir
174
175 ds = m_zero
176 message(1) = 'Surface element with briggs curvilinear coordinates not implemented'
177 call messages_fatal(1)
178
179 end function curv_briggs_surface_element
180
181 ! ---------------------------------------------------------
182 function curv_briggs_jacobian(this, chi) result(jacobian)
183 class(curv_briggs_t), intent(in) :: this
184 real(real64), intent(in) :: chi(:)
185 real(real64) :: jacobian(1:this%dim, 1:this%dim)
186
187 integer :: idim
188
189 ! Jacobian is diagonal in this method
190 jacobian(:, :) = m_zero
191 do idim = 1, this%dim
192 jacobian(idim, idim) = m_one - this%beta*cos(m_two*m_pi*chi(idim)/this%lsize(idim))
193 end do
194
195 end function curv_briggs_jacobian
196
197 ! ---------------------------------------------------------
198 function curv_briggs_jacobian_inverse(this, chi) result(jacobian_inverse)
199 class(curv_briggs_t), intent(in) :: this
200 real(real64), intent(in) :: chi(:)
201 real(real64) :: jacobian_inverse(1:this%dim, 1:this%dim)
202
203 integer :: idim
204
205 jacobian_inverse = this%jacobian(chi)
206 do idim = 1, this%dim
207 jacobian_inverse(idim, idim) = m_one/jacobian_inverse(idim, idim)
208 end do
209
212end module curv_briggs_oct_m
213
214!! Local Variables:
215!! mode: f90
216!! coding: utf-8
217!! End:
double sin(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
This module implements the curvilinear coordinates given in E.L. Briggs, D.J. Sullivan,...
real(real64) function, dimension(1:this%dim, 1:this%dim) curv_briggs_jacobian(this, chi)
class(curv_briggs_t) function, pointer curv_briggs_constructor(namespace, dim, lsize, spacing)
subroutine, public curv_briggs_copy(this_out, this_in)
real(real64) function curv_briggs_surface_element(this, idir)
real(real64) function, dimension(1:this%dim, 1:this%dim) curv_briggs_jacobian_inverse(this, chi)
subroutine curv_briggs_finalize(this)
pure real(real64) function, dimension(1:this%dim) curv_briggs_to_cartesian(this, chi)
subroutine curv_briggs_write_info(this, iunit, namespace)
real(real64) function, dimension(1:this%dim) curv_briggs_from_cartesian(this, xx)
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
abstract class to describe coordinate systems
int true(void)