WIP: Split define_transport_coefficients

unstructured/M3Dmodules.f90

@@ -459,6 +459,33 @@ module basic
 
   type(spline1d) :: q_spline
 
+contains
+
+  !   Added 1/1/2016 to get consistency between 2D,3D,Cyl,Tor
+  !   BCL moved 12/6/2019
+  subroutine tpi_factors(tpifac,tpirzero)
+    use math
+    implicit none
+    real, intent(out) :: tpifac, tpirzero
+    if(nplanes.eq.1) then
+       if(itor.eq.1) then
+          tpifac = 1.
+          tpirzero = 1.
+       else
+          tpifac = 1./rzero
+          tpirzero = 1.
+       endif
+    else
+       if(itor.eq.1) then
+          tpifac = twopi
+          tpirzero = twopi
+       else
+          tpifac = twopi
+          tpirzero = twopi*rzero
+       endif
+    endif
+  end subroutine tpi_factors
+
 end module basic
 
 module arrays

unstructured/adapt.f90

@@ -233,8 +233,7 @@ subroutine adapt_by_psi
     i_control%err_p_old = 0.
     n_control%err_i = 0.
     n_control%err_p_old = 0.
-    if(myrank.eq.0 .and. iprint.ge.2) print *, "  transport coefficients"
-    call define_transport_coefficients
+    call define_transport_coefficients(1)
     call derived_quantities(1)
     !ke_previous = ekin
   end subroutine adapt_by_psi
@@ -406,11 +405,11 @@ subroutine adapt_by_error
        n_control%err_i = 0.
        n_control%err_p_old = 0.
        call reset_scalars
-       if(myrank.eq.0 .and. iprint.ge.2) print *, "  transport coefficients"
-       call define_transport_coefficients
        if(eqsubtract.eq.1) then
-         call derived_quantities(0)
+          call define_transport_coefficients(0)
+          call derived_quantities(0)
        end if
+       call define_transport_coefficients(1)
        call derived_quantities(1)
        meshAdapted =1
     end if

unstructured/diagnostics.f90

@@ -80,9 +80,6 @@ module diagnostics
   real, dimension(iflux_loops_max) :: flux_loop_val
   integer, dimension(iflux_loops_max) :: flux_loop_itri
 
-  vectype, dimension(MAX_PTS) :: rhop79
-  vectype, dimension(MAX_PTS) :: Lorentz_pel
-
 contains
 
   ! ======================================================================
@@ -164,7 +161,6 @@ end subroutine distribute_timings
   ! resets diagnostic energy and scalar quantities to zero
   ! ======================================================================
   subroutine reset_scalars()
-    use pellet
 
     implicit none
 
@@ -235,10 +231,6 @@ subroutine reset_scalars()
     nfluxd = 0.
     nfluxv = 0.
     nsource = 0.
-    if(ipellet_abl.gt.0) then
-       nsource_pel = 0. ! this is an array
-       temp_pel = 0.    ! this is an array
-    end if
 
     bwb2 = 0.
 
@@ -272,7 +264,6 @@ end subroutine reset_scalars
   ! ======================================================================
   subroutine distribute_scalars()    
     use basic
-    use pellet
 
     implicit none
 
@@ -281,7 +272,6 @@ subroutine distribute_scalars()
     integer, parameter :: num_scalars = 73
     integer :: ier
     double precision, dimension(num_scalars) :: temp, temp2
-    double precision, allocatable  :: ptemp(:)
 
     ! Allreduce energy terms
     if(maxrank .gt. 1) then
@@ -437,19 +427,6 @@ subroutine distribute_scalars()
        wall_force_n0_z_halo = temp2(72)
        helicity        = temp2(73)
 
-       if(ipellet_abl.gt.0) then
-          allocate(ptemp(npellets))
-
-          ptemp = nsource_pel
-          call mpi_allreduce(ptemp, nsource_pel, npellets, MPI_DOUBLE_PRECISION,  &
-                             MPI_SUM, MPI_COMM_WORLD, ier)
-          ptemp = temp_pel
-          call mpi_allreduce(ptemp, temp_pel, npellets, MPI_DOUBLE_PRECISION,  &
-                            MPI_SUM, MPI_COMM_WORLD, ier)
-
-          deallocate(ptemp)
-       end if
-
     endif
 
   end subroutine distribute_scalars
