33 use,
intrinsic :: iso_fortran_env
52 real(real64),
allocatable :: lsize(:)
56 real(real64),
allocatable :: Jlocal(:)
57 real(real64),
allocatable :: Jrange(:)
59 real(real64),
allocatable :: xi(:,:)
60 real(real64),
allocatable :: Q(:,:,:)
63 type(root_solver_t) :: rs, rs_init
75 procedure curv_modine_constructor
79 class(curv_modine_t),
pointer :: modine_p
80 real(real64),
allocatable :: x_p(:)
86 type(namespace_t),
intent(in) :: namespace
87 integer,
intent(in) :: dim
88 integer,
intent(in) :: npos
89 real(real64),
intent(in) :: pos(1:dim,1:npos)
90 real(real64),
intent(in) :: lsize(1:dim)
91 real(real64),
intent(in) :: spacing(1:dim)
92 class(curv_modine_t),
pointer :: modine
95 real(real64) :: kappa(1:npos)
100 safe_allocate(modine)
103 modine%local_basis = .
true.
104 modine%orthogonal = .
true.
108 safe_allocate(modine%lsize(1:dim))
131 if (modine%xbar <
m_zero .or. modine%xbar >
m_one)
then
132 message(1) =
'The parameter "CurvModineXBar" must lie between 0 and 1.'
135 if (modine%Jbar <= modine%xbar)
then
136 message(1) =
'The parameter "CurvModineJBar" must be larger than the parameter "CurvModineXBar".'
140 safe_allocate(modine%Jlocal(1:npos))
141 safe_allocate(modine%Jrange(1:npos))
156 modine%Jlocal(:) = modine%Jlocal(1)
157 else if (n == npos)
then
162 write(
message(1),
'(3a,i2)')
'Error in block CurvModineJlocal: number of lines must be 1 or ', npos
167 modine%Jlocal(:) = 0.125_real64
171 if (modine%Jlocal(iatom)<
m_zero.or.modine%Jlocal(iatom)>
m_one)
then
172 message(1) =
'The parameter "CurvModineJlocal" must lie between 0 and 1.'
186 if (
parse_block(namespace,
"CurvModineJrange", blk) == 0)
then
191 else if (n == npos)
then
196 write(
message(1),
'(3a,i2)')
'Error in block CurvModineJrange: number of lines must be 1 or ', npos
201 kappa(:) = 1.2_real64
206 if (modine%Jlocal(iatom) >
m_one-1e-5_real64)
then
207 modine%Jrange(iatom) = kappa(iatom)
209 modine%Jrange(iatom) = kappa(iatom) / &
217 solver_type =
root_newton, maxiter = 500, abs_tolerance = 1.0e-13_real64)
219 safe_allocate(modine%xi(1:modine%dim, 1:modine%npos))
220 safe_allocate(modine%Q(1:modine%dim, 1:modine%dim, 1:modine%npos))
223 modine%min_mesh_scaling_product = modine%Jbar**modine%dim
230 real(real64),
intent(in) :: x
239 integer :: iatom, idim
240 real(real64),
allocatable :: xi_tmp(:, :), q_tmp(:, :, :)
241 real(real64) :: residue, rr2, dd
242 integer :: iter, nu, mu, idim1, idim2
243 integer,
parameter :: maxiter = 1000
244 real(real64),
parameter :: tolerance = 1.0e-10_real64
248 safe_allocate(xi_tmp(1:modine%dim, 1:modine%npos))
249 safe_allocate(q_tmp(1:modine%dim, 1:modine%dim, 1:modine%npos))
252 modine%Q(:, :, :) =
m_zero
253 do iatom = 1, modine%npos
254 modine%xi(:, iatom) = pos(:, iatom)
255 do idim = 1, modine%dim
256 modine%Q(idim, idim, iatom) =
m_one-modine%Jlocal(iatom)**(
m_one/real(modine%dim, real64))
263 xi_tmp(:, :) = modine%xi(:, :)
264 q_tmp(:, :, :) = modine%Q(:, :, :)
266 modine%Q(:, :, :) =
m_zero
267 do mu = 1, modine%npos
269 modine%xi(:, mu) = pos(:, mu)
270 do idim = 1, modine%dim
272 modine%Q(idim, idim, mu) =
m_one-modine%Jlocal(mu)**(
m_one/real(modine%dim, real64))
274 do nu = 1, modine%npos
276 rr2 = sum((xi_tmp(:, mu) - xi_tmp(:, nu))**2)
277 dd =
exp(-rr2/(
m_two*modine%Jrange(nu)**2))
279 modine%xi(:, mu) = modine%xi(:, mu) + &
280 matmul(modine%Q(:, :, mu), xi_tmp(:, mu) - xi_tmp(:, nu)) * dd
282 do idim1 = 1, modine%dim
283 do idim2 = 1, modine%dim
284 modine%Q(idim1, idim2, mu) = modine%Q(idim1, idim2, mu) - q_tmp(idim1, idim2, nu) * dd
285 modine%Q(idim1, idim2, mu) = modine%Q(idim1, idim2, mu) + &
286 sum(q_tmp(idim1, :, nu)*(xi_tmp(:, mu) - xi_tmp(:, nu)))*(xi_tmp(idim2, mu) - xi_tmp(idim2, nu)) * &
287 m_one/modine%Jrange(nu) * dd
294 residue = norm2(xi_tmp(:, :) - modine%xi(:, :)) / real(modine%dim*modine%npos, real64)
295 residue = residue + norm2(q_tmp(:, :, :) - modine%Q(:, :, :)) / real(modine%dim*modine%dim*modine%npos, real64)
296 if (residue <= tolerance)
exit
300 if (residue > tolerance)
then
301 message(1) =
"During the construction of the adaptive grid, the relaxation"
302 message(2) =
"method did not converge."
303 message(3) =
"This might be due to values of CurvModineJrange that are too large"
304 message(4) =
"or CurvModineJlocal that are too small."
309 safe_deallocate_a(xi_tmp)
310 safe_deallocate_a(q_tmp)
324 this_out%dim = this_in%dim
325 this_out%local_basis = this_in%local_basis
326 safe_allocate_source_a(this_out%lsize, this_in%lsize)
327 this_out%xbar = this_in%xbar
328 this_out%Jbar = this_in%Jbar
329 safe_allocate_source_a(this_out%Jlocal, this_in%Jlocal)
330 safe_allocate_source_a(this_out%Jrange, this_in%Jrange)
331 safe_allocate_source_a(this_out%xi, this_in%xi)
332 safe_allocate_source_a(this_out%Q, this_in%Q)
333 this_out%npos = this_in%npos
344 safe_deallocate_a(this%lsize)
345 safe_deallocate_a(this%Jlocal)
346 safe_deallocate_a(this%Jrange)
347 safe_deallocate_a(this%xi)
348 safe_deallocate_a(this%Q)
356 real(real64),
intent(in) :: chi(:)
357 real(real64) :: xx(1:this%dim)
359 real(real64) :: chi2(this%dim), rr2, dd
369 do iatom = 1, this%npos
370 rr2 = sum((chi2(:) - this%xi(:, iatom))**2)
371 dd =
exp(-rr2/(m_two*this%Jrange(iatom)**2))
373 xx(:) = xx(:) - matmul(this%Q(:, :, iatom), chi2(:) - this%xi(:, iatom)) * dd
381 real(real64),
intent(in) :: xx(:)
382 real(real64) :: chi(1:this%dim)
385 real(real64) :: start(1:this%dim)
391 safe_allocate(
x_p(1:this%dim))
396 if (abs(xx(i)) > this%lsize(i))
then
398 else if (abs(xx(i)) > this%xbar*this%lsize(i))
then
400 start(i) = sign(m_one, xx(i)) * (abs(xx(i))/this%Jbar + abs(xx(i)) *(m_one - m_one/this%Jbar) * &
401 (abs(xx(i)) - this%xbar*this%lsize(i)) / (this%lsize(i) - this%xbar*this%lsize(i)))
403 if (abs(start(i)) >= this%lsize(i))
then
404 start(i) = sign(this%lsize(i), xx(i))
407 start(i) = xx(i)/this%Jbar
411 call droot_solver_run(this%rs,
getf, chi, conv, startval=start)
417 message(1) =
"During the construction of the adaptive grid, the Newton-Raphson"
418 message(2) =
"method did not converge for point:"
419 write(message(3),
'(3f14.6)') xx(1:this%dim)
420 call messages_fatal(3)
428 integer,
optional,
intent(in) :: iunit
429 type(namespace_t),
optional,
intent(in) :: namespace
433 write(message(1),
'(a)')
' Curvilinear Method = modine'
434 call messages_info(1, iunit=iunit, namespace=namespace)
440 real(real64) function curv_modine_surface_element(this, idir) result(ds)
442 integer,
intent(in) :: idir
445 message(1) =
'Surface element with modine curvilinear coordinates not implemented'
446 call messages_fatal(1)
453 real(real64),
intent(in) :: chi_(:)
454 real(real64),
intent(out) :: chi2(:)
455 real(real64),
optional,
intent(out) :: jac(:)
457 integer,
parameter :: qq = 3
458 real(real64) :: chibar(this%dim), rr, chi
464 chibar = this%xbar*this%lsize/this%Jbar
470 chi2(i) = this%Jbar * chi
471 if (
present(jac)) jac(i) = this%Jbar
473 if (chi > this%lsize(i))
then
476 if (
present(jac)) jac(i) = m_one
477 else if (chi > chibar(i))
then
478 rr = (chi - chibar(i))/(this%lsize(i) - chibar(i))
480 chi2(i) = chi2(i) + this%lsize(i)/m_two*(1 - this%Jbar) * rr**qq * &
481 (qq + m_one - (qq - m_one)*rr)
483 if (
present(jac))
then
484 jac(i) = jac(i) + this%lsize(i)/m_two*(m_one - this%Jbar) * rr**(qq - 1)/(this%lsize(i) - chibar(i)) * &
485 (qq*(qq + m_one) - (qq**2 - m_one)*rr)
489 if (neg) chi2(i) = -chi2(i)
499 real(real64),
intent(in) :: chi(:)
500 real(real64),
intent(out) :: xx(:)
501 real(real64),
intent(out) :: jac(:, :)
503 real(real64) :: chi2(this%dim), rr2, dd, j2(this%dim)
504 integer :: iatom, idim, idim2
513 do idim = 1, this%dim
514 jac(idim, idim) = m_one
517 do iatom = 1, this%npos
518 rr2 = sum((chi2(:) - this%xi(:, iatom))**2)
519 dd =
exp(-rr2/(m_two*this%Jrange(iatom)**2))
521 xx(:) = xx(:) - matmul(this%Q(:, :, iatom), chi2(:) - this%xi(:, iatom)) * dd
523 do idim = 1, this%dim
524 do idim2 = 1, this%dim
525 jac(idim, idim2) = jac(idim, idim2) - this%Q(idim, idim2, iatom) * dd
526 jac(idim, idim2) = jac(idim, idim2) + &
527 sum(this%Q(idim, :, iatom)*(chi2(:) - this%xi(:, iatom)))*(chi2(idim2) - this%xi(idim2, iatom)) * &
528 m_one/this%Jrange(iatom) * dd
533 do idim = 1, this%dim
534 jac(idim, :) = jac(idim, :) * j2(:)
540 pure subroutine getf(yy, ff, jf)
541 real(real64),
intent(in) :: yy(:)
542 real(real64),
intent(out) :: ff(:), jf(:, :)
547 ff(:) = ff(:) -
x_p(:)
554 real(real64),
intent(in) :: chi(:)
555 real(real64) :: jacobian(1:this%dim, 1:this%dim)
557 real(real64) :: dummy(this%dim)
566 real(real64),
intent(in) :: chi(:)
567 real(real64) :: jacobian_inverse(1:this%dim, 1:this%dim)
569 real(real64) :: dummy(this%dim)
572 call lalg_inverse(this%dim, jacobian_inverse,
"dir")
real(real64) pure function product_log(x)
double exp(double __x) __attribute__((__nothrow__
This module implements the curvilinear coordinates given in N. A. Modine, G. Zumbach,...
class(curv_modine_t), pointer modine_p
real(real64) function, dimension(1:this%dim) curv_modine_from_cartesian(this, xx)
real(real64) function curv_modine_surface_element(this, idir)
real(real64) function, dimension(1:this%dim, 1:this%dim) curv_modine_jacobian_inverse(this, chi)
pure subroutine curv_modine_jacobian_inv(this, chi, xx, Jac)
real(real64) function, dimension(1:this%dim, 1:this%dim) curv_modine_jacobian(this, chi)
class(curv_modine_t) function, pointer curv_modine_constructor(namespace, dim, npos, pos, lsize, spacing)
pure subroutine getf(yy, ff, jf)
subroutine curv_modine_finalize(this)
real(real64), dimension(:), allocatable x_p
pure subroutine curv_modine_chi2chi2(this, chi_, chi2, Jac)
subroutine, public curv_modine_copy(this_out, this_in)
subroutine curv_modine_write_info(this, iunit, namespace)
real(real64) function, dimension(1:this%dim) curv_modine_to_cartesian(this, chi)
real(real64), parameter, public m_two
real(real64), parameter, public m_zero
real(real64), parameter, public m_four
real(real64), parameter, public m_third
real(real64), parameter, public m_half
real(real64), parameter, public m_one
real(real64), parameter, public m_three
real(real64), parameter, public m_e
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
integer function, public parse_block(namespace, name, blk, check_varinfo_)
integer, parameter, public root_newton
subroutine, public root_solver_init(rs, namespace, dimensionality, solver_type, maxiter, rel_tolerance, abs_tolerance)
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
This module defines the unit system, used for input and output.
type(unit_system_t), public units_inp
the units systems for reading and writing
abstract class to describe coordinate systems