Skip to content

Fix memory and initialization error in RE simulation with method of characteristics #90

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
changliu777 wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Fix memory and initialization error in RE simulation with method of characteristics #90

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions unstructured/electric_field.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ subroutine electric_field_R(ilin,o,izone)

! eta J
! ~~~~~
if (.not.((kinetic.eq.1).or.((irunaway.ge.1).and.(ra_characteristics.eq.1)))) then
if(ilin.eq.1) then
o = o - temp79b*bz179(:,OP_DZ)*ri_79
#if defined(USE3D) || defined(USECOMPLEX)
Expand All @@ -35,6 +36,7 @@ subroutine electric_field_R(ilin,o,izone)
(ri2_79*pst79(:,OP_DRP) - ri_79*bfpt79(:,OP_DZP))
#endif
end if
endif
if(izone.eq.ZONE_CONDUCTOR) return

! -VxB
Expand Down Expand Up @@ -150,6 +152,7 @@ subroutine electric_field_Z(ilin,o,izone)

! eta J
! ~~~~~
if (.not.((kinetic.eq.1).or.((irunaway.ge.1).and.(ra_characteristics.eq.1)))) then
if(ilin.eq.1) then
o = o + temp79b*bz179(:,OP_DR)*ri_79
#if defined(USE3D) || defined(USECOMPLEX)
Expand All @@ -163,6 +166,7 @@ subroutine electric_field_Z(ilin,o,izone)
(ri2_79*pst79(:,OP_DZP) + ri_79*bfpt79(:,OP_DRP))
#endif
end if
endif
if(izone.eq.ZONE_CONDUCTOR) return

! -VxB
Expand Down Expand Up @@ -273,11 +277,13 @@ subroutine electric_field_phi(ilin,o, izone)

! eta J
! ~~~~~
if (.not.((kinetic.eq.1).or.((irunaway.ge.1).and.(ra_characteristics.eq.1)))) then
if(ilin.eq.1) then
o = o - ri_79*eta79(:,OP_1)*ps179(:,OP_GS)
else
o = o - ri_79*eta79(:,OP_1)*pst79(:,OP_GS)
end if
endif
if(izone.eq.ZONE_CONDUCTOR) return

! VxB
Expand Down
14 changes: 11 additions & 3 deletions unstructured/hdf5_output.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -196,13 +196,21 @@ subroutine read_real_attr(parent_id, name, value, error)
character(LEN=*), intent(in) :: name
real, intent(out) :: value
integer, intent(out) :: error
logical :: attr_exists

integer(HID_T) :: attr_id
integer(HSIZE_T), dimension(1) :: dims = 1

call h5aopen_name_f(parent_id, name, attr_id, error)
call h5aread_f(attr_id, H5T_NATIVE_DOUBLE, value, dims, error)
call h5aclose_f(attr_id, error)
! call h5aopen_name_f(parent_id, name, attr_id, error)
! call h5aread_f(attr_id, H5T_NATIVE_DOUBLE, value, dims, error)
! call h5aclose_f(attr_id, error)

call h5aexists_f(parent_id, name, attr_exists, error)
if (error == 0 .and. attr_exists) then
call h5aopen_name_f(parent_id, name, attr_id, error)
call h5aread_f(attr_id, H5T_NATIVE_DOUBLE, value, dims, error)
call h5aclose_f(attr_id, error)
end if

end subroutine read_real_attr

Expand Down
1 change: 1 addition & 0 deletions unstructured/output.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -316,6 +316,7 @@ subroutine hdf5_write_parameters(error)
call write_real_attr(root_id, "zlim2" , zlim2, error)
call write_int_attr (root_id, "ikprad" , ikprad, error)
call write_int_attr (root_id, "kprad_z" , kprad_z, error)
call write_real_attr(root_id, "bzsign" , bzsign, error)

call h5gclose_f(root_id, error)

Expand Down
2 changes: 2 additions & 0 deletions unstructured/restart_hdf5.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,6 +136,8 @@ subroutine rdrestart_hdf5()
ikprad_in = 0
end if

if ((irunaway.gt.0).or.(kinetic.eq.1)) call read_real_attr(root_id, "bzsign", bzsign, error)

! Read Scalars
call h5gopen_f(root_id, "scalars", scalar_group_id, error)

Expand Down
21 changes: 15 additions & 6 deletions unstructured/runaway_advection.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -450,7 +450,7 @@ subroutine init_particles(lrestart, ierr)
call get_field_coefs(1)
!call MPI_Win_fence(0, win_elfieldcoefs)

npar=nelms_global*npoints*1.1/nrows
npar=nelms_global*npoints*2.1/nrows

!Allocate local storage for particle data
disp_unit = 1
Expand Down Expand Up @@ -571,6 +571,12 @@ subroutine init_particles(lrestart, ierr)
if (gid_max<locparts2) gid_max=locparts2
locparts = locparts + 1
pdata(locparts2) = dpar
else
locparts2 = (itri-1)*npoints+ipar
dpar%gid = locparts2
locparts2 = locparts2-(ielm_min-1)*npoints
dpar%deleted = .true.
pdata(locparts2) = dpar
endif
end do !iz
end do
Expand All @@ -585,13 +591,16 @@ subroutine init_particles(lrestart, ierr)
call mpi_allreduce(gid_min, ipart_begin, 1, MPI_INTEGER, MPI_MIN, hostcomm, ierr)
call mpi_allreduce(gid_max, ipart_end, 1, MPI_INTEGER, MPI_MAX, hostcomm, ierr)
call mpi_barrier(mpi_comm_world, ierr)
particle_map(:,:)=1
if (hostrank == 0) then
write(0,*) ipart_begin, ipart_end
do ipart = ipart_begin, ipart_end
itri=int((pdata(ipart)%gid-1)/npoints)+1
ipar=mod((pdata(ipart)%gid-1),npoints)+1
! write(0,*) ipar,itri,ipart
particle_map(ipar,itri)=ipart
if (.not.pdata(ipart)%deleted) then
itri=int((pdata(ipart)%gid-1)/npoints)+1
ipar=mod((pdata(ipart)%gid-1),npoints)+1
! write(0,*) ipar,itri,ipart
particle_map(ipar,itri)=ipart
endif
end do
end if
call mpi_barrier(mpi_comm_world, ierr)
Expand All @@ -611,7 +620,7 @@ subroutine init_particles(lrestart, ierr)
#endif

call create_field(nre_vec)
call set_matrix_index(diff2_mat, 75)
call set_matrix_index(diff2_mat, 175)
call create_mat(diff2_mat, 1, 1, icomplex, 1)
!call create_newvar_matrix(diff_mat, NV_NOBOUND, NV_I_MATRIX, 1)
do itri=1,nelms
Expand Down