@@ -626,31 +603,6 @@ subroutine second(tcpu)
   end subroutine second
 
 
-!   Added 1/1/2016 to get consistency between 2D,3D,Cyl,Tor
-subroutine tpi_factors(tpifac,tpirzero)
-  use basic
-  use math
-  implicit none
-  real, intent(out) :: tpifac, tpirzero
-  if(nplanes.eq.1) then
-     if(itor.eq.1) then
-        tpifac = 1.
-        tpirzero = 1.
-     else
-        tpifac = 1./rzero
-        tpirzero = 1.
-     endif
-  else
-     if(itor.eq.1) then
-        tpifac = twopi
-        tpirzero = twopi
-     else
-        tpifac = twopi
-        tpirzero = twopi*rzero
-     endif
-  endif
-end subroutine tpi_factors
-
 ! ======================================================================
 ! calculate scalars
 ! -----------------
@@ -669,7 +621,6 @@ subroutine calculate_scalars()
   use boundary_conditions
   use math
   use gyroviscosity
-  use pellet
 
   implicit none
 
@@ -755,13 +706,9 @@ subroutine calculate_scalars()
   call finalize(field_vec)
 
   numelms = local_elements()
-
-  if(ipellet.ne.0) call calculate_Lor_vol
 
-  ! BCL Warning: nsource_pel and temp_pel are now vectors
-  !              this compiles, but may break at runtime for OpenMP (OMP=1)
 !$OMP PARALLEL DO PRIVATE(mr,dum1,ier,is_edge,n,iedge,idim,izone,izonedim,i) &
-!$OMP& REDUCTION(+:ekinp,ekinpd,ekinph,ekint,ekintd,ekinth,ekin3,ekin3d,ekin3h,wallcur,emagp,emagpd,emagph,emagt,emagtd,emagth,emag3,area,parea,totcur,pcur,m_iz,tflux,pflux,tvor,volume,pvol,totden,pden,totrad,linerad,bremrad,ionrad,reckrad,recprad,totre,nsource,epotg,tmom,pmom,bwb2,efluxp,efluxt,efluxs,efluxk,tau_em,tau_sol,tau_com,tau_visc,tau_gyro,tau_parvisc,nfluxd,nfluxv,xray_signal,Lor_vol,nsource_pel,temp_pel,wall_force_n0_x,wall_force_n0_y,wall_force_n0_z,wall_force_n1_x,wall_force_n1_y,wall_force_n1_z,totne,w_pe,pcur_co,pcur_sn,m_iz_co,m_iz_sn,w_m,w_p,wall_force_n0_x_halo,wall_force_n0_z_halo,helicity)
+!$OMP& REDUCTION(+:ekinp,ekinpd,ekinph,ekint,ekintd,ekinth,ekin3,ekin3d,ekin3h,wallcur,emagp,emagpd,emagph,emagt,emagtd,emagth,emag3,area,parea,totcur,pcur,m_iz,tflux,pflux,tvor,volume,pvol,totden,pden,totrad,linerad,bremrad,ionrad,reckrad,recprad,totre,nsource,epotg,tmom,pmom,bwb2,efluxp,efluxt,efluxs,efluxk,tau_em,tau_sol,tau_com,tau_visc,tau_gyro,tau_parvisc,nfluxd,nfluxv,xray_signal,wall_force_n0_x,wall_force_n0_y,wall_force_n0_z,wall_force_n1_x,wall_force_n1_y,wall_force_n1_z,totne,w_pe,pcur_co,pcur_sn,m_iz_co,m_iz_sn,w_m,w_p,wall_force_n0_x_halo,wall_force_n0_z_halo,helicity)
   do itri=1,numelms
 
      !call zonfac(itri, izone, izonedim)
@@ -912,25 +859,8 @@ subroutine calculate_scalars()
 #endif
 
      ! particle source
-     if(idens.eq.1) then        
-        nsource = nsource - twopi*int1(sig79)/tpifac
+     if(idens.eq.1) nsource = nsource - twopi*int1(sig79)/tpifac
 
