[cig-commits] [commit] master: Fixing a bug to find radial interpolation points in assemble_sph (3b19a61)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Apr 8 03:04:23 PDT 2014
Repository : ssh://geoshell/calypso
On branch : master
Link : https://github.com/geodynamics/calypso/compare/bf5dcaf71a4089a4c2f22940f4edea71b9abedd1...9730b061d69d156b271dfe841a55ae371b4e1c03
>---------------------------------------------------------------
commit 3b19a61cccc4ffdaf014edf2e3b10bfd9c4dc1a4
Author: Hiroaki Matsui <h_kemono at mac.com>
Date: Mon Mar 17 12:17:52 2014 -0700
Fixing a bug to find radial interpolation points in assemble_sph
>---------------------------------------------------------------
3b19a61cccc4ffdaf014edf2e3b10bfd9c4dc1a4
.../dynamobench_case_1/control_MHD | 1 -
.../COMM_src/parallel_ucd_IO_select.F90 | 13 ++--
.../SPH_SPECTR_src/extend_potential_field.f90 | 58 ++++++++++++++
.../UTILS_src/MERGE/merge_sph_step_spectr.f90 | 9 +++
.../UTILS_src/MERGE/r_interpolate_marged_sph.f90 | 88 +++++++++++++++++-----
.../INITIAL_FIELD/const_sph_initial_spectr.f90 | 12 +--
6 files changed, 149 insertions(+), 32 deletions(-)
diff --git a/examples/dynamo_benchmark/dynamobench_case_1/control_MHD b/examples/dynamo_benchmark/dynamobench_case_1/control_MHD
index dcb4184..2de0c14 100644
--- a/examples/dynamo_benchmark/dynamobench_case_1/control_MHD
+++ b/examples/dynamo_benchmark/dynamobench_case_1/control_MHD
@@ -28,7 +28,6 @@ begin MHD_control
field_file_prefix 'field/out'
!
field_file_fmt_ctl 'merged_VTK'
-! field_file_fmt_ctl 'VTK'
end data_files_def
!
begin model
diff --git a/src/Fortran_libraries/PARALLEL_src/COMM_src/parallel_ucd_IO_select.F90 b/src/Fortran_libraries/PARALLEL_src/COMM_src/parallel_ucd_IO_select.F90
index 4667784..efb847a 100644
--- a/src/Fortran_libraries/PARALLEL_src/COMM_src/parallel_ucd_IO_select.F90
+++ b/src/Fortran_libraries/PARALLEL_src/COMM_src/parallel_ucd_IO_select.F90
@@ -106,6 +106,8 @@
subroutine choose_para_fld_file_format(file_fmt_ctl, i_file_fmt, &
& id_field_file_format)
!
+ use skip_comment_f
+!
integer(kind= kint), intent(in) :: i_file_fmt
character(len=kchara), intent(in) :: file_fmt_ctl
integer(kind= kint), intent(inout) :: id_field_file_format
@@ -116,15 +118,10 @@
return
end if
!
- if( file_fmt_ctl.eq.'merged_vtk' &
- & .or. file_fmt_ctl.eq.'MERGED_VTK' &
- & .or. file_fmt_ctl.eq.'merged_vtk_ascii' &
- & .or. file_fmt_ctl.eq.'MERGED_VTK_ASCII') then
+ if( cmp_no_case(file_fmt_ctl, 'merged_VTK').gt.0 &
+ & .or. cmp_no_case(file_fmt_ctl, 'merged_vtk_ascii').gt.0) then
id_field_file_format = iflag_sgl_vtk
- else if(file_fmt_ctl.eq.'merged_hdf5' &
- & .or. file_fmt_ctl.eq.'merged_HDF5' &
- & .or. file_fmt_ctl.eq.'Merged_HDF5' &
- & .or. file_fmt_ctl.eq.'MERGED_HDF5') then
+ else if(cmp_no_case(file_fmt_ctl, 'merged_HDF5').gt. 0) then
id_field_file_format = iflag_sgl_hdf5
else
call choose_ucd_file_format(file_fmt_ctl, i_file_fmt, &
diff --git a/src/Fortran_libraries/SERIAL_src/SPH_SPECTR_src/extend_potential_field.f90 b/src/Fortran_libraries/SERIAL_src/SPH_SPECTR_src/extend_potential_field.f90
index 225e93d..740549d 100644
--- a/src/Fortran_libraries/SERIAL_src/SPH_SPECTR_src/extend_potential_field.f90
+++ b/src/Fortran_libraries/SERIAL_src/SPH_SPECTR_src/extend_potential_field.f90
@@ -3,6 +3,9 @@
!
! modified by H. Matsuiui on Apr., 2009
!
+! subroutine ext_outside_scalar(is_fld, kr_out)
+! subroutine ext_inside_scalar(is_fld, kr_in)
+!
! subroutine ext_outside_potential(is_fld, kr_out)
! Output: d_rj(kr_out+1:,is_fld:is_fld+2)
! subroutine ext_inside_potential(is_fld, kr_in)
@@ -36,6 +39,61 @@
!
! -------------------------------------------------------------------
!
+ subroutine ext_outside_scalar(is_fld, kr_out)
+!
+ integer(kind = kint), intent(in) :: is_fld, kr_out
+ real(kind = kreal) :: ratio
+ integer(kind = kint) :: inod, inod_cmb
+ integer(kind = kint) :: j, l_gl, k
+!
+!
+!$omp parallel private(k,ratio)
+ do k = kr_out+1, nidx_rj(1)
+ ratio = radius_1d_rj_r(kr_out) * a_r_1d_rj_r(k)
+!$omp do private(j,l_gl,inod,inod_cmb)
+ do j = 1, nidx_rj(2)
+ inod = j + (k-1) * nidx_rj(2)
+ inod_cmb = j + (kr_out-1) * nidx_rj(2)
+ l_gl = idx_gl_1d_rj_j(j,2)
+!
+ d_rj(inod,is_fld ) = d_rj(inod_cmb,is_fld) * ratio**l_gl
+ end do
+!$omp end do nowait
+ end do
+!$omp end parallel
+!
+ end subroutine ext_outside_scalar
+!
+! -------------------------------------------------------------------
+!
+ subroutine ext_inside_scalar(is_fld, kr_in)
+!
+ integer(kind = kint), intent(in) :: is_fld, kr_in
+!
+ real(kind = kreal) :: ratio
+ integer(kind = kint) :: inod, inod_icb
+ integer(kind = kint) :: j, l_gl, k
+!
+!
+!$omp parallel private(k,ratio)
+ do k = 1, kr_in-1
+ ratio = radius_1d_rj_r(k) * a_r_1d_rj_r(kr_in)
+!$omp do private(j,l_gl,inod,inod_icb)
+ do j = 1, nidx_rj(2)
+ inod = j + (k-1) * nidx_rj(2)
+ inod_icb = j + (kr_in-1) * nidx_rj(2)
+ l_gl = idx_gl_1d_rj_j(j,2)
+!
+ d_rj(inod,is_fld ) = d_rj(inod_icb,is_fld) * ratio**l_gl
+ end do
+!$omp end do nowait
+ end do
+!$omp end parallel
+!
+ end subroutine ext_inside_scalar
+!
+! -------------------------------------------------------------------
+!
subroutine ext_outside_potential(is_fld, kr_out)
!
integer(kind = kint), intent(in) :: is_fld, kr_out
diff --git a/src/Fortran_libraries/UTILS_src/MERGE/merge_sph_step_spectr.f90 b/src/Fortran_libraries/UTILS_src/MERGE/merge_sph_step_spectr.f90
index f5b4817..78aaf9f 100644
--- a/src/Fortran_libraries/UTILS_src/MERGE/merge_sph_step_spectr.f90
+++ b/src/Fortran_libraries/UTILS_src/MERGE/merge_sph_step_spectr.f90
@@ -16,6 +16,7 @@
module merge_sph_step_spectr
!
use m_precision
+ use m_machine_parameter
use m_constants
use m_field_data_IO
!
@@ -104,18 +105,26 @@
!
phys_file_head = org_sph_fst_head
!
+
do ip = 1, np_sph_org
my_rank = ip - 1
call sel_read_alloc_step_SPH_file(my_rank, istep)
call deallocate_phys_data_name_IO
!
if(iflag_same_rgrid .eq. 0) then
+! write(*,*) 'itp_rj_merged_phys_from_IO'
call itp_rj_merged_phys_from_IO &
& (org_sph_mesh(ip)%sph_mesh%sph_rj%nnod_rj, &
& org_sph_mesh(ip)%sph_mesh%sph_rj%nidx_rj(2), &
& org_sph_mesh(ip)%sph_mesh%sph_rj%idx_gl_1d_rj_j, &
& phys_data_IO)
+
+! write(*,*) 'extend_potential_magne'
call extend_potential_magne
+! write(*,*) 'extend_inner_core_temp'
+ call extend_inner_core_temp
+! write(*,*) 'extend_inner_core_light'
+ call extend_inner_core_light
else
call copy_rj_merged_phys_from_IO &
& (org_sph_mesh(ip)%sph_mesh%sph_rj%nidx_rj(2), &
diff --git a/src/Fortran_libraries/UTILS_src/MERGE/r_interpolate_marged_sph.f90 b/src/Fortran_libraries/UTILS_src/MERGE/r_interpolate_marged_sph.f90
index 8d1ade4..6d23da1 100644
--- a/src/Fortran_libraries/UTILS_src/MERGE/r_interpolate_marged_sph.f90
+++ b/src/Fortran_libraries/UTILS_src/MERGE/r_interpolate_marged_sph.f90
@@ -12,6 +12,8 @@
!! subroutine itp_rj_merged_phys_from_IO(nnod_org, &
!! & nri_org, jmax_org, idx_gl_1d_j_org, d_rj_IO)
!! subroutine extend_potential_magne
+!! subroutine extend_inner_core_temp
+!! subroutine extend_inner_core_light
!!@endverbatim
!!
!!@param nnod_org Number of spectr data for original data
@@ -48,7 +50,7 @@
!
private :: nri_old2new, k_old2new_in, k_old2new_out
private :: coef_old2new_in
- private :: allocate_radial_itp_tbl
+ private :: allocate_radial_itp_tbl, extend_inner_core_scalar
!
! -----------------------------------------------------------------------
!
@@ -149,7 +151,7 @@
exit
end if
end do
- kr_outer_domain = nidx_rj(1) + 1
+ kr_outer_domain = nidx_rj(1)
do k = 1, nidx_rj(1)
if(radius_1d_rj_r(k) .gt. r_org(nri_org)) then
kr_outer_domain = k - 1
@@ -160,9 +162,9 @@
write(*,*) 'kr_inner_domain', kr_inner_domain
write(*,*) 'kr_outer_domain', kr_outer_domain
! do k = 1, nidx_rj(1)
-! write(*,'(i5,1pe16.8,2i5,1p3e16.8)') k, radius_1d_rj_r(k), &
-! & k_old2new_in(k), k_old2new_out(k), &
-! & r_org(k_old2new_in(k)), r_org(k_old2new_out(k)), &
+! write(*,'(i5,1pe16.8,2i5,1p3e16.8)') k, radius_1d_rj_r(k), &
+! & k_old2new_in(k), k_old2new_out(k), &
+! & r_org(k_old2new_in(k)), r_org(k_old2new_out(k)), &
! & coef_old2new_in(k)
! end do
!
@@ -183,26 +185,27 @@
integer(kind = kint) :: inod_gl, inod_in, inod_out
!
!
-!$omp parallel
+!!$omp parallel private(nd)
do nd = 1, ntot_phys_rj
-!$omp do private(j,j_gl,kr,inod_in,inod_out,inod_gl)
+!!$omp do private(j,j_gl,kr,inod_in,inod_out,inod_gl)
do j = 1, jmax_org
j_gl = idx_gl_1d_j_org(j,1)
!
- if(j_gl .ge. nidx_rj(2)) cycle
+ if(j_gl .lt. nidx_rj(2)) then
!
- do kr = kr_inner_domain, kr_outer_domain
- inod_gl = 1 + j_gl + (kr - 1) * nidx_rj(2)
- inod_in = j + (k_old2new_in(kr) - 1) * jmax_org
- inod_out = j + (k_old2new_out(kr) - 1) * jmax_org
- d_rj(inod_gl,nd) &
- & = coef_old2new_in(kr) * d_rj_IO(inod_in,nd) &
+ do kr = kr_inner_domain, kr_outer_domain
+ inod_gl = 1 + j_gl + (kr - 1) * nidx_rj(2)
+ inod_in = j + (k_old2new_in(kr) - 1) * jmax_org
+ inod_out = j + (k_old2new_out(kr) - 1) * jmax_org
+ d_rj(inod_gl,nd) &
+ & = coef_old2new_in(kr) * d_rj_IO(inod_in,nd) &
& + (1.0d0 - coef_old2new_in(kr)) * d_rj_IO(inod_out,nd)
- end do
+ end do
+ end if
end do
-!$omp end do
+!!$omp end do
end do
-!$omp end parallel
+!!$omp end parallel
!
end subroutine itp_rj_merged_phys_from_IO
!
@@ -239,4 +242,55 @@
!
! -----------------------------------------------------------------------
!
+ subroutine extend_inner_core_temp
+!
+ use m_phys_labels
+!
+!
+ call extend_inner_core_scalar(fhd_temp)
+!
+ end subroutine extend_inner_core_temp
+!
+! -----------------------------------------------------------------------
+!
+ subroutine extend_inner_core_light
+!
+ use m_phys_labels
+!
+!
+ call extend_inner_core_scalar(fhd_light)
+!
+ end subroutine extend_inner_core_light
+!
+! -----------------------------------------------------------------------
+!
+ subroutine extend_inner_core_scalar(field_name)
+!
+ use extend_potential_field
+!
+ use m_sph_spectr_data
+!
+ character(len = kchara), intent(in) :: field_name
+!
+ integer(kind = kint) :: is_field
+ integer(kind = kint) :: i
+!
+!
+ is_field = 0
+ do i = 1, num_phys_rj
+ if(phys_name_rj(i) .eq. field_name) then
+ is_field = istack_phys_comp_rj(i-1) + 1
+ exit
+ end if
+ end do
+ if(is_field .eq. 0) return
+!
+ if(kr_inner_domain .gt. 1) then
+ call ext_inside_scalar(is_field, kr_inner_domain)
+ end if
+!
+ end subroutine extend_inner_core_scalar
+!
+! -----------------------------------------------------------------------
+!
end module r_interpolate_marged_sph
diff --git a/src/programs/data_utilities/INITIAL_FIELD/const_sph_initial_spectr.f90 b/src/programs/data_utilities/INITIAL_FIELD/const_sph_initial_spectr.f90
index 5a4915a..7fe6fba 100644
--- a/src/programs/data_utilities/INITIAL_FIELD/const_sph_initial_spectr.f90
+++ b/src/programs/data_utilities/INITIAL_FIELD/const_sph_initial_spectr.f90
@@ -78,16 +78,16 @@
time = time_init
!
! Set initial velocity if velocity is exist
-! if(ipol%i_velo .gt. izero) call set_initial_velocity
+ if(ipol%i_velo .gt. izero) call set_initial_velocity
!
! Set initial temperature if temperature is exist
-! if(ipol%i_temp .gt. izero) call set_initial_temperature
+ if(ipol%i_temp .gt. izero) call set_initial_temperature
!
! Set initial composition if composition is exist
-! if(ipol%i_light .gt. izero) call set_initial_composition
+ if(ipol%i_light .gt. izero) call set_initial_composition
!
! Set initial magnetic field if magnetic field is exist
-! if(ipol%i_magne .gt. izero) call set_initial_magne_sph
+ if(ipol%i_magne .gt. izero) call set_initial_magne_sph
!
! Set heat source if heat source is exist
if(ipol%i_heat_source .gt. izero) then
@@ -333,8 +333,8 @@
ii = local_sph_data_address(k,jj)
rr = radius_1d_rj_r(k)
! Substitute initial heat source
- d_rj(ii,ipol%i_heat_source) = 0.35 * four*r_CMB**2 &
- & / (four * r_ICB**3 / three)
+ d_rj(ii,ipol%i_heat_source) = 0.893 * r_CMB**2 &
+ & / (r_ICB**3 / three)
end do
end if
!
More information about the CIG-COMMITS
mailing list