From e4a019393466159eaa0bb4e68e83bf6a61a12c94 Mon Sep 17 00:00:00 2001 From: "Brendan C. Lyons" Date: Mon, 9 Sep 2019 20:45:53 -0400 Subject: [PATCH 1/9] Split define_transport_coefficients 1. define_pellet_source before kprad_ionize - Advances pellet location - Calculates volume of pellet distribution for normalization - Reads LP source from file - Calculate ne and temperatures - Calculates ne and Te at pellet - Calculates ablation rate 2. define_transport_coefficients & define_sources before ludefall - Transport coefficients (only resistivity, kappa, viscosity, and visc_e) - Sources (density, momentum, & heat that used to be in define_transport_coefficients) 3. derived_quantities after kprad_advect - Calls define_transport_coefficients again - Calculates other derived quantities - Calculates scalar diagnostics (not ablation rate anymore) --- unstructured/adapt.f90 | 4 - unstructured/diagnostics.f90 | 114 +------- unstructured/hdf5_output.f90 | 1 - unstructured/makefile | 10 +- unstructured/newpar.f90 | 25 +- unstructured/pellet.f90 | 47 ++++ unstructured/time_step.f90 | 19 +- unstructured/time_step_split.f90 | 17 +- unstructured/transport.f90 | 435 ++++++++++++++++++++++--------- 9 files changed, 397 insertions(+), 275 deletions(-) diff --git a/unstructured/adapt.f90 b/unstructured/adapt.f90 index 348e7cd73..e956448ee 100644 --- a/unstructured/adapt.f90 +++ b/unstructured/adapt.f90 @@ -233,8 +233,6 @@ 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 derived_quantities(1) !ke_previous = ekin end subroutine adapt_by_psi @@ -406,8 +404,6 @@ 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) end if diff --git a/unstructured/diagnostics.f90 b/unstructured/diagnostics.f90 index 0b592f141..80e60e8da 100644 --- a/unstructured/diagnostics.f90 +++ b/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 @@ -669,7 +646,6 @@ subroutine calculate_scalars() use boundary_conditions use math use gyroviscosity - use pellet implicit none @@ -755,13 +731,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) @@ -911,25 +883,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() @@ -1018,8 +973,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 @@ -1065,73 +1018,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 ! ~~~~~~~~~~~~~~~ diff --git a/unstructured/hdf5_output.f90 b/unstructured/hdf5_output.f90 index ad5f4d5dd..7dff3d1f2 100644 --- a/unstructured/hdf5_output.f90 +++ b/unstructured/hdf5_output.f90 @@ -862,7 +862,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 diff --git a/unstructured/makefile b/unstructured/makefile index 5da2ee3ea..753288bde 100644 --- a/unstructured/makefile +++ b/unstructured/makefile @@ -143,12 +143,11 @@ OBJS := fftw_fortran.o read_namelist.o gsl_wrapper.o \ iterdb.o read_gyro.o read_neo.o radiation.o \ element.o $(V_OBJ) field.o nintegrate_mod.o \ M3Dmodules.o resistive_wall.o kprad.o \ - m3dc1_nint.o boundary.o gyroviscosity.o bootstrap.o \ + m3dc1_nint.o boundary.o gyroviscosity.o bootstrap.o \ metricterms_new.o two_fluid.o harned_mikic.o biharmonic.o \ - electric_field.o pellet.o parallel_heat_flux.o \ - temperature_plots.o \ + electric_field.o parallel_heat_flux.o temperature_plots.o \ electrostatic_potential.o newvar.o runaway.o diagnostics.o \ - read_schaffer_field.o neutral_beam.o \ + pellet.o read_schaffer_field.o neutral_beam.o \ coils.o coil_sets.o kprad_m3dc1.o model.o \ fit_magnetics.o wall.o \ init_common.o \ @@ -160,8 +159,7 @@ OBJS := fftw_fortran.o read_namelist.o gsl_wrapper.o \ init_eqdsk.o init_dskbal.o init_jsolver.o \ init_3dwave.o init_3ddiffusion.o init_frs.o init_ftz.o init_eigen.o \ init_intkink.o init_lz.o init_kstar.o init_basicq.o init_cyl.o \ - transport.o \ - auxiliary_fields.o \ + auxiliary_fields.o transport.o \ time_step_split.o time_step_unsplit.o \ hdf5_output.o particle.o \ time_step.o output.o restart_hdf5.o \ diff --git a/unstructured/newpar.f90 b/unstructured/newpar.f90 index 1b2755502..870d69042 100644 --- a/unstructured/newpar.f90 +++ b/unstructured/newpar.f90 @@ -275,8 +275,6 @@ Program Reducedquintic if(ntime.eq.0 .or. (ntime.eq.ntime0 .and. eqsubtract.eq.1)) then if(eqsubtract.eq.1) then - if(myrank.eq.0 .and. iprint.ge.2) print *, " transport coefficients" - call define_transport_coefficients call derived_quantities(0) if(iwrite_aux_vars.eq.1) call calculate_auxiliary_fields(0) end if @@ -295,8 +293,6 @@ Program Reducedquintic ! Calculate all quantities derived from basic fields ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - if(myrank.eq.0 .and. iprint.ge.2) print *, " transport coefficients" - call define_transport_coefficients call derived_quantities(1) @@ -677,8 +673,7 @@ end subroutine smooth_fields ! ====================================================================== ! derived_quantities ! ~~~~~~~~~~~~~~~~~~ -! calculates all derived quantities, including auxiliary fields -! and scalars +! calculates all derived quantities, including scalars ! ====================================================================== subroutine derived_quantities(ilin) use basic @@ -687,7 +682,6 @@ subroutine derived_quantities(ilin) use diagnostics use sparse use transport_coefficients - use auxiliary_fields implicit none @@ -698,6 +692,10 @@ subroutine derived_quantities(ilin) vectype :: temp + ! Density and temperature are now updated within define_transport_coefficients + if(myrank.eq.0 .and. iprint.ge.2) print *, " transport coefficients" + call define_transport_coefficients(ilin) + ! Find lcfs ! ~~~~~~~~~ if(myrank.eq.0 .and. iprint.ge.2) print *, " finding lcfs" @@ -738,25 +736,14 @@ subroutine derived_quantities(ilin) endif else if(myrank.eq.0 .and. iprint.ge.1) then - write(*,'(A,2e12.4)') ' no temperatue maximum found near ',xmag,zmag + write(*,'(A,2e12.4)') ' no temperature maximum found near ',xmag,zmag endif endif - ! Electron density - call calculate_ne(ilin, den_field(ilin), ne_field(ilin), eqsubtract) - - ! Define auxiliary fields ! ~~~~~~~~~~~~~~~~~~~~~~~ if(myrank.eq.0 .and. itimer.eq.1) call second(tstart) - if(itemp.eq.0 .and. (numvar.eq.3 .or. ipres.gt.0) .and. imp_temp.eq.0) then - if(myrank.eq.0 .and. iprint.ge.2) print *, " temperatures" - call calculate_temperatures(ilin, te_field(ilin), ti_field(ilin), & - pe_field(ilin), p_field(ilin), ne_field(ilin), den_field(ilin), & - eqsubtract) - end if - ! toroidal current if(myrank.eq.0 .and. iprint.ge.2) print *, " toroidal current" if(inocurrent_tor.eq.1) then diff --git a/unstructured/pellet.f90 b/unstructured/pellet.f90 index 31d4035d8..eccd29adf 100644 --- a/unstructured/pellet.f90 +++ b/unstructured/pellet.f90 @@ -114,6 +114,53 @@ subroutine pellet_init() end subroutine pellet_init + subroutine calculate_Lor_vol() + + use basic + use mesh_mod + use m3dc1_nint + use math + use diagnostics + + implicit none + + include 'mpif.h' + + integer :: itri, numelms, ier + 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 + + vectype elemental function pellet_distribution(ip, r, phi, z, pres, inorm) use math use basic diff --git a/unstructured/time_step.f90 b/unstructured/time_step.f90 index 17d69414e..f851cda32 100644 --- a/unstructured/time_step.f90 +++ b/unstructured/time_step.f90 @@ -107,6 +107,11 @@ subroutine onestep calc_matrices = 0 endif + if(ipellet.ne.0) then + if(myrank.eq.0 .and. iprint.ge.2) print *, " define pellet source" + call define_pellet_source(1) + end if + ! Advance impurity charge states if(myrank.eq.0 .and. itimer.eq.1) call second(tstart) call kprad_ionize(dt) @@ -116,7 +121,10 @@ subroutine onestep endif if(myrank.eq.0 .and. iprint.ge.2) print *, " transport coefficients" - call define_transport_coefficients + call define_transport_coefficients(1) + + if(myrank.eq.0 .and. iprint.ge.2) print *, " sources" + call define_sources(1) ! start of loop to repeat timestep if max iterations exceeded in 3D do irepeat = 1, max_repeat @@ -177,8 +185,6 @@ subroutine onestep dtold = dt if(ntime.gt.1 .and. linear.eq.0) call variable_timestep - call pellet_advance - call runaway_advance ! copy time advance vectors to field data @@ -187,9 +193,6 @@ subroutine onestep ! if(eqsubtract.eq.0) call subtract_axi !DEBUG call export_time_advance_vectors - ! Calculate all quantities derived from basic fields - call derived_quantities(1) - ! Advect impurity charge states if(myrank.eq.0 .and. itimer.eq.1) call second(tstart) call kprad_advect(dt) @@ -198,6 +201,10 @@ subroutine onestep t_kprad = t_kprad + tend - tstart endif + ! Calculate all quantities derived from basic fields + call derived_quantities(1) + + ! Conserve toroidal flux if(iconstflux.eq.1 .and. numvar.ge.2) then call conserve_flux diff --git a/unstructured/time_step_split.f90 b/unstructured/time_step_split.f90 index 764e6d75e..4e93d08f7 100644 --- a/unstructured/time_step_split.f90 +++ b/unstructured/time_step_split.f90 @@ -864,16 +864,17 @@ subroutine step_split(calc_matrices) call destroy_vector(temp) call destroy_vector(temp2) - if(irecalc_eta.eq.1) then - call export_time_advance_vectors_split - call define_transport_coefficients - end if - - ! Electron temperature + ! Electron density neold_vec = ne_vec call calculate_ne(1, den_v, ne_v, eqsubtract) ne_field(1) = ne_v ! This is needed so that boundary_te works properly den_field(1) = den_v + + if(irecalc_eta.eq.1) then + call export_time_advance_vectors_split + call define_transport_coefficients(1) + end if + endif ! on idens=1 ! @@ -1103,7 +1104,7 @@ subroutine step_split(calc_matrices) call matvecmult(d2_mat,phi_vec,b2_phi) call add(b1_phi, b2_phi) - ! Inculde density terms + ! Include density terms if(idens.eq.1) then call matvecmult(r42_mat,ne_vec,b2_phi) call add(b1_phi, b2_phi) @@ -1204,7 +1205,7 @@ subroutine step_split(calc_matrices) call export_time_advance_vectors_split ! redefine transport coefficients with new den/pe values call lcfs(psi_field(1)) - call define_transport_coefficients + call define_transport_coefficients(1) ! revert fields to old values phi_vec = b2_phi diff --git a/unstructured/transport.f90 b/unstructured/transport.f90 index c3963c7a9..1a065be9e 100644 --- a/unstructured/transport.f90 +++ b/unstructured/transport.f90 @@ -1112,7 +1112,7 @@ end function bs_func ! define_transport_coefficients ! ============================= -subroutine define_transport_coefficients() +subroutine define_transport_coefficients(ilin) use basic use arrays @@ -1123,21 +1123,20 @@ subroutine define_transport_coefficients() use pellet use diagnostics use kprad_m3dc1 + use auxiliary_fields implicit none include 'mpif.h' + integer, intent(in) :: ilin integer :: itri, izone integer :: numelms, def_fields,ier logical, save :: first_time = .true. - logical :: solve_sigma, solve_kappa, solve_visc, solve_resistivity, & - solve_visc_e, solve_q, solve_totrad, solve_linerad, solve_bremrad, & - solve_ionrad, solve_reckrad, solve_recprad, solve_cd, solve_f, & - solve_fp + logical :: solve_kappa, solve_visc, solve_resistivity, solve_visc_e - integer, parameter :: num_scalars = 15 + integer, parameter :: num_scalars = 4 integer, dimension(num_scalars) :: temp, temp2 vectype, dimension(dofs_per_element) :: dofs @@ -1145,15 +1144,185 @@ subroutine define_transport_coefficients() if((linear.eq.1).and.(.not.first_time)) return first_time = .false. - if(myrank.eq.0 .and. iprint.ge.1) & - print *, "Calculating transport coefficients" + if(myrank.eq.0 .and. iprint.ge.1) print *, "Calculating transport coefficients" + + if(myrank.eq.0 .and. iprint.ge.2) print *, " Calculate electron density and temperatures" + call calculate_ne(ilin, den_field(ilin), ne_field(ilin), eqsubtract) + if(itemp.eq.0 .and. (numvar.eq.3 .or. ipres.gt.0) .and. imp_temp.eq.0) & + call calculate_temperatures(ilin, te_field(ilin), ti_field(ilin), & + pe_field(ilin), p_field(ilin), ne_field(ilin), & + den_field(ilin),eqsubtract) ! which transport coefficients need matrix solve solve_resistivity = .false. solve_visc = .false. solve_kappa = .false. - solve_sigma = .false. solve_visc_e = .false. + + ! clear variables + resistivity_field = 0. + kappa_field = 0. + visc_field = 0. + if(ibootstrap.ne.0) visc_e_field = 0. + + call finalize(field0_vec) + call finalize(field_vec) + + ! specify which primitive fields are to be evalulated + def_fields = FIELD_N + FIELD_PE + FIELD_P + FIELD_PSI + FIELD_I + FIELD_B2I + if(itemp.ge.1) def_fields = def_fields + FIELD_TE + if(iresfunc.eq.2 .or. iresfunc.eq.3 .or. iresfunc.eq.4) & + def_fields = def_fields + FIELD_ETA + if(ivisfunc.eq.3) def_fields = def_fields + FIELD_MU + + if(myrank.eq.0 .and. iprint.ge.2) print *, ' defining...' + + ! Calculate RHS + numelms = local_elements() +!$OMP PARALLEL DO & +!$OMP& PRIVATE(dofs) + do itri=1,numelms + + call define_element_quadrature(itri, int_pts_aux, 5) + call define_fields(itri, def_fields, 1, linear) + + call get_zone(itri, izone) + + + dofs = resistivity_func() + if(.not.solve_resistivity) solve_resistivity = any(dofs.ne.0.) + +!$OMP CRITICAL + if(solve_resistivity) & + call vector_insert_block(resistivity_field%vec,itri,1,dofs,VEC_ADD) +!$OMP END CRITICAL + + dofs = kappa_func() + if(.not.solve_kappa) solve_kappa = any(dofs.ne.0.) +!$OMP CRITICAL + if(solve_kappa) & + call vector_insert_block(kappa_field%vec,itri,1,dofs,VEC_ADD) +!$OMP END CRITICAL + + dofs = viscosity_func() + if(.not.solve_visc) solve_visc = any(dofs.ne.0.) +!$OMP CRITICAL + if(solve_visc) & + call vector_insert_block(visc_field%vec,itri,1,dofs,VEC_ADD) +!$OMP END CRITICAL + + if(ibootstrap.ne.0) then + dofs = electron_viscosity_func() + if(.not.solve_visc_e) solve_visc_e = any(dofs.ne.0.) +!$OMP CRITICAL + if(solve_visc_e) & + call vector_insert_block(visc_e_field%vec,itri,1,dofs,VEC_ADD) +!$OMP END CRITICAL + end if + end do +!$OMP END PARALLEL DO + + + ! Solve all the variables that have been defined + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! make sure all processes agree on what needs to be solved + if(maxrank.gt.1) then + temp = 0 + temp2 = 0 + if(solve_resistivity) temp(1) = 1 + if(solve_kappa) temp(2) = 1 + if(solve_visc) temp(3) = 1 + if(solve_visc_e) temp(4) = 1 + + call mpi_allreduce(temp, temp2, num_scalars, MPI_INTEGER, & + MPI_MAX, MPI_COMM_WORLD, ier) + + solve_resistivity = temp2(1).eq.1 + solve_kappa = temp2(2).eq.1 + solve_visc = temp2(3).eq.1 + solve_visc_e = temp2(4).eq.1 + end if + + if(myrank.eq.0 .and. iprint.ge.1) print *, ' solving...' + + if(solve_resistivity) then + if(myrank.eq.0 .and. iprint.ge.1) print *, ' resistivity' + call newvar_solve(resistivity_field%vec, mass_mat_lhs) + end if + + if(solve_kappa) then + if(myrank.eq.0 .and. iprint.ge.1) print *, ' kappa' + call newvar_solve(kappa_field%vec, mass_mat_lhs) + endif + + if(solve_visc) then + if(myrank.eq.0 .and. iprint.ge.1) print *, ' viscosity' + call newvar_solve(visc_field%vec, mass_mat_lhs) + endif + + if(solve_visc_e) then + if(myrank.eq.0 .and. iprint.ge.1) print *, ' electron viscosity' + call newvar_solve(visc_e_field%vec, mass_mat_lhs) + endif + + ! the "compressible" viscosity is the same as the "incompressible" + ! viscosity up to a constant + visc_c_field = visc_field + + ! add in constant components + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~ + call add(resistivity_field, etar*eta_fac) + call add(visc_field, amu) + call add(visc_c_field, amuc) + + if(myrank.eq.0 .and. iprint.ge.2) print *, 'done define_transport_coefficients' + +end subroutine define_transport_coefficients + + +! define_sources +! ============================= +subroutine define_sources(ilin) + + use basic + use arrays + use m3dc1_nint + use newvar_mod + use sparse + use neutral_beam + use pellet + use diagnostics + use kprad_m3dc1 + use auxiliary_fields + + implicit none + + include 'mpif.h' + + integer, intent(in) :: ilin + integer :: itri, izone + integer :: numelms, def_fields,ier + + logical :: solve_sigma, solve_q, solve_totrad, solve_linerad, solve_bremrad, & + solve_ionrad, solve_reckrad, solve_recprad, solve_cd, solve_f, solve_fp + + integer, parameter :: num_scalars = 11 + integer, dimension(num_scalars) :: temp, temp2 + vectype, dimension(dofs_per_element) :: dofs + + if(myrank.eq.0 .and. iprint.ge.1) print *, "Calculating sources" + + if(myrank.eq.0 .and. iprint.ge.2) print *, " Calculate electron density and temperatures" + call calculate_ne(ilin, den_field(ilin), ne_field(ilin), eqsubtract) + if(itemp.eq.0 .and. (numvar.eq.3 .or. ipres.gt.0) .and. imp_temp.eq.0) & + call calculate_temperatures(ilin, te_field(ilin), ti_field(ilin), & + pe_field(ilin), p_field(ilin), ne_field(ilin), & + den_field(ilin),eqsubtract) + + + ! which transport coefficients need matrix solve + solve_sigma = .false. solve_f = .false. solve_q = .false. solve_totrad = .false. @@ -1166,10 +1335,6 @@ subroutine define_transport_coefficients() solve_fp = .false. ! clear variables - resistivity_field = 0. - kappa_field = 0. - - visc_field = 0. if(density_source) sigma_field = 0. if(momentum_source) Fphi_field = 0. if(heat_source) Q_field = 0. @@ -1182,7 +1347,6 @@ subroutine define_transport_coefficients() Recprad_field = 0. end if if(icd_source .gt. 0) cd_field = 0. - if(ibootstrap.ne.0) visc_e_field = 0. if(ipforce.gt.0) pforce_field = 0. if(ipforce.gt.0) pmach_field = 0. @@ -1192,20 +1356,12 @@ subroutine define_transport_coefficients() ! specify which primitive fields are to be evalulated def_fields = FIELD_N + FIELD_PE + FIELD_P + FIELD_PSI + FIELD_I + FIELD_B2I if(itemp.ge.1) def_fields = def_fields + FIELD_TE - if(iresfunc.eq.2 .or. iresfunc.eq.3 .or. iresfunc.eq.4) & - def_fields = def_fields + FIELD_ETA - if(ivisfunc.eq.3) def_fields = def_fields + FIELD_MU if(ibeam.ge.1) def_fields = def_fields + FIELD_V if(ipforce.gt.0) def_fields = def_fields + FIELD_PHI + FIELD_CHI + FIELD_NI if(iarc_source.ne.0) def_fields = def_fields + FIELD_WALL if(myrank.eq.0 .and. iprint.ge.2) print *, ' defining...' - - if(ipellet.ne.0) then - ! make sure normalization for pellet_distribution defined - call calculate_Lor_vol - end if ! Calculate RHS numelms = local_elements() @@ -1218,23 +1374,6 @@ subroutine define_transport_coefficients() call get_zone(itri, izone) - - dofs = resistivity_func() - if(.not.solve_resistivity) solve_resistivity = any(dofs.ne.0.) - -!$OMP CRITICAL - if(solve_resistivity) & - call vector_insert_block(resistivity_field%vec,itri,1,dofs,VEC_ADD) -!$OMP END CRITICAL - - dofs = kappa_func() - if(.not.solve_kappa) solve_kappa = any(dofs.ne.0.) -!$OMP CRITICAL - if(solve_kappa) & - call vector_insert_block(kappa_field%vec,itri,1,dofs,VEC_ADD) -!$OMP END CRITICAL - - if(density_source) then dofs = sigma_func(izone) if(.not.solve_sigma) solve_sigma = any(dofs.ne.0.) @@ -1244,13 +1383,6 @@ subroutine define_transport_coefficients() !$OMP END CRITICAL end if - dofs = viscosity_func() - if(.not.solve_visc) solve_visc = any(dofs.ne.0.) -!$OMP CRITICAL - if(solve_visc) & - call vector_insert_block(visc_field%vec,itri,1,dofs,VEC_ADD) -!$OMP END CRITICAL - if(momentum_source) then dofs = force_func(izone) if(.not.solve_f) solve_f = any(dofs.ne.0.) @@ -1335,14 +1467,6 @@ subroutine define_transport_coefficients() call vector_insert_block(cd_field%vec,itri,1,dofs,VEC_ADD) end if - if(ibootstrap.ne.0) then - dofs = electron_viscosity_func() - if(.not.solve_visc_e) solve_visc_e = any(dofs.ne.0.) -!$OMP CRITICAL - if(solve_visc_e) & - call vector_insert_block(visc_e_field%vec,itri,1,dofs,VEC_ADD) -!$OMP END CRITICAL - end if end do !$OMP END PARALLEL DO @@ -1354,69 +1478,41 @@ subroutine define_transport_coefficients() if(maxrank.gt.1) then temp = 0 temp2 = 0 - if(solve_resistivity) temp(1) = 1 - if(solve_kappa) temp(2) = 1 - if(solve_sigma) temp(3) = 1 - if(solve_visc) temp(4) = 1 - if(solve_visc_e) temp(5) = 1 - if(solve_f) temp(6) = 1 - if(solve_q) temp(7) = 1 - if(solve_fp) temp(8) = 1 - if(solve_cd) temp(9) = 1 - if(solve_totrad) temp(10) = 1 - if(solve_linerad) temp(11) = 1 - if(solve_bremrad) temp(12) = 1 - if(solve_ionrad) temp(13) = 1 - if(solve_reckrad) temp(14) = 1 - if(solve_recprad) temp(15) = 1 + if(solve_sigma) temp(1) = 1 + if(solve_f) temp(2) = 1 + if(solve_q) temp(3) = 1 + if(solve_fp) temp(4) = 1 + if(solve_cd) temp(5) = 1 + if(solve_totrad) temp(6) = 1 + if(solve_linerad) temp(7) = 1 + if(solve_bremrad) temp(8) = 1 + if(solve_ionrad) temp(9) = 1 + if(solve_reckrad) temp(10) = 1 + if(solve_recprad) temp(11) = 1 call mpi_allreduce(temp, temp2, num_scalars, MPI_INTEGER, & MPI_MAX, MPI_COMM_WORLD, ier) - solve_resistivity = temp2(1).eq.1 - solve_kappa = temp2(2).eq.1 - solve_sigma = temp2(3).eq.1 - solve_visc = temp2(4).eq.1 - solve_visc_e = temp2(5).eq.1 - solve_f = temp2(6).eq.1 - solve_q = temp2(7).eq.1 - solve_fp = temp2(8).eq.1 - solve_cd = temp2(9).eq.1 - solve_totrad = temp2(10).eq.1 - solve_linerad = temp2(11).eq.1 - solve_bremrad = temp2(12).eq.1 - solve_ionrad = temp2(13).eq.1 - solve_reckrad = temp2(14).eq.1 - solve_recprad = temp2(15).eq.1 + solve_sigma = temp2(1).eq.1 + solve_f = temp2(2).eq.1 + solve_q = temp2(3).eq.1 + solve_fp = temp2(4).eq.1 + solve_cd = temp2(5).eq.1 + solve_totrad = temp2(6).eq.1 + solve_linerad = temp2(7).eq.1 + solve_bremrad = temp2(8).eq.1 + solve_ionrad = temp2(9).eq.1 + solve_reckrad = temp2(10).eq.1 + solve_recprad = temp2(11).eq.1 end if if(myrank.eq.0 .and. iprint.ge.1) print *, ' solving...' - if(solve_resistivity) then - if(myrank.eq.0 .and. iprint.ge.1) print *, ' resistivity' - call newvar_solve(resistivity_field%vec, mass_mat_lhs) - end if - - if(solve_kappa) then - if(myrank.eq.0 .and. iprint.ge.1) print *, ' kappa' - call newvar_solve(kappa_field%vec, mass_mat_lhs) - endif - - if(solve_sigma) then + if(solve_sigma) then if(myrank.eq.0 .and. iprint.ge.1) print *, ' sigma' call newvar_solve(sigma_field%vec, mass_mat_lhs_dc) endif - if(solve_visc) then - if(myrank.eq.0 .and. iprint.ge.1) print *, ' viscosity' - call newvar_solve(visc_field%vec, mass_mat_lhs) - endif - - if(solve_visc_e) then - if(myrank.eq.0 .and. iprint.ge.1) print *, ' electron viscosity' - call newvar_solve(visc_e_field%vec, mass_mat_lhs) - endif - if(solve_f) then if(myrank.eq.0 .and. iprint.ge.1) print *, ' fphi' call newvar_solve(Fphi_field%vec, mass_mat_lhs_dc) @@ -1465,27 +1561,54 @@ subroutine define_transport_coefficients() if(solve_fp) then if(myrank.eq.0 .and. iprint.ge.1) print *, ' pforce' call newvar_solve(pforce_field%vec, mass_mat_lhs) + endif - if(myrank.eq.0 .and. iprint.ge.1) print *, ' pmach' - call newvar_solve(pmach_field%vec, mass_mat_lhs) + if(myrank.eq.0 .and. iprint.ge.2) & + print *, 'done define_sources' - endif +end subroutine define_sources - ! the "compressible" viscosity is the same as the "incompressible" - ! viscosity up to a constant - visc_c_field = visc_field +! define pellet source +! ============================= +subroutine define_pellet_source(ilin) + + use basic + use mesh_mod + use arrays + use m3dc1_nint + use diagnostics + use metricterms_new + use kprad_m3dc1 + use pellet + use math + use auxiliary_fields + + implicit none + + include 'mpif.h' + + integer, intent(in) :: ilin ! 0 for equilibrium fields, 1 for perturbed + integer :: itri, numelms, def_fields, ier + integer :: ip + integer :: izone, izonedim + real :: tpifac, tpirzero + double precision, allocatable :: ptemp(:) + if(myrank.eq.0 .and. iprint.ge.1) & + print *, "begin define_pellet_source" - ! add in constant components - ! ~~~~~~~~~~~~~~~~~~~~~~~~~~ - call add(resistivity_field, etar*eta_fac) - call add(visc_field, amu) - call add(visc_c_field, amuc) + if(ntime.ne.0) then + if(myrank.eq.0 .and. iprint.ge.2) print *, " Advance pellet position" + call pellet_advance + end if + ! Define normalization for pellet_distribution + call calculate_Lor_vol ! Read LP data if(iread_lp_source.eq.1) then + if(myrank.eq.0 .and. iprint.ge.2) print *, " Read LP source" call read_lp_source('cloud.txt', ier) if(ier.ne.0) then if(myrank.eq.0) print *, 'Error reading LP source ', 'cloud.txt' @@ -1493,9 +1616,83 @@ subroutine define_transport_coefficients() end if end if - if(myrank.eq.0 .and. iprint.ge.2) & - print *, 'done define_transport_coefficients' + if(myrank.eq.0 .and. iprint.ge.2) print *, " Calculate electron density and temperatures" + call calculate_ne(ilin, den_field(ilin), ne_field(ilin), eqsubtract) + if(itemp.eq.0 .and. (numvar.eq.3 .or. ipres.gt.0) .and. imp_temp.eq.0) & + call calculate_temperatures(ilin, te_field(ilin), ti_field(ilin), & + pe_field(ilin), p_field(ilin), ne_field(ilin), & + den_field(ilin),eqsubtract) + + if(ipellet_abl.gt.0) then + if(myrank.eq.0 .and. iprint.ge.2) print *, " Calculate ablation inputs" + + call tpi_factors(tpifac,tpirzero) + nsource_pel = 0. + temp_pel = 0. + + def_fields = FIELD_P + FIELD_N + call finalize(field0_vec) + call finalize(field_vec) + numelms = local_elements() + +!$OMP PARALLEL DO PRIVATE(itri,ier,izone,izonedim,ip) & +!$OMP& REDUCTION(+:nsource_pel,temp_pel) + 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, def_fields, 0, 0) + + ! Pellet radius and density/temperature at the pellet surface + 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 + end do +!$OMP END PARALLEL DO + + 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) + + if(myrank.eq.0 .and. iprint.ge.2) print *, " Calculate ablation rate" + call calculate_ablation + + end if + + if(myrank.eq.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 *, " 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 + if(ipellet_abl.gt.0) then + 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 + end if + end do + end if + + if(myrank.eq.0 .and. iprint.ge.1) print *, 'done define_pellet_source' + +end subroutine define_pellet_source -end subroutine define_transport_coefficients end module transport_coefficients From 0181456b26b383a5968c851fba35a4d8efc35f37 Mon Sep 17 00:00:00 2001 From: "Brendan C. Lyons" Date: Tue, 10 Sep 2019 17:57:48 -0400 Subject: [PATCH 2/9] Fix pellet advance for first time step --- unstructured/init_conds.f90 | 1 - unstructured/transport.f90 | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/unstructured/init_conds.f90 b/unstructured/init_conds.f90 index b098dd725..2ce17b526 100644 --- a/unstructured/init_conds.f90 +++ b/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 diff --git a/unstructured/transport.f90 b/unstructured/transport.f90 index 1a065be9e..4a696e62b 100644 --- a/unstructured/transport.f90 +++ b/unstructured/transport.f90 @@ -1598,7 +1598,7 @@ subroutine define_pellet_source(ilin) if(myrank.eq.0 .and. iprint.ge.1) & print *, "begin define_pellet_source" - if(ntime.ne.0) then + if(ntime.gt.1) then if(myrank.eq.0 .and. iprint.ge.2) print *, " Advance pellet position" call pellet_advance end if From 9dc194b3ab7b448c4dc7102936eb688e4e9f4693 Mon Sep 17 00:00:00 2001 From: "Brendan C. Lyons" Date: Fri, 13 Sep 2019 17:33:37 -0400 Subject: [PATCH 3/9] Remove transport coefficients from derived quantities --- unstructured/adapt.f90 | 5 ++++- unstructured/newpar.f90 | 15 +++++++++++---- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/unstructured/adapt.f90 b/unstructured/adapt.f90 index e956448ee..41827ad78 100644 --- a/unstructured/adapt.f90 +++ b/unstructured/adapt.f90 @@ -233,6 +233,7 @@ subroutine adapt_by_psi i_control%err_p_old = 0. n_control%err_i = 0. n_control%err_p_old = 0. + call define_transport_coefficients(1) call derived_quantities(1) !ke_previous = ekin end subroutine adapt_by_psi @@ -405,8 +406,10 @@ subroutine adapt_by_error n_control%err_p_old = 0. call reset_scalars 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 diff --git a/unstructured/newpar.f90 b/unstructured/newpar.f90 index 870d69042..db24fac28 100644 --- a/unstructured/newpar.f90 +++ b/unstructured/newpar.f90 @@ -275,6 +275,7 @@ Program Reducedquintic if(ntime.eq.0 .or. (ntime.eq.ntime0 .and. eqsubtract.eq.1)) then if(eqsubtract.eq.1) then + call define_transport_coefficients(0) call derived_quantities(0) if(iwrite_aux_vars.eq.1) call calculate_auxiliary_fields(0) end if @@ -293,6 +294,7 @@ Program Reducedquintic ! Calculate all quantities derived from basic fields ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + call define_transport_coefficients(1) call derived_quantities(1) @@ -682,6 +684,7 @@ subroutine derived_quantities(ilin) use diagnostics use sparse use transport_coefficients + use auxiliary_fields implicit none @@ -692,10 +695,14 @@ subroutine derived_quantities(ilin) vectype :: temp - ! Density and temperature are now updated within define_transport_coefficients - if(myrank.eq.0 .and. iprint.ge.2) print *, " transport coefficients" - call define_transport_coefficients(ilin) - + ! Update density and temperature + if(myrank.eq.0 .and. iprint.ge.2) print *, " Calculate electron density and temperatures" + call calculate_ne(ilin, den_field(ilin), ne_field(ilin), eqsubtract) + if(itemp.eq.0 .and. (numvar.eq.3 .or. ipres.gt.0) .and. imp_temp.eq.0) & + call calculate_temperatures(ilin, te_field(ilin), ti_field(ilin), & + pe_field(ilin), p_field(ilin), ne_field(ilin), & + den_field(ilin),eqsubtract) + ! Find lcfs ! ~~~~~~~~~ if(myrank.eq.0 .and. iprint.ge.2) print *, " finding lcfs" From ec2e8a8a256925359cc0cc541eaa69d33cb3a25b Mon Sep 17 00:00:00 2001 From: "Brendan C. Lyons" Date: Tue, 1 Oct 2019 01:24:00 -0400 Subject: [PATCH 4/9] Move ikappafunc=5 definition to define_fields --- unstructured/m3dc1_nint.f90 | 104 +++++++++++++++++++++++++++++++++++- unstructured/transport.f90 | 22 +++----- 2 files changed, 109 insertions(+), 17 deletions(-) diff --git a/unstructured/m3dc1_nint.f90 b/unstructured/m3dc1_nint.f90 index 7cdfe4f84..2352908a9 100644 --- a/unstructured/m3dc1_nint.f90 +++ b/unstructured/m3dc1_nint.f90 @@ -493,7 +493,16 @@ subroutine define_fields(itri, fieldi, gdef, ilin, ieqs) end if if(iand(fields, FIELD_ETA).eq.FIELD_ETA) then if(iresfunc.eq.2 .or. iresfunc.eq.3) fields = ior(fields,FIELD_PSI) - if(iresfunc.eq.4) fields = ior(ior(fields,FIELD_N),FIELD_P) + if(iresfunc.eq.4) then + fields = ior(ior(fields,FIELD_N),FIELD_P) + if(itemp.eq.1) fields = ior(fields,FIELD_TE) + end if + end if + if(iand(fields, FIELD_KAP).eq.FIELD_KAP) then + if(ikappafunc.eq.5) then + fields = ior(ior(ior(fields,FIELD_N),FIELD_P),FIELD_PSI) + if(itemp.eq.1) fields = ior(fields,FIELD_TE) + end if end if if(iand(fields, FIELD_MU).eq.FIELD_MU) then if(ivisfunc.eq.3) fields = ior(fields,FIELD_PSI) @@ -1098,7 +1107,98 @@ subroutine define_fields(itri, fieldi, gdef, ilin, ieqs) if(iand(fields, FIELD_KAP).eq.FIELD_KAP) then if(itri.eq.1 .and. myrank.eq.0 .and. iprint.ge.2) print *, " kappa..." - call eval_ops(itri, kappa_field, kap79) + if(ikappafunc.eq.5) then + kap79 = 0. + if(izone.eq.1. .and. kappa0.gt.0) then + ! Te + if(itemp.eq.1) then + temp79b = tet79(:,OP_1) + else + temp79b = pet79(:,OP_1)/net79(:,OP_1) + end if + + kap79 = 0. + where(real(temp79b).gt.(kappa0/(kappa_max-kappat))) + kap79(:,OP_1) = net79(:,OP_1)/pet79(:,OP_1) + kap79(:,OP_DR) = net79(:,OP_DR)/pet79(:,OP_1) & + - net79(:,OP_1)*pet79(:,OP_DR)/pet79(:,OP_1)**2 + kap79(:,OP_DZ) = net79(:,OP_DZ)/pet79(:,OP_1) & + - net79(:,OP_1)*pet79(:,OP_DZ)/pet79(:,OP_1)**2 + kap79(:,OP_DRR) = net79(:,OP_DRR)/pet79(:,OP_1) & + - 2.*net79(:,OP_DR)*pet79(:,OP_DR)/pet79(:,OP_1)**2 & + + 2.*net79(:,OP_1)*pet79(:,OP_DR)**2/pet79(:,OP_1)**3 & + - net79(:,OP_1)*pet79(:,OP_DRR)/pet79(:,OP_1)**2 + kap79(:,OP_DRZ) = net79(:,OP_DRZ)/pet79(:,OP_1) & + - (net79(:,OP_DR)*pet79(:,OP_DZ) + & + net79(:,OP_DZ)*pet79(:,OP_DR))/pet79(:,OP_1)**2 & + + 2.*net79(:,OP_1)*pet79(:,OP_DR)*pet79(:,OP_DZ)/pet79(:,OP_1)**3 & + - net79(:,OP_1)*pet79(:,OP_DRZ)/pet79(:,OP_1)**2 + kap79(:,OP_DZZ) = net79(:,OP_DZZ)/pet79(:,OP_1) & + - 2.*net79(:,OP_DZ)*pet79(:,OP_DZ)/pet79(:,OP_1)**2 & + + 2.*net79(:,OP_1)*pet79(:,OP_DZ)**2/pet79(:,OP_1)**3 & + - net79(:,OP_1)*pet79(:,OP_DZZ)/pet79(:,OP_1)**2 +#ifdef USE3D + kap79(:,OP_DP) = net79(:,OP_DP)/pet79(:,OP_1) & + - net79(:,OP_1)*pet79(:,OP_DP)/pet79(:,OP_1)**2 +#endif + end where + + kap79 = kap79*kappa0 + kap79(:,OP_1) = kap79(:,OP_1) + kappat + + where(real(temp79b).le.(kappa0/(kappa_max-kappat))) + kap79(:,OP_1) = kappa_max + end where + + if(kappaf.ge.0. .and. gradp_crit.ne.0) then + temp79a = pt79(:,OP_DR)**2 + pt79(:,OP_DZ)**2 +#ifdef USE3D + temp79a = temp79a + ri2_79*pt79(:,OP_DP)**2 +#endif + do i=1, OP_NUM + where(real(temp79a).lt.gradp_crit**2) kap79(:,i) = kap79(:,i) * kappaf + end do + end if + + if(kappah.ne.0.) then + tm79 = (pst79 - psimin)/(psibound - psimin) + temp79a = tanh((real(tm79(:,OP_1))-1.)/.2) + temp79b = cosh((real(tm79(:,OP_1))-1.)/.2)**-1 + kap79(:,OP_1) = kap79(:,OP_1) + kappah*temp79a**2 + kap79(:,OP_DR) = kap79(:,OP_DR) & + + kappah*temp79a*temp79b**2*tm79(:,OP_DR) + kap79(:,OP_DZ) = kap79(:,OP_DZ) & + + kappah*temp79a*temp79b**2*tm79(:,OP_DZ) + kap79(:,OP_DRR) = kap79(:,OP_DRR) & + + kappah*(temp79a*temp79b**2*tm79(:,OP_DRR) & + + 0.5*(temp79b**4)*(tm79(:,OP_DR)**2) & + - (temp79a**2)*(temp79b**2)*(tm79(:,OP_DR)**2)) + kap79(:,OP_DRZ) = kap79(:,OP_DRZ) & + + kappah*(temp79a*temp79b**2*tm79(:,OP_DRZ) & + + 0.5*(temp79b**4)*(tm79(:,OP_DR)**2) & + - (temp79a**2)*(temp79b**2)*(tm79(:,OP_DR)**2)) + kap79(:,OP_DZZ) = kap79(:,OP_DZZ) & + + kappah*(temp79a*temp79b**2*tm79(:,OP_DZZ) & + + 0.5*(temp79b**4)*tm79(:,OP_DR)*tm79(:,OP_DZ) & + - (temp79a**2)*(temp79b**2)*tm79(:,OP_DR)*tm79(:,OP_DZ)) +#ifdef USE3D + kap79(:,OP_DP) = kap79(:,OP_DP) & + + kappah*temp79a*temp79b**2*tm79(:,OP_DP) +#endif + end if + + where(kap79.ne.kap79) kap79 = 0. + temp79a = kap79(:,OP_1) + do i=1, OP_NUM + where(real(temp79a).lt.0.) kap79(:,i) = 0. + where(real(temp79a).gt.kappa_max) kap79(:,i) = 0. + end do + where(real(temp79a).gt.kappa_max) kap79(:,OP_1) = kappa_max + + end if + else + call eval_ops(itri, kappa_field, kap79) + end if if(ikapscale.eq.1) then kar79 = kappar*kap79 diff --git a/unstructured/transport.f90 b/unstructured/transport.f90 index 4a696e62b..4f2c25b3b 100644 --- a/unstructured/transport.f90 +++ b/unstructured/transport.f90 @@ -939,19 +939,9 @@ function kappa_func() end if case(5) - ! kappa ~ 1/|Te| + ! kappa ~ 1/Te (with a maximum) - if(kappa0.eq.0) then - temp79a = 0. - else - if(itemp.eq.0) then - temp79a = kappa0*abs(real(net79(:,OP_1))/real(pet79(:,OP_1))) - else - temp79a = kappa0/abs(real(tet79(:,OP_1))) - end if - end if - - where(real(temp79a).gt.(kappa_max-kappat)) temp79a = kappa_max - kappat + temp79a = kap79(:,OP_1) - kappat case(10,11) if(.not.allocated(kappa_spline%x)) then @@ -997,7 +987,8 @@ function kappa_func() temp79a = temp79a + kappat - if(kappaf.ge.0. .and. gradp_crit.ne.0) then + ! BCL 9/30/19: ikappafunc condition added here since defined in m3dc1_nint.f90 + if(kappaf.ge.0. .and. gradp_crit.ne.0 .and. ikappafunc.ne.5) then temp79b = pt79(:,OP_DR)**2 + pt79(:,OP_DZ)**2 #ifdef USE3D temp79b = temp79b + ri2_79*pt79(:,OP_DP)**2 @@ -1009,12 +1000,13 @@ function kappa_func() temp = temp + intx2(mu79(:,:,OP_1),temp79a) - if(kappah.ne.0.) then + ! BCL 9/30/19: ikappafunc condition added here since defined in m3dc1_nint.f90 + if(kappah.ne.0. .and. ikappafunc.ne.5) then temp79b = (pst79(:,OP_1) - psimin)/(psibound - psimin) temp79a = kappah*tanh((real(temp79b) - 1.)/.2)**2 temp = temp + intx2(mu79(:,:,OP_1),temp79a) end if - + kappa_func = temp end function kappa_func From 706af29d2c1c9f1f7487eca46e17ff8b368e5c8c Mon Sep 17 00:00:00 2001 From: "Brendan C. Lyons" Date: Fri, 4 Oct 2019 20:08:43 -0400 Subject: [PATCH 5/9] Update volume integral to work in 3D runs --- unstructured/idl/mask_bound.pro | 26 ------------------------- unstructured/idl/volume_integral.pro | 29 ++++++++++++++++++++++------ 2 files changed, 23 insertions(+), 32 deletions(-) delete mode 100644 unstructured/idl/mask_bound.pro diff --git a/unstructured/idl/mask_bound.pro b/unstructured/idl/mask_bound.pro deleted file mode 100644 index 8fe3dec22..000000000 --- a/unstructured/idl/mask_bound.pro +++ /dev/null @@ -1,26 +0,0 @@ -; Given R and Z matrices and boundary points B (from get_boundary_path) -; return a mask matrix for R,Z points inside the boundary - -function mask_bound, R, Z, B - - sr = size(R) - mask = fltarr(sr[1], sr[2], sr[3])+1. - - sb = size(B) - nb = sb[2] - - for k=0,nb-1 do begin - - - bx2 = B[0,(k+1) mod nb]-R - by2 = B[1,(k+1) mod nb]-Z - bx1 = B[0,k]-R - by1 = B[1,k]-Z - - mask = mask and ((bx2*by1 - bx1*by2) le 0.) - - endfor - - return, mask - -end \ No newline at end of file diff --git a/unstructured/idl/volume_integral.pro b/unstructured/idl/volume_integral.pro index 52e8bcac7..72324cd8a 100644 --- a/unstructured/idl/volume_integral.pro +++ b/unstructured/idl/volume_integral.pro @@ -1,23 +1,40 @@ ; Compute the toroidal volume integral of field inside the boundary ; /core: volume inside LCFS instead of inside boundary -function volume_integral, field, x, y, t, slice=slice, filename=filename, points=pts, core=core,_EXTRA=ex +function volume_integral, name, x, y, t, slice=slice, filename=filename, points=pts, core=core,_EXTRA=ex - if(size(field, /type) eq 7) then out = read_field(field,x,y,t,slices=slice,filename=filename,points=pts,_EXTRA=ex) + i3d = read_parameter('3d', filename=filename) + nplanes = read_parameter('nplanes', filename=filename) + if(size(name, /type) eq 7) then begin + if(i3d eq 0) then begin + field = read_field(name,x,y,t,slices=slice,filename=filename,points=pts,_EXTRA=ex) + endif else begin + field = read_field(name,x,y,t,slices=slice,filename=filename,points=pts,taverage=5*nplanes,_EXTRA=ex) + endelse + endif else begin + field = name + endelse R = radius_matrix(x,y,t) Z = z_matrix(x,y,t) if keyword_set(core) then begin - psin = read_field('psi_norm',x,y,t,slices=slice,filename=filename,points=pts,_EXTRA=ex) + if(i3d eq 0) then begin + psin = read_field('psi_norm',x,y,t,slices=slice,filename=filename,points=pts,_EXTRA=ex) + endif else begin + psin = read_field('psi_norm',x,y,t,slices=slice,filename=filename,points=pts,taverage=5*nplanes,_EXTRA=ex) + endelse lcfs = find_lcfs(psi,x,y,xpoint=xp,filename=filename,points=pts,_EXTRA=ex) mask = psin le 1.0 ; exclude private flux region if (xp[1] lt 0.) then mask = mask and (Z ge xp[1]) else mask = mask and (Z le xp[1]) endif else begin B = get_boundary_path(filename=filename,points=pts) - mask = mask_bound(R,Z,B) + P = obj_new('IDLanROI',B[0,*],B[1,*]) + mask = P->ContainsPoints(R,Z) + mask = reform(mask,[1,pts,pts]) gt 0 + R = reform(R,[1,pts,pts]) + Z = reform(Z,[1,pts,pts]) endelse - ; assumes uniform x and uniform y grids - return, 2.0*!pi*(x[1]-x[0])*(y[1]-y[0])*total(mask*R*out,/nan) + return, 2.0*!pi*(x[1]-x[0])*(y[1]-y[0])*total(mask*R*field,/nan) end From a029bb1d66f71e0dd25e1487817775cf3ef785fd Mon Sep 17 00:00:00 2001 From: "Brendan C. Lyons" Date: Fri, 4 Oct 2019 20:10:25 -0400 Subject: [PATCH 6/9] Bug fix for toroidal average in read_field.pro --- unstructured/idl/read_field.pro | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unstructured/idl/read_field.pro b/unstructured/idl/read_field.pro index 087f2e60c..14f49eec0 100644 --- a/unstructured/idl/read_field.pro +++ b/unstructured/idl/read_field.pro @@ -33,7 +33,7 @@ function read_field, name, x, y, t, slices=slices, mesh=mesh, $ data = 0 if(taverage eq 1) then taverage=16 phi = period*findgen(taverage) / taverage - if(itor eq 1) then period = period*180./!pi + if(itor eq 1) then phi = phi*180./!pi for i=0, taverage-1 do begin data = data + $ read_field(name, x, y, t, slices=time, mesh=mesh, $ From d44b098599c67c078fb92c7ddb83c5a30134200a Mon Sep 17 00:00:00 2001 From: "Brendan C. Lyons" Date: Thu, 10 Oct 2019 15:29:46 -0400 Subject: [PATCH 7/9] Fix merging error for ikappafunc=5 --- unstructured/transport.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/unstructured/transport.f90 b/unstructured/transport.f90 index 4f2c25b3b..3cbf0f9ef 100644 --- a/unstructured/transport.f90 +++ b/unstructured/transport.f90 @@ -1165,6 +1165,7 @@ subroutine define_transport_coefficients(ilin) if(itemp.ge.1) def_fields = def_fields + FIELD_TE if(iresfunc.eq.2 .or. iresfunc.eq.3 .or. iresfunc.eq.4) & def_fields = def_fields + FIELD_ETA + if(ikappafunc.eq.5) def_fields = def_fields + FIELD_KAP if(ivisfunc.eq.3) def_fields = def_fields + FIELD_MU if(myrank.eq.0 .and. iprint.ge.2) print *, ' defining...' From 4d945906890ceedca3f05f0ac50b2880a7246db1 Mon Sep 17 00:00:00 2001 From: blyons_nersc Date: Wed, 16 Oct 2019 03:41:21 -0700 Subject: [PATCH 8/9] Prevent output_scalars from shrinking datasets --- unstructured/hdf5_output.f90 | 44 ++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 12 deletions(-) diff --git a/unstructured/hdf5_output.f90 b/unstructured/hdf5_output.f90 index 7dff3d1f2..313554eac 100644 --- a/unstructured/hdf5_output.f90 +++ b/unstructured/hdf5_output.f90 @@ -541,11 +541,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 @@ -555,12 +557,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) @@ -573,20 +575,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) @@ -597,6 +614,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 From 675ef151b50d20e16ba4d9ba57f62dce14f0a86e Mon Sep 17 00:00:00 2001 From: blyons_nersc Date: Fri, 6 Dec 2019 16:10:15 -0800 Subject: [PATCH 9/9] Move tpi_factors to basic to prevent circular dependency --- unstructured/M3Dmodules.f90 | 27 +++++++++++++++++++++++++++ unstructured/diagnostics.f90 | 25 ------------------------- unstructured/pellet.f90 | 7 +++---- unstructured/transport.f90 | 1 - 4 files changed, 30 insertions(+), 30 deletions(-) diff --git a/unstructured/M3Dmodules.f90 b/unstructured/M3Dmodules.f90 index 993e911de..6a7bdd9c7 100644 --- a/unstructured/M3Dmodules.f90 +++ b/unstructured/M3Dmodules.f90 @@ -458,6 +458,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 diff --git a/unstructured/diagnostics.f90 b/unstructured/diagnostics.f90 index 881faf1a8..5464578e6 100644 --- a/unstructured/diagnostics.f90 +++ b/unstructured/diagnostics.f90 @@ -603,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 ! ----------------- diff --git a/unstructured/pellet.f90 b/unstructured/pellet.f90 index ab1ca0acd..4189a6438 100644 --- a/unstructured/pellet.f90 +++ b/unstructured/pellet.f90 @@ -42,7 +42,7 @@ module pellet subroutine pellet_init() use basic use read_ascii -! use diagnostics + implicit none character(LEN=10), parameter :: pellet_filename = 'pellet.dat' @@ -120,7 +120,6 @@ subroutine calculate_Lor_vol() use mesh_mod use m3dc1_nint use math - use diagnostics implicit none @@ -164,7 +163,7 @@ end subroutine calculate_Lor_vol vectype elemental function pellet_distribution(ip, r, phi, z, pres, inorm) use math use basic -! use diagnostics + implicit none integer, intent(in) :: ip real, intent(in) :: r, phi, z, pres @@ -262,7 +261,7 @@ end function pellet_distribution subroutine pellet_advance use basic -! use diagnostics + implicit none real, allocatable :: x(:), y(:) diff --git a/unstructured/transport.f90 b/unstructured/transport.f90 index aa5c97a6d..b1c728882 100644 --- a/unstructured/transport.f90 +++ b/unstructured/transport.f90 @@ -1574,7 +1574,6 @@ subroutine define_pellet_source(ilin) use mesh_mod use arrays use m3dc1_nint - use diagnostics use metricterms_new use kprad_m3dc1 use pellet