-        ! Pellet radius and density/temperature at the pellet surface
-        if(ipellet_abl.gt.0) then
-           do ip=1,npellets
-              if(r_p(ip).ge.1e-8) then
-                 ! weight density/temp by pellet distribution (normalized)
-                 temp79a = pellet_distribution(ip, x_79, phi_79, z_79, real(pt79(:,OP_1)), 1)
-                 nsource_pel(ip) = nsource_pel(ip) + twopi*int2(net79(:,OP_1),temp79a)/tpifac
-                 temp_pel(ip) = temp_pel(ip) + twopi*int2(pet79(:,OP_1)/net79(:,OP_1),temp79a)*p0_norm/(1.6022e-12*n0_norm*tpifac)
-              else
-                 nsource_pel(ip) = 0.
-                 temp_pel(ip) = 0.
-              end if
-           end do
-       endif
-     endif
-
      ! gravitational potential energy
      epotg = epotg + grav_pot()
 
@@ -1020,8 +950,6 @@ subroutine calculate_scalars()
 
   call distribute_scalars
 
-  if(ipellet_abl.gt.0) call calculate_ablation
-
   ekin = ekinp + ekint + ekin3
   emag = emagp + emagt + emag3
   ekind = ekinpd + ekintd + ekin3d
@@ -1067,73 +995,10 @@ subroutine calculate_scalars()
      print *, "  Ionization loss = ", ionrad
      print *, "  Recombination radiation (kinetic) = ", reckrad
      print *, "  Recombination radiation (potential) = ", recprad
-     if(ipellet_abl.gt.0 .and. iprint.ge.2) then
-        do ip=1,npellets
-           print *, "  Pellet #", ip
-           print *, "    particles injected = ",pellet_rate(ip)*dt*(n0_norm*l0_norm**3)
-           print *, "    radius (in cm) = ", r_p(ip)*l0_norm
-           print *, "    local electron temperature (in eV) = ", temp_pel(ip)
-           print *, "    local electron density (in ne14) = ", nsource_pel(ip)
-           print *, "    rpdot (in cm/s) = ", rpdot(ip)*l0_norm/t0_norm
-           print *, "    Lor_vol = ", Lor_vol(ip)
-           print *, "    R position: ", pellet_r(ip)*l0_norm
-           print *, "    phi position: ", pellet_phi(ip)
-           print *, "    Z position: ", pellet_z(ip)*l0_norm
-        end do
-     endif
   endif
 
 end subroutine calculate_scalars
 
-
-subroutine calculate_Lor_vol()
-
-  use basic
-  use mesh_mod
-  use m3dc1_nint
-  use math
-  use pellet
-
-  implicit none
-
-  include 'mpif.h'
-
-  integer :: itri, numelms, ier
-  integer :: is_edge(3)  ! is inode on boundary
-  real :: tpifac,tpirzero
-  integer :: izone, izonedim
-  real, allocatable :: temp(:)
-  integer :: ip
-
-  call tpi_factors(tpifac,tpirzero)
-
-  numelms = local_elements()
-
-  allocate(temp(npellets))
-  temp = 0.
-
-  do itri=1,numelms
-
-     call m3dc1_ent_getgeomclass(2, itri-1,izonedim,izone)
-     if(izone.ne.1) cycle
-     call define_element_quadrature(itri, int_pts_diag, int_pts_tor)
-     call define_fields(itri, FIELD_P, 0, 0)
-
-     ! perform volume integral of pellet cloud (without normalization)
-     do ip=1,npellets
-        temp79a  = pellet_distribution(ip, x_79, phi_79, z_79, real(pt79(:,OP_1)), 0)
-        temp(ip) = temp(ip) + twopi*int1(temp79a)/tpifac
-     end do
-
-  end do
-
-  call mpi_allreduce(temp, Lor_vol, npellets, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, ier )
-
-  deallocate(temp)
-
-end subroutine calculate_Lor_vol
-
-
 !======================================================================
 ! magnetic_region
 ! ~~~~~~~~~~~~~~~

unstructured/hdf5_output.f90

