[cig-commits] [commit] master: Fix a treatment of non-solenoidal vector averages over a sphere (9b73e1f)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Apr 8 03:04:36 PDT 2014
Repository : ssh://geoshell/calypso
On branch : master
Link : https://github.com/geodynamics/calypso/compare/bf5dcaf71a4089a4c2f22940f4edea71b9abedd1...9730b061d69d156b271dfe841a55ae371b4e1c03
>---------------------------------------------------------------
commit 9b73e1f4a80f9c21bec13ad5c3cf68d9f110c58c
Author: Hiroaki Matsui <h_kemono at mac.com>
Date: Thu Mar 20 13:27:16 2014 -0700
Fix a treatment of non-solenoidal vector averages over a sphere
>---------------------------------------------------------------
9b73e1f4a80f9c21bec13ad5c3cf68d9f110c58c
src/Fortran_libraries/MHD_src/Makefile | 2 +-
.../MHD_src/sph_MHD/cal_div_buoyancies_sph_MHD.f90 | 16 ++++++++++++
.../MHD_src/sph_MHD/const_sph_radial_grad.f90 | 2 ++
.../MHD_src/sph_MHD/set_reference_sph_mhd.f90 | 4 +--
.../SPH_SHELL_src/cal_rms_by_sph_spectr.f90 | 23 ++++++++++++++---
.../SPH_SHELL_src/cal_sph_exp_1st_diff.f90 | 30 +++++++++++++---------
.../SPH_SHELL_src/pickup_sph_coefs.f90 | 11 +++++---
.../SPH_SHELL_src/pickup_sph_spectr.f90 | 26 ++++++++++++++++++-
.../SERIAL_src/IO/m_pickup_sph_spectr_data.f90 | 6 ++++-
.../spherical_harmonics/spherical_harmonics.f90 | 4 ++-
10 files changed, 99 insertions(+), 25 deletions(-)
diff --git a/src/Fortran_libraries/MHD_src/Makefile b/src/Fortran_libraries/MHD_src/Makefile
index de4a2f2..fc821fa 100644
--- a/src/Fortran_libraries/MHD_src/Makefile
+++ b/src/Fortran_libraries/MHD_src/Makefile
@@ -55,7 +55,7 @@ mod_list:
@echo MOD_MHD= \\ >> $(MAKENAME)
@echo '$$(MOD_MHD_IO)' \\ >> $(MAKENAME)
@echo '$$(MOD_SPH_MHD)' \\ >> $(MAKENAME)
- @echo '$$(MOD_SPH_SNAPSHOT)' \\ >> $(MAKENAME)
+ @echo '$$(MOD_SPH_SNAPSHOT)' >> $(MAKENAME)
@echo >> $(MAKENAME)
@for dir in $(SUBDIRS); do \
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_div_buoyancies_sph_MHD.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_div_buoyancies_sph_MHD.f90
index 3aca9cf..3e0272d 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_div_buoyancies_sph_MHD.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_div_buoyancies_sph_MHD.f90
@@ -128,6 +128,14 @@
end do
!$omp end parallel do
!
+ if (idx_rj_degree_zero .eq. 0) return
+!$omp parallel do private (inod,k)
+ do k = 1, nidx_rj(1)
+ inod = idx_rj_degree_zero + (k-1)*nidx_rj(2)
+ d_rj(inod,is_div) = half * d_rj(inod,is_div)
+ end do
+!$omp end parallel do
+!
end subroutine cal_div_double_buoyancy_sph_MHD
!
!-----------------------------------------------------------------------
@@ -153,6 +161,14 @@
end do
!$omp end parallel do
!
+ if (idx_rj_degree_zero .eq. 0) return
+!$omp parallel do private (inod,k)
+ do k = 1, nidx_rj(1)
+ inod = idx_rj_degree_zero + (k-1)*nidx_rj(2)
+ d_rj(inod,is_div) = half * d_rj(inod,is_div)
+ end do
+!$omp end parallel do
+!
end subroutine cal_div_buoyancy_sph_MHD
!
!-----------------------------------------------------------------------
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_radial_grad.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_radial_grad.f90
index 22410f5..9db3137 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_radial_grad.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/const_sph_radial_grad.f90
@@ -83,6 +83,7 @@
call cal_sph_nod_gradient_2(sph_bc%kr_in, sph_bc%kr_out, &
& d_rj(1,is_fld), d_rj(1,is_grad) )
call sel_bc_radial_grad_scalar(sph_bc, is_fld, is_grad)
+ call normalize_sph_average_grad(d_rj(1,is_grad))
!
end subroutine const_radial_grad_scalar
!
@@ -193,6 +194,7 @@
!
call cal_sph_nod_gradient_2(sph_bc_U%kr_in, sph_bc_U%kr_out, &
& d_rj(1,is_press), d_rj(1,is_grad) )
+ call normalize_sph_average_grad(d_rj(1,is_grad))
!
call delete_bc_rj_vector(nidx_rj(2), sph_bc_U%kr_in, is_grad)
call delete_bc_rj_vector(nidx_rj(2), sph_bc_U%kr_out, is_grad)
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_reference_sph_mhd.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_reference_sph_mhd.f90
index ffe3a86..0717af1 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_reference_sph_mhd.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_reference_sph_mhd.f90
@@ -160,7 +160,7 @@
d_rj(inod,ipol%i_temp) = d_rj(inod,ipol%i_temp) &
& - reftemp_rj(k,0)
d_rj(inod,ipol%i_grad_t) = d_rj(inod,ipol%i_grad_t) &
- & - reftemp_rj(k,1) * radius_1d_rj_r(k)**2
+ & - two*reftemp_rj(k,1) * radius_1d_rj_r(k)**2
d_rj(inod,idpdr%i_grad_t) = d_rj(inod,ipol%i_temp)
end do
end if
@@ -199,7 +199,7 @@
d_rj(inod,ipol%i_temp) = d_rj(inod,ipol%i_temp) &
& + reftemp_rj(k,0)
d_rj(inod,ipol%i_grad_t) = d_rj(inod,ipol%i_grad_part_t) &
- & + reftemp_rj(k,1) * radius_1d_rj_r(k)**2
+ & + two*reftemp_rj(k,1) * radius_1d_rj_r(k)**2
d_rj(inod,idpdr%i_grad_t) = d_rj(inod,ipol%i_temp)
end do
end if
diff --git a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/cal_rms_by_sph_spectr.f90 b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/cal_rms_by_sph_spectr.f90
index 326e768..c136f42 100644
--- a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/cal_rms_by_sph_spectr.f90
+++ b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/cal_rms_by_sph_spectr.f90
@@ -84,6 +84,23 @@
end do
!$omp end parallel do
!
+ j = idx_rj_degree_zero
+ if(idx_rj_degree_zero .eq. izero) then
+ do k = 1, nidx_rj(1)
+ rms_sph_dat(j,k,jcomp ) = zero
+ rms_sph_dat(j,k,jcomp+1) = zero
+ rms_sph_dat(j,k,jcomp+2) = zero
+ end do
+ else
+ do k = 1, nidx_rj(1)
+ idx = idx_rj_degree_zero + (k-1) * nidx_rj(2)
+ rms_sph_dat(j,k,jcomp ) = (half * d_rj(idx,icomp))**2 &
+ & * a_r_1d_rj_r(k)*a_r_1d_rj_r(k)
+ rms_sph_dat(j,k,jcomp+1) = zero
+ rms_sph_dat(j,k,jcomp+2) = rms_sph_dat(j,k,jcomp )
+ end do
+ end if
+!
end subroutine cal_rms_each_vector_sph_spec
!
! -----------------------------------------------------------------------
@@ -150,11 +167,9 @@
do k = 1, nidx_rj(1)
kg = idx_gl_1d_rj_r(k)
inod = idx_rj_degree_zero + (k-1) * nidx_rj(2)
- ave_sph_lc(kg,jcomp ) = d_rj(inod,icomp) &
- & * radius_1d_rj_r(kg) * radius_1d_rj_r(kg)
+ ave_sph_lc(kg,jcomp ) = half * d_rj(inod,icomp)
ave_sph_lc(kg,jcomp+1) = zero
- ave_sph_lc(kg,jcomp+2) = d_rj(inod,icomp) &
- & * radius_1d_rj_r(kg) * radius_1d_rj_r(kg)
+ ave_sph_lc(kg,jcomp+2) = ave_sph_lc(kg,jcomp )
end do
end if
!
diff --git a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/cal_sph_exp_1st_diff.f90 b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/cal_sph_exp_1st_diff.f90
index 2222c63..5da2e73 100644
--- a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/cal_sph_exp_1st_diff.f90
+++ b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/cal_sph_exp_1st_diff.f90
@@ -9,6 +9,7 @@
!!@verbatim
!! subroutine cal_sph_nod_gradient_2(kr_in, kr_out, &
!! & dnod_rj, dnod_dr)
+!! subroutine normalize_sph_average_grad(dnod_dr)
!! subroutine cal_sph_nod_vect_dr_2(kr_in, kr_out, dnod_rj, dnod_dr)
!!@endverbatim
!!
@@ -71,24 +72,29 @@
end do
!$omp end parallel do
!
- if(idx_rj_degree_zero .eq. 0) return
+ end subroutine cal_sph_nod_gradient_2
!
- j = idx_rj_degree_zero
-!$omp parallel do private(inod,i_p1,i_n1,k,d1sdr)
- do k = kr_in+1, kr_out-1
- inod = (k-1) * nidx_rj(2) + idx_rj_degree_zero
- i_p1 = inod + nidx_rj(2)
- i_n1 = inod - nidx_rj(2)
+! -----------------------------------------------------------------------
!
- d1sdr = d1nod_mat_fdm_2(k,-1) * dnod_rj(i_n1) &
- & + d1nod_mat_fdm_2(k, 0) * dnod_rj(inod) &
- & + d1nod_mat_fdm_2(k, 1) * dnod_rj(i_p1)
+ subroutine normalize_sph_average_grad(dnod_dr)
!
- dnod_dr(inod,1) = d1sdr * radius_1d_rj_r(k)**2
+ real(kind = kreal), intent(inout) :: dnod_dr(nnod_rj,3)
+!
+ integer(kind = kint) :: inod, k
+!
+!
+ if(idx_rj_degree_zero .eq. 0) return
+!
+!$omp parallel do private(inod,k)
+ do k = 1, nidx_rj(1)
+ inod = (k-1) * nidx_rj(2) + idx_rj_degree_zero
+ dnod_dr(inod,1) = two * dnod_dr(inod,1)
+ dnod_dr(inod,2) = zero
+ dnod_dr(inod,3) = zero
end do
!$omp end parallel do
!
- end subroutine cal_sph_nod_gradient_2
+ end subroutine normalize_sph_average_grad
!
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
diff --git a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/pickup_sph_coefs.f90 b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/pickup_sph_coefs.f90
index ea51244..7a880fd 100644
--- a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/pickup_sph_coefs.f90
+++ b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/pickup_sph_coefs.f90
@@ -66,6 +66,8 @@
& idx_pick_sph_mode, idx_pick_sph_l, idx_pick_sph_m, &
& ntot_pick_sph_mode, num_pick_sph_mode, idx_pick_sph_gl, &
& idx_pick_sph_lc)
+ call set_scale_4_vect_l0 &
+ & (num_pick_sph_mode, idx_pick_sph_gl(1), scale_for_zelo(1))
call deallocate_iflag_pick_sph
call deallocate_pick_sph_mode
!
@@ -105,9 +107,12 @@
icou = istack_phys_comp_rj(i_fld-1)
jcou = istack_comp_pick_sph(j_fld-1)
if(num_phys_comp_rj(i_fld) .eq. 3) then
- d_rj_pick_sph_lc(jcou+1,ipick)= d_rj(inod,icou+1)
- d_rj_pick_sph_lc(jcou+2,ipick)= d_rj(inod,icou+3)
- d_rj_pick_sph_lc(jcou+3,ipick)= d_rj(inod,icou+2)
+ d_rj_pick_sph_lc(jcou+1,ipick) &
+ & = scale_for_zelo(inum) * d_rj(inod,icou+1)
+ d_rj_pick_sph_lc(jcou+2,ipick) &
+ & = scale_for_zelo(inum) * d_rj(inod,icou+3)
+ d_rj_pick_sph_lc(jcou+3,ipick) &
+ & = scale_for_zelo(inum) * d_rj(inod,icou+2)
else
do nd = 1, num_phys_comp_rj(i_fld)
d_rj_pick_sph_lc(jcou+nd,ipick) = d_rj(inod,icou+nd)
diff --git a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/pickup_sph_spectr.f90 b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/pickup_sph_spectr.f90
index 22a9ea3..571a5f2 100644
--- a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/pickup_sph_spectr.f90
+++ b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/pickup_sph_spectr.f90
@@ -12,7 +12,9 @@
!! subroutine set_picked_sph_address(l_truncation, jmax, idx_gl_j, &
!! & num_pick_sph, num_pick_sph_l, num_pick_sph_m, &
!! & idx_pick_sph, idx_pick_sph_l, idx_pick_sph_m, &
-!! & ntot_pickup, num_pickup, idx_pick_gl, idx_pick_lc)
+!! & ntot_pickup, num_pickup, idx_pick_gl, idx_pick_lc, &
+!! subroutine set_scale_4_vect_l0(num_pickup, &
+!! & idx_pick_gl, scale_for_zelo)
!
module pickup_sph_spectr
!
@@ -165,4 +167,26 @@
!
! -----------------------------------------------------------------------
!
+ subroutine set_scale_4_vect_l0(num_pickup, &
+ & idx_pick_gl, scale_for_zelo)
+!
+ integer(kind = kint), intent(in) :: num_pickup
+ integer(kind = kint), intent(in) :: idx_pick_gl(num_pickup)
+ real(kind = kreal), intent(inout) :: scale_for_zelo(num_pickup)
+!
+ integer(kind = kint) :: inum
+!
+!
+ do inum = 1, num_pickup
+ if(idx_pick_gl(inum) .eq. 0) then
+ scale_for_zelo(inum) = half
+ else
+ scale_for_zelo(inum) = one
+ end if
+ end do
+!
+ end subroutine set_scale_4_vect_l0
+!
+! -----------------------------------------------------------------------
+!
end module pickup_sph_spectr
diff --git a/src/Fortran_libraries/SERIAL_src/IO/m_pickup_sph_spectr_data.f90 b/src/Fortran_libraries/SERIAL_src/IO/m_pickup_sph_spectr_data.f90
index 7e6c859..e20ed79 100644
--- a/src/Fortran_libraries/SERIAL_src/IO/m_pickup_sph_spectr_data.f90
+++ b/src/Fortran_libraries/SERIAL_src/IO/m_pickup_sph_spectr_data.f90
@@ -83,6 +83,8 @@
real(kind = kreal), allocatable :: d_rj_pick_sph_gl(:,:)
!> Localy evaluated monitoring spectrum
real(kind = kreal), allocatable :: d_rj_pick_sph_lc(:,:)
+!> Scale factor for vector at l=m=0
+ real(kind = kreal), allocatable :: scale_for_zelo(:)
!> Name of monitoring spectrum
character(len=kchara), allocatable :: pick_sph_spec_name(:)
!
@@ -148,6 +150,7 @@
allocate( idx_pick_sph_lc(ntot_pick_sph_mode) )
allocate( d_rj_pick_sph_lc(ntot_comp_pick_sph,num) )
allocate( d_rj_pick_sph_gl(ntot_comp_pick_sph,num) )
+ allocate( scale_for_zelo(ntot_comp_pick_sph) )
allocate( pick_sph_spec_name(ntot_comp_pick_sph) )
allocate( istack_comp_pick_sph(0:num_fld_pick_sph) )
allocate( ifield_monitor_rj(num_fld_pick_sph) )
@@ -161,6 +164,7 @@
idx_pick_sph_lc = 0
d_rj_pick_sph_lc = 0.0d0
d_rj_pick_sph_gl = 0.0d0
+ scale_for_zelo = 1.0d0
end if
!
end subroutine allocate_pick_sph_monitor
@@ -190,7 +194,7 @@
!
!
deallocate(idx_pick_sph_gl, d_rj_pick_sph_gl)
- deallocate(idx_pick_sph_lc, d_rj_pick_sph_lc)
+ deallocate(idx_pick_sph_lc, d_rj_pick_sph_lc, scale_for_zelo)
deallocate(pick_sph_spec_name)
deallocate(istack_comp_pick_sph, ifield_monitor_rj)
!
diff --git a/src/Fortran_libraries/SERIAL_src/spherical_harmonics/spherical_harmonics.f90 b/src/Fortran_libraries/SERIAL_src/spherical_harmonics/spherical_harmonics.f90
index 3b4eb91..bad188f 100644
--- a/src/Fortran_libraries/SERIAL_src/spherical_harmonics/spherical_harmonics.f90
+++ b/src/Fortran_libraries/SERIAL_src/spherical_harmonics/spherical_harmonics.f90
@@ -30,6 +30,7 @@
!* g(j,11) : 1 / (2*l+1)
!* g(j,12) : l*(l+1) / (2*l+1)
!* g(j,13) : 1 / (l*(l+1))
+!* g(j,13) : 1 (l=m=0)
!
!* g(j,16) : (2*l+1) / 4
!* g(j,17) : (2*l+1) / ( 4*l*(l+1) )
@@ -37,6 +38,7 @@
!
!
! Note: g(0,3) = 1/2 for spherical harmonics transform
+! Note: g(0,13) = 2 for gradient of scalar
! (See g_sph_rlm in schmidt_poly_on_rtm_grid.f90)
!
!*
@@ -120,7 +122,7 @@
!
if ( j .eq. 0 ) then
g(j,9) = one / (four*pi)
- g(j,13) = zero
+ g(j,13) = one
else
g(j,9) = dble(2*l+1) / (four*pi*dble( l*(l+1) ))
g(j,13) = one / dble( l*(l+1) )
More information about the CIG-COMMITS
mailing list