@@ -552,11 +552,13 @@ subroutine output_scalar(parent_id, name, value, t, error)
     integer(HSIZE_T) :: chunk_size(1) = (/ 100 /)
     integer(HSIZE_T) :: dims(1)
     integer(HSIZE_T) :: maxdims(1)
-    integer(HSIZE_T), parameter :: local_dims(1) = (/ 1 /)
-    integer(HSIZE_T), dimension(1,1) :: coord
-    integer(SIZE_T), parameter :: num_elements = 1
+    integer(HSIZE_T) :: local_dims(1)
+    integer(HSIZE_T), allocatable :: coord(:,:)
+    integer(SIZE_T)  :: num_elements
     integer(HID_T) :: memspace, filespace, dset_id, p_id, plist_id
-    real :: values(1)
+    integer(HSIZE_T) :: curdim
+    integer :: i
+    real, allocatable :: values(:)
     logical :: exists
 
 #ifdef USETAU
@@ -566,12 +568,12 @@ subroutine output_scalar(parent_id, name, value, t, error)
 
     dims(1) = t+1
     maxdims(1) = H5S_UNLIMITED_F
-    values(1) = value
-    coord(1,1) = t + 1
-
+
+    !@if(myrank.eq.0 .and. iprint.ge.2) print *, "   ", name
     call h5lexists_f(parent_id, name, exists, error)
 
     if(.not.exists) then
+       num_elements = 1
        call h5screate_simple_f(1, dims, filespace, error, maxdims)
        call h5pcreate_f(H5P_DATASET_CREATE_F, p_id, error)
        call h5pset_chunk_f(p_id, 1, chunk_size, error)
@@ -584,20 +586,35 @@ subroutine output_scalar(parent_id, name, value, t, error)
        end if
        call h5pclose_f(p_id, error)
        call h5sclose_f(filespace, error)
+       curdim = dims(1)
     else
        call h5dopen_f(parent_id, name, dset_id, error)
-       call h5dset_extent_f(dset_id, dims, error)
+       call h5dget_space_f(dset_id, filespace, error)
+       call h5sget_simple_extent_npoints_f(filespace, curdim, error)
+       call h5sclose_f(filespace, error)
+       if(dims(1).gt.curdim) then
+          num_elements = 1
+          call h5dset_extent_f(dset_id, dims, error)
+          curdim = dims(1)
+       else
+          num_elements = curdim - dims(1) + 1
+       end if
     endif
 
+    local_dims(1) = num_elements
+    allocate(values(num_elements))
+    allocate(coord(1,num_elements))
+    values = sqrt(-1.)
+    values(1) = value
+    coord(1,:) = (/ (i, i = dims(1), curdim) /)
+
     if(myrank.eq.0) then
        call h5screate_simple_f(1, local_dims, memspace, error)
        call h5dget_space_f(dset_id, filespace, error)
-       call h5sselect_elements_f(filespace, H5S_SELECT_SET_F, 1, &
-            num_elements, coord, error)
-
        call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error)
        call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_INDEPENDENT_F, error)
-
+       call h5sselect_elements_f(filespace, H5S_SELECT_SET_F, 1, &
+            num_elements, coord, error)
        call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, values, local_dims, error, &
             file_space_id=filespace, mem_space_id=memspace, xfer_prp=plist_id)
 
@@ -608,6 +625,9 @@ subroutine output_scalar(parent_id, name, value, t, error)
     endif
     call h5dclose_f(dset_id, error)
 
+    deallocate(values)
+    deallocate(coord)
+
   end subroutine output_scalar
 
   ! read_scalar
@@ -873,7 +893,6 @@ subroutine read_1dextendarr(parent_id, name, value, NMAX, t, error)
     integer(HSIZE_T) :: dims(2), maxdims(2), local_dims(2), off(2)
     integer(SIZE_T) :: num_elements
     integer(HID_T) :: memspace, filespace, dset_id, p_id, plist_id
-    logical :: exists
 
 #ifdef USETAU
     integer :: dummy     ! this is necessary to prevent TAU from

unstructured/init_conds.f90

@@ -514,7 +514,6 @@ subroutine initial_conditions()
   use init_common
   use kprad_m3dc1
   use pellet
-  use diagnostics
   use cylinder
 
   implicit none