diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md index f6b3f67a..fe7105ff 100644 --- a/.github/ISSUE_TEMPLATE/bug_report.md +++ b/.github/ISSUE_TEMPLATE/bug_report.md @@ -2,7 +2,7 @@ name: Bug report about: Create a report to help us improve title: '' -labels: '' +labels: bug assignees: '' --- diff --git a/.github/ISSUE_TEMPLATE/code-enhancement.md b/.github/ISSUE_TEMPLATE/code-enhancement.md new file mode 100644 index 00000000..47d34f74 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/code-enhancement.md @@ -0,0 +1,20 @@ +--- +name: Code enhancement +about: Suggest an improvement for this project +title: '' +labels: enhancement +assignees: '' + +--- + +**Is your feature request related to a problem? Please describe.** +A clear and concise description of what the problem is. Ex. I'm always frustrated when [...] + +**Describe the solution you'd like** +A clear and concise description of what you want to happen. + +**Describe alternatives you've considered** +A clear and concise description of any alternative solutions or features you've considered. + +**Additional context** +Add any other context or screenshots about the feature request here. diff --git a/.gitignore b/.gitignore index 32f8cd5b..0465e71b 100644 --- a/.gitignore +++ b/.gitignore @@ -15,5 +15,4 @@ uedge.egg-info/ dist/uedge-8.0.0-py3.7-macosx-10.9-x86_64.egg uedge.egg-info *.c -*.f - +*.f \ No newline at end of file diff --git a/a.out b/a.out new file mode 100755 index 00000000..1c8e93d1 Binary files /dev/null and b/a.out differ diff --git a/aph/aph.v b/aph/aph.v index c6960b02..207e62ae 100755 --- a/aph/aph.v +++ b/aph/aph.v @@ -1,7 +1,7 @@ aph # Atomic physics (hydrogenic) ***** Physical_constants: -ev real /1.6022e-19/ # 1 electron volt in Joules +ev_aph real /1.6022e-19/ # 1 electron volt in Joules m_prot real /1.67e-27/ # proton mass ***** Data_input: @@ -15,7 +15,7 @@ aphdir character*120 # name of directory containing data files data_directory character*120 # another dirname containing data files. This is to be be passed in ***** Ionization_energy: -erad real [eV] /25./ # tot elec engy loss/ioniz (rad+binding) if istabon=0 +erad real [eV] /25./ +input # tot elec engy loss/ioniz (rad+binding) if istabon=0 ***** Rtdata: # hydrogenic rate table data from ADPAK via Braams' rate code @@ -145,23 +145,23 @@ readehr2(fname:string) subroutine ***** Aphwrk: # working arrays for 2-d spline interpolation -nxdata integer -nydata integer -xdata(1:nxdata) _real -ydata(1:nydata) _real -fdata(1:nxdata,1:nydata) _real -ldf integer -iflag integer -kxords integer /4/ # order of spline fit versus log(te) - # kxords=4 (default) is cubic interpolation -kyords integer /4/ # order of spline fit versus log10(ne) - # kyords=4 (default) is cubic interpolation -xknots(1:nxdata+kxords) _real -yknots(1:nydata+kyords) _real -workh(1:nxdata*nydata+2*kxords*(nxdata+1)) _real # work array -rsacoef(1:nxdata,1:nydata) _real # spline coeff's for ionization -rracoef(1:nxdata,1:nydata) _real # spline coeff's for recombination -rqacoef(1:nxdata,1:nydata) _real # spline coeff's for line emission +nxdata_aph integer +nydata_aph integer +xdata_aph(1:nxdata_aph) _real +ydata_aph(1:nydata_aph) _real +fdata_aph(1:nxdata_aph,1:nydata_aph) _real +ldf_aph integer +iflag_aph integer +kxords_aph integer /4/ # order of spline fit versus log(te) + # kxords_aph=4 (default) is cubic interpolation +kyords_aph integer /4/ # order of spline fit versus log10(ne) + # kyords_aph=4 (default) is cubic interpolation +xknots_aph(1:nxdata_aph+kxords_aph) _real +yknots_aph(1:nydata_aph+kyords_aph) _real +workh(1:nxdata_aph*nydata_aph+2*kxords_aph*(nxdata_aph+1)) _real # work array +rsacoef(1:nxdata_aph,1:nydata_aph) _real # spline coeff's for ionization +rracoef(1:nxdata_aph,1:nydata_aph) _real # spline coeff's for recombination +rqacoef(1:nxdata_aph,1:nydata_aph) _real # spline coeff's for line emission ***** Subs: # Subroutines that can be called from the parser diff --git a/aph/aphrates.m b/aph/aphrates.m index 33b581a2..abd76397 100755 --- a/aph/aphrates.m +++ b/aph/aphrates.m @@ -27,7 +27,7 @@ real function erl1 (te, ne, tau) c----------------------------------------------------------------------c if (istabon .le. 7) then # various older models - erl1 = (rqa(te,ne,0)-13.6*ev*rsa(te,ne,0.,0))*ne + erl1 = (rqa(te,ne,0)-13.6*ev_aph*rsa(te,ne,0.,0))*ne c----------------------------------------------------------------------c elseif (istabon .eq. 8 .or. istabon .eq. 9) then # linear interpolation @@ -35,7 +35,7 @@ real function erl1 (te, ne, tau) jr = 1 c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -63,7 +63,7 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- jr = 1 c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -98,7 +98,7 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- endif c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -174,7 +174,7 @@ real function erl2 (te, ne, tau) c----------------------------------------------------------------------c if (istabon .le. 7) then # various older models - erl2 = (13.6*ev+1.5*te)*ne*rra(te,ne,0.,1) + erl2 = (13.6*ev_aph+1.5*te)*ne*rra(te,ne,0.,1) c----------------------------------------------------------------------c elseif (istabon .eq. 8 .or. istabon .eq. 9) then # linear interpolation @@ -182,7 +182,7 @@ real function erl2 (te, ne, tau) jr = 1 c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -210,7 +210,7 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- jr = 1 c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -245,7 +245,7 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- endif c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -329,7 +329,7 @@ real function rcx (t0, n0, k) iti = 0 ini = 0 c compute abscissae -- - rlti = log(t0/ev) + rlti = log(t0/ev_aph) rlni = max(htln(0),min(htln(htnn),log(n0))) c find iti -- 51 if (iti .lt. htnt-1) then @@ -370,7 +370,7 @@ real function rcx (t0, n0, k) c use DEGAS table look-up c compute abscissa -- - zloge=log(t0/ev) + zloge=log(t0/ev_aph) rle=max(rlemin, min(zloge,rlemax)) c table index for interpolation -- je=int((rle-rlemin)/delekpt) + 1 @@ -394,7 +394,7 @@ call mcrates(n0,t0,t0,za,zamax,zn,kdum,kdum,rcxcopy) c----------------------------------------------------------------------c else # use analytic model (hydrogen) for all other istabon - a = 3*t0 / (10*ev) + a = 3*t0 / (10*ev_aph) rcx = 1.7e-14 * a**0.333 if (issgvcxc.eq.1) rcx = sgvcxc # use fixed sig-v if (issgvcxc.eq.2) rcx = sgvcxc*sqrt(t0/m_prot) # fixed sig @@ -443,8 +443,8 @@ real function rqa (te, ne, k) c----------------------------------------------------------------------c if (istabon .eq. 0) then # use analytic model (hydrogen) with c # constant energy loss per ionization - a = te / (10*ev) - rqa = erad * ev * 3.0e-14 * a*a / (3.0 + a*a) + a = te / (10*ev_aph) + rqa = erad * ev_aph * 3.0e-14 * a*a / (3.0 + a*a) c----------------------------------------------------------------------c elseif ((istabon .eq. 1) .or. (istabon .eq. 2)) then @@ -454,7 +454,7 @@ real function rqa (te, ne, k) ite = 0 ine = 0 c compute abscissae -- - rlte = log(te/ev) + rlte = log(te/ev_aph) rlne = max(htln(0),min(htln(htnn),log(ne))) c find ite -- 51 if (ite .lt. htnt-1) then @@ -488,7 +488,7 @@ real function rqa (te, ne, k) c compute electron energy loss rate parameter -- t0 = (1-fxne)*htlqa(ite,ine,k) + fxne*htlqa(ite,ine+1,k) t1 = (1-fxne)*htlqa(ite+1,ine,k) + fxne*htlqa(ite+1,ine+1,k) - rqa = ev*exp((1-fxte)*t0+fxte*t1) + rqa = ev_aph*exp((1-fxte)*t0+fxte*t1) c----------------------------------------------------------------------c elseif (istabon .eq. 3) then @@ -497,7 +497,7 @@ real function rqa (te, ne, k) jr = 1 c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -518,14 +518,14 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- w2=w21+fjd*(w22-w21) w = w1 + fje*(w2-w1) c electron energy loss rate parameter -- - rqa = ev * w * rsa(te,ne,0.,k) + rqa = ev_aph * w * rsa(te,ne,0.,k) c----------------------------------------------------------------------c elseif (istabon .eq. 4) then c use POST93 table look-up c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -546,16 +546,16 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- w2=w21+fjd*(w22-w21) w = w1 + fje*(w2-w1) c electron energy loss rate parameter -- - rqa = w + 13.6 * ev * rsa(te,ne,0.,k) + rqa = w + 13.6 * ev_aph * rsa(te,ne,0.,k) c----------------------------------------------------------------------c elseif (istabon .eq. 5) then c use spline fit to POST93 table data -c xuse=min(max(xdata(1),log(te/ev)),xdata(nxdata)) -c yuse=min(max(ydata(1),log10(ne)),ydata(nydata)) +c xuse=min(max(xdata_aph(1),log(te/ev_aph)),xdata_aph(nxdata_aph)) +c yuse=min(max(ydata_aph(1),log10(ne)),ydata_aph(nydata_aph)) c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -563,23 +563,23 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- xuse=rle yuse=rld - nxcoef=nxdata - nycoef=nydata + nxcoef=nxdata_aph + nycoef=nydata_aph - vlogw = B2VAhL(xuse, yuse, 0, 0, xknots, yknots, nxcoef, - . nycoef, kxords, kyords, rqacoef, ldf, workh, iflag) + vlogw = B2VAhL(xuse, yuse, 0, 0, xknots_aph, yknots_aph, nxcoef, + . nycoef, kxords_aph, kyords_aph, rqacoef, ldf_aph, workh, iflag_aph) w=10**vlogw c electron energy loss rate parameter -- - rqa = w + 13.6 * ev * rsa(te,ne,0.,k) + rqa = w + 13.6 * ev_aph * rsa(te,ne,0.,k) c----------------------------------------------------------------------c elseif (istabon .eq. 6) then c use spline fit to POST93 table data -c xuse=min(max(xdata(1),log(te/ev)),xdata(nxdata)) -c yuse=min(max(ydata(1),log10(ne)),ydata(nydata)) +c xuse=min(max(xdata_aph(1),log(te/ev_aph)),xdata_aph(nxdata_aph)) +c yuse=min(max(ydata_aph(1),log10(ne)),ydata_aph(nydata_aph)) c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -587,22 +587,22 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- xuse=rle yuse=rld - nxcoef=nxdata - nycoef=nydata + nxcoef=nxdata_aph + nycoef=nydata_aph tsval = gettime(sec4) - w = B2VAhL(xuse, yuse, 0, 0, xknots, yknots, nxcoef, nycoef, - . kxords, kyords, rqacoef, ldf, workh, iflag) + w = B2VAhL(xuse, yuse, 0, 0, xknots_aph, yknots_aph, nxcoef, nycoef, + . kxords_aph, kyords_aph, rqacoef, ldf_aph, workh, iflag_aph) totb2val = totb2val + gettime(sec4) - tsval c electron energy loss rate parameter -- - rqa = w + 13.6 * ev * rsa(te,ne,0.,k) + rqa = w + 13.6 * ev_aph * rsa(te,ne,0.,k) c----------------------------------------------------------------------c elseif (istabon .eq. 7) then c use polynomial fit from Bob Campbell - 8/93 -c Note that the 13.6 * ev * rsa(te,ne,k) is omitted here as Campbell +c Note that the 13.6 * ev_aph * rsa(te,ne,k) is omitted here as Campbell c has already added it in - rqa = svradp(te/ev,ne) + rqa = svradp(te/ev_aph,ne) c----------------------------------------------------------------------c elseif (istabon .gt. 7) then # write error message @@ -662,7 +662,7 @@ real function rra (te, ne, tau, k) ite = 0 ine = 0 c compute abscissae -- - rlte = log(te/ev) + rlte = log(te/ev_aph) rlne = max(htln(0),min(htln(htnn),log(ne))) c find ite -- 51 if (ite .lt. htnt-1) then @@ -706,7 +706,7 @@ real function rra (te, ne, tau, k) jr = 1 c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -730,11 +730,11 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- c----------------------------------------------------------------------c elseif (istabon .eq. 5) then c use spline fit to POST93 table data -c xuse=min(max(xdata(1),log(te/ev)),xdata(nxdata)) -c yuse=min(max(ydata(1),log10(ne)),ydata(nydata)) +c xuse=min(max(xdata_aph(1),log(te/ev_aph)),xdata_aph(nxdata_aph)) +c yuse=min(max(ydata_aph(1),log10(ne)),ydata_aph(nydata_aph)) c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -742,21 +742,21 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- xuse=rle yuse=rld - nxcoef=nxdata - nycoef=nydata + nxcoef=nxdata_aph + nycoef=nydata_aph - vlog10rra = B2VAhL(xuse, yuse, 0, 0, xknots, yknots, nxcoef, - . nycoef, kxords, kyords, rracoef, ldf, workh, iflag) + vlog10rra = B2VAhL(xuse, yuse, 0, 0, xknots_aph, yknots_aph, nxcoef, + . nycoef, kxords_aph, kyords_aph, rracoef, ldf_aph, workh, iflag_aph) rra=10**vlog10rra c----------------------------------------------------------------------c elseif (istabon .eq. 6) then c use spline fit to POST93 table data -c xuse=min(max(xdata(1),log(te/ev)),xdata(nxdata)) -c yuse=min(max(ydata(1),log10(ne)),ydata(nydata)) +c xuse=min(max(xdata_aph(1),log(te/ev_aph)),xdata_aph(nxdata_aph)) +c yuse=min(max(ydata_aph(1),log10(ne)),ydata_aph(nydata_aph)) c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -764,24 +764,24 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- xuse=rle yuse=rld - nxcoef=nxdata - nycoef=nydata + nxcoef=nxdata_aph + nycoef=nydata_aph tsval = gettime(sec4) - rra = B2VAhL(xuse, yuse, 0, 0, xknots, yknots, nxcoef, nycoef, - . kxords, kyords, rracoef, ldf, workh, iflag) + rra = B2VAhL(xuse, yuse, 0, 0, xknots_aph, yknots_aph, nxcoef, nycoef, + . kxords_aph, kyords_aph, rracoef, ldf_aph, workh, iflag_aph) totb2val = totb2val + gettime(sec4) - tsval c----------------------------------------------------------------------c elseif (istabon .eq. 7) then c use polynomial fit from Bob Campbell - 8/93 - rra = srecf(te/ev,ne) + rra = srecf(te/ev_aph,ne) c----------------------------------------------------------------------c elseif ((istabon>9 .and. istabon<14) .or. istabon .eq. 17) then c log-log interp on Stotler DEGAS2 tables jr = 1 c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -816,7 +816,7 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- endif c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -909,7 +909,7 @@ real function rsa (te, ne, tau, k) c----------------------------------------------------------------------c if (istabon .eq. 0) then # use analytic model (hydrogen) - a = te / (10*ev) + a = te / (10*ev_aph) rsa = 3.0e-14 * a*a / (3.0 + a*a) c----------------------------------------------------------------------c @@ -920,7 +920,7 @@ real function rsa (te, ne, tau, k) ite = 0 ine = 0 c compute abscissae -- - rlte = log(te/ev) + rlte = log(te/ev_aph) rlne = max(htln(0),min(htln(htnn),log(ne))) c find ite -- 51 if (ite .lt. htnt-1) then @@ -964,7 +964,7 @@ real function rsa (te, ne, tau, k) jr = 1 c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -988,11 +988,11 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- c----------------------------------------------------------------------c elseif (istabon .eq. 5) then c use spline fit to POST93 table data -c xuse=min(max(xdata(1),log(te/ev)),xdata(nxdata)) -c yuse=min(max(ydata(1),log10(ne)),ydata(nydata)) +c xuse=min(max(xdata_aph(1),log(te/ev_aph)),xdata_aph(nxdata_aph)) +c yuse=min(max(ydata_aph(1),log10(ne)),ydata_aph(nydata_aph)) c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -1000,21 +1000,21 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- xuse=rle yuse=rld - nxcoef=nxdata - nycoef=nydata + nxcoef=nxdata_aph + nycoef=nydata_aph - vlog10rsa = B2VAhL(xuse, yuse, 0, 0, xknots, yknots, nxcoef, - . nycoef, kxords, kyords, rsacoef, ldf, workh, iflag) + vlog10rsa = B2VAhL(xuse, yuse, 0, 0, xknots_aph, yknots_aph, nxcoef, + . nycoef, kxords_aph, kyords_aph, rsacoef, ldf_aph, workh, iflag_aph) rsa=10**vlog10rsa c----------------------------------------------------------------------c elseif (istabon .eq. 6) then c use spline fit to POST93 table data -c xuse=min(max(xdata(1),log(te/ev)),xdata(nxdata)) -c yuse=min(max(ydata(1),log10(ne)),ydata(nydata)) +c xuse=min(max(xdata_aph(1),log(te/ev_aph)),xdata_aph(nxdata_aph)) +c yuse=min(max(ydata_aph(1),log10(ne)),ydata_aph(nydata_aph)) c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -1022,17 +1022,17 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- xuse=rle yuse=rld - nxcoef=nxdata - nycoef=nydata + nxcoef=nxdata_aph + nycoef=nydata_aph tsval = gettime(sec4) - rsa = B2VAhL(xuse, yuse, 0, 0, xknots, yknots, nxcoef, nycoef, - . kxords, kyords, rsacoef, ldf, workh, iflag) + rsa = B2VAhL(xuse, yuse, 0, 0, xknots_aph, yknots_aph, nxcoef, nycoef, + . kxords_aph, kyords_aph, rsacoef, ldf_aph, workh, iflag_aph) totb2val = totb2val + gettime(sec4) - tsval c----------------------------------------------------------------------c elseif (istabon .eq. 7) then c use polynomial fit from Bob Campbell - 8/93 - rsa = sionf(te/ev,ne) + rsa = sionf(te/ev_aph,ne) c----------------------------------------------------------------------c elseif ((istabon>9 .and. istabon<14) .or. istabon .eq. 17) then @@ -1041,7 +1041,7 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- jr = 1 c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -1076,7 +1076,7 @@ c fractional parts of intervals (je,je+1) and (jd,jd+1) -- endif c compute abscissae -- - zloge=log(te/ev) + zloge=log(te/ev_aph) rle=max(rlemin, min(zloge,rlemax)) zlogd=log10(ne) rld=max(rldmin, min(zlogd,rldmax)) @@ -1328,7 +1328,7 @@ real function svdiss (te) b7 = 2.159670289222e-04 b8 = -4.928545325189e-06 - logt = log (te/ev) + logt = log (te/ev_aph) logsv = b0+logt*(b1+logt*(b2+logt*(b3+logt*(b4 . +logt*(b5+logt*(b6+logt*(b7+logt*b8))))))) svdiss = (1e-6)*exp(logsv) diff --git a/aph/aphread.m b/aph/aphread.m index 54a12c58..a9507070 100755 --- a/aph/aphread.m +++ b/aph/aphread.m @@ -753,8 +753,8 @@ c charge exchange rate (cm**3/sec): c (data from POST 93 tables) c Allocate arrays for spline fitting -- - nxdata=mpe # temperature - nydata=mpd # density + nxdata_aph=mpe # temperature + nydata_aph=mpd # density call gchange("Aphwrk",0) call splined1 @@ -774,61 +774,61 @@ call gchange("Aphwrk",0) integer i,j c Define data arrays -- - do i=1,nxdata - xdata(i)=ekpt(i) + do i=1,nxdata_aph + xdata_aph(i)=ekpt(i) enddo - do j=1,nydata - ydata(j)=dkpt(j) + do j=1,nydata_aph + ydata_aph(j)=dkpt(j) enddo - ldf=nxdata + ldf_aph=nxdata_aph c Define the order of the spline fit -c kxords=4 # cubic in x=log(temperature) -c kyords=4 # cubic in y=log10(density) +c kxords_aph=4 # cubic in x=log(temperature) +c kyords_aph=4 # cubic in y=log10(density) c Compute the coefficients -- c first, for ionization: - do i=1,nxdata - do j=1,nydata + do i=1,nxdata_aph + do j=1,nydata_aph if (istabon .eq. 5) then - fdata(i,j)=log10( wsveh(i,j,1) ) + fdata_aph(i,j)=log10( wsveh(i,j,1) ) elseif (istabon .eq. 6) then - fdata(i,j)=wsveh(i,j,1) + fdata_aph(i,j)=wsveh(i,j,1) endif enddo enddo - iflag = 1 - call s2copy (nxdata,nydata,fdata,1,nxdata,rsacoef,1,nxdata) - call B2INhT (xdata, nxdata, ydata, nydata, kxords, kyords, - . xknots, yknots, rsacoef, ldf, workh, iflag) + iflag_aph = 1 + call s2copy (nxdata_aph,nydata_aph,fdata_aph,1,nxdata_aph,rsacoef,1,nxdata_aph) + call B2INhT (xdata_aph, nxdata_aph, ydata_aph, nydata_aph, kxords_aph, kyords_aph, + . xknots_aph, yknots_aph, rsacoef, ldf_aph, workh, iflag_aph) c next, for recombination: - do i=1,nxdata - do j=1,nydata + do i=1,nxdata_aph + do j=1,nydata_aph if (istabon .eq. 5) then - fdata(i,j)=log10( wsveh0(i,j,1) ) + fdata_aph(i,j)=log10( wsveh0(i,j,1) ) elseif (istabon .eq. 6) then - fdata(i,j)=wsveh0(i,j,1) + fdata_aph(i,j)=wsveh0(i,j,1) endif enddo enddo - iflag = 1 - call s2copy (nxdata,nydata,fdata,1,nxdata,rracoef,1,nxdata) - call B2INhT (xdata, nxdata, ydata, nydata, kxords, kyords, - . xknots, yknots, rracoef, ldf, workh, iflag) + iflag_aph = 1 + call s2copy (nxdata_aph,nydata_aph,fdata_aph,1,nxdata_aph,rracoef,1,nxdata_aph) + call B2INhT (xdata_aph, nxdata_aph, ydata_aph, nydata_aph, kxords_aph, kyords_aph, + . xknots_aph, yknots_aph, rracoef, ldf_aph, workh, iflag_aph) c next, for hydrogen line radiation: - do i=1,nxdata - do j=1,nydata + do i=1,nxdata_aph + do j=1,nydata_aph if (istabon .eq. 5) then - fdata(i,j)=log10( wlemiss(i,j) ) + fdata_aph(i,j)=log10( wlemiss(i,j) ) elseif (istabon .eq. 6) then - fdata(i,j)=wlemiss(i,j) + fdata_aph(i,j)=wlemiss(i,j) endif enddo enddo - iflag = 1 - call s2copy (nxdata,nydata,fdata,1,nxdata,rqacoef,1,nxdata) - call B2INhT (xdata, nxdata, ydata, nydata, kxords, kyords, - . xknots, yknots, rqacoef, ldf, workh, iflag) + iflag_aph = 1 + call s2copy (nxdata_aph,nydata_aph,fdata_aph,1,nxdata_aph,rqacoef,1,nxdata_aph) + call B2INhT (xdata_aph, nxdata_aph, ydata_aph, nydata_aph, kxords_aph, kyords_aph, + . xknots_aph, yknots_aph, rqacoef, ldf_aph, workh, iflag_aph) return end diff --git a/api/api.v b/api/api.v index ec3cac98..a4f44abb 100755 --- a/api/api.v +++ b/api/api.v @@ -40,7 +40,7 @@ getprad(nx,ny,ngsp,te,ne:real,ng:real,afrac,atau, ***** Impurity_transport: # variable specifying radial transport rate of impurities dnimp real /1./ [m**2/s] -methimp integer /33/ # specifies interp. for finite diff. for (y,x) +methimp integer /33/ +input # specifies interp. for finite diff. for (y,x) # 66 is log interp., 77 inverse interp., # otherwise linear interp. csexpn real /0./ # exponent for reducing impurity || vel from @@ -138,14 +138,14 @@ z2datm(nt,nr,nn) _real # average Z**2 ***** Imslwrk: # working arrays for 3-d spline interpolation -nxdata integer -nydata integer +nxdata_api integer +nydata_api integer nzdata integer -xdata(1:nxdata) _real -ydata(1:nydata) _real +xdata_api(1:nxdata_api) _real +ydata_api(1:nydata_api) _real zdata(1:nzdata) _real -fdata(1:nxdata,1:nydata,1:nzdata) _real -ldf integer # first dimension of 3-d data array +fdata_api(1:nxdata_api,1:nydata_api,1:nzdata) _real +ldf_api integer # first dimension of 3-d data array mdf integer # second dimension of 3-d data array iflagi integer # input/output flag for 3-d spline routines nwork2 integer # size of array work2 @@ -154,18 +154,18 @@ nwork3 integer # size of array work3 work3(nwork3) _real # work array for B3INT iworki(10) integer # work array for B3VAL icont integer /0/ # input flag for B3VAL -kxords integer /4/ # order of spline fit versus x -# kxords=4 (default) is cubic interpolation -kyords integer /4/ # order of spline fit versus y -# kyords=4 (default) is cubic interpolation +kxords_api integer /4/ # order of spline fit versus x +# kxords_api=4 (default) is cubic interpolation +kyords_api integer /4/ # order of spline fit versus y +# kyords_api=4 (default) is cubic interpolation kzords integer /4/ # order of spline fit versus z # kzords=4 (default) is cubic interpolation -xknots(1:nxdata+kxords) _real -yknots(1:nydata+kyords) _real +xknots_api(1:nxdata_api+kxords_api) _real +yknots_api(1:nydata_api+kyords_api) _real zknots(1:nzdata+kzords) _real -emcoef(1:nxdata,1:nydata,1:nzdata) _real # spline coeff's for emissivity -z1coef(1:nxdata,1:nydata,1:nzdata) _real # spline coeff's for average-Z -z2coef(1:nxdata,1:nydata,1:nzdata) _real # spline coeff's for average-Z**2 +emcoef(1:nxdata_api,1:nydata_api,1:nzdata) _real # spline coeff's for emissivity +z1coef(1:nxdata_api,1:nydata_api,1:nzdata) _real # spline coeff's for average-Z +z2coef(1:nxdata_api,1:nydata_api,1:nzdata) _real # spline coeff's for average-Z**2 ***** P93fcn: # setup and evaluation routines for spline representation of data on impurity @@ -228,7 +228,7 @@ z2avgbs(te:real,nratio:real,ntau:real) real function mntau(MXMISO*MXMISO) real # usol(KXA*MXNZCH*MXMISO) real # sbar(KXA*MXMISO1) real # - zi(MXMISO*MXNZCH) real # + zi_api(MXMISO*MXNZCH) real # ***** Cyield: # Variables used for DIVIMP physical sputtering models diff --git a/api/apip93.m b/api/apip93.m index 6e77d379..d4ae720c 100755 --- a/api/apip93.m +++ b/api/apip93.m @@ -91,7 +91,7 @@ subroutine readpost1 (nget) subroutine splinem implicit none Use(P93dat) # nt,nr,nn -Use(Imslwrk) # nxdata,nydata,nzdata +Use(Imslwrk) # nxdata_api,nydata_api,nzdata external gallot, splinem1 c Construct 3-dimensional B-spline representation for impurity @@ -99,12 +99,12 @@ subroutine readpost1 (nget) c n*tau (data from POST '93 tables) c Allocate arrays for spline fitting -- - nxdata=nt # temperature - nydata=nr # density ratio + nxdata_api=nt # temperature + nydata_api=nr # density ratio nzdata=nn # n*tau - nwork2 = kyords*kzords + 3*max(kxords,kyords,kzords) + kzords + 2 - nwork3 = nxdata*nydata*nzdata + 2*max( kxords*(nxdata+1), - & kyords*(nydata+1), kzords*(nzdata+1) ) + nwork2 = kyords_api*kzords + 3*max(kxords_api,kyords_api,kzords) + kzords + 2 + nwork3 = nxdata_api*nydata_api*nzdata + 2*max( kxords_api*(nxdata_api+1), + & kyords_api*(nydata_api+1), kzords*(nzdata+1) ) call gallot("Imslwrk",0) call splinem1 @@ -117,69 +117,69 @@ call gallot("Imslwrk",0) subroutine splinem1 implicit none Use(P93dat) # tdatm,rdatm,ndatm,emdatm,z1datm,z2datm -Use(Imslwrk) # nxdata,nydata,nzdata,xdata,ydata,zdata,fdata,ldf,mdf, - # kxords,kyords,kzords,xknots,yknots,zknots,emcoef, +Use(Imslwrk) # nxdata_api,nydata_api,nzdata,xdata_api,ydata_api,zdata,fdata_api,ldf_api,mdf, + # kxords_api,kyords_api,kzords,xknots_api,yknots_api,zknots,emcoef, # z1coef,z2coef integer i,j,k external B3INT c Define data arrays -- - do i=1,nxdata - xdata(i)=log10( tdatm(i,1,1) ) + do i=1,nxdata_api + xdata_api(i)=log10( tdatm(i,1,1) ) enddo - do j=1,nydata - ydata(j)=log10( rdatm(1,j,1) ) + do j=1,nydata_api + ydata_api(j)=log10( rdatm(1,j,1) ) enddo do k=1,nzdata zdata(k)=log10( ndatm(1,1,k) ) enddo c Define the order of the spline fit -c kxords=4 # cubic in x=log10(temperature) -c kyords=4 # cubic in y=log10(density ratio) +c kxords_api=4 # cubic in x=log10(temperature) +c kyords_api=4 # cubic in y=log10(density ratio) c kzords=4 # cubic in z=log10(n*tau) - ldf=nxdata - mdf=nydata + ldf_api=nxdata_api + mdf=nydata_api iflagi = 1 # let B3INT choose knots c Compute the coefficients -- c first, for emissivity: - do i=1,nxdata - do j=1,nydata + do i=1,nxdata_api + do j=1,nydata_api do k=1,nzdata - fdata(i,j,k)=log10( emdatm(i,j,k) ) - emcoef(i,j,k)=fdata(i,j,k) + fdata_api(i,j,k)=log10( emdatm(i,j,k) ) + emcoef(i,j,k)=fdata_api(i,j,k) enddo enddo enddo - call B3INT (xdata, nxdata, ydata, nydata, zdata, nzdata, - & kxords, kyords, kzords, xknots, yknots, zknots, - & emcoef, ldf, mdf, work3, iflagi) + call B3INT (xdata_api, nxdata_api, ydata_api, nydata_api, zdata, nzdata, + & kxords_api, kyords_api, kzords, xknots_api, yknots_api, zknots, + & emcoef, ldf_api, mdf, work3, iflagi) c next, for average Z: - do i=1,nxdata - do j=1,nydata + do i=1,nxdata_api + do j=1,nydata_api do k=1,nzdata - fdata(i,j,k)=z1datm(i,j,k) - z1coef(i,j,k)=fdata(i,j,k) + fdata_api(i,j,k)=z1datm(i,j,k) + z1coef(i,j,k)=fdata_api(i,j,k) enddo enddo enddo - call B3INT (xdata, nxdata, ydata, nydata, zdata, nzdata, - & kxords, kyords, kzords, xknots, yknots, zknots, - & z1coef, ldf, mdf, work3, iflagi) + call B3INT (xdata_api, nxdata_api, ydata_api, nydata_api, zdata, nzdata, + & kxords_api, kyords_api, kzords, xknots_api, yknots_api, zknots, + & z1coef, ldf_api, mdf, work3, iflagi) c next, for average Z**2: - do i=1,nxdata - do j=1,nydata + do i=1,nxdata_api + do j=1,nydata_api do k=1,nzdata - fdata(i,j,k)=z2datm(i,j,k) - z2coef(i,j,k)=fdata(i,j,k) + fdata_api(i,j,k)=z2datm(i,j,k) + z2coef(i,j,k)=fdata_api(i,j,k) enddo enddo enddo - call B3INT (xdata, nxdata, ydata, nydata, zdata, nzdata, - & kxords, kyords, kzords, xknots, yknots, zknots, - & z2coef, ldf, mdf, work3, iflagi) + call B3INT (xdata_api, nxdata_api, ydata_api, nydata_api, zdata, nzdata, + & kxords_api, kyords_api, kzords, xknots_api, yknots_api, zknots, + & z2coef, ldf_api, mdf, work3, iflagi) return end @@ -189,8 +189,8 @@ call B3INT (xdata, nxdata, ydata, nydata, zdata, nzdata, real function emissbs (vte,vnr,vnt) implicit none real vnt,vnr,vte -Use(Imslwrk) # nxdata,nydata,nzdata,xdata,ydata,zdata, - # kxords,kyords,kzords,xknots,yknots,zknots,emcoef +Use(Imslwrk) # nxdata_api,nydata_api,nzdata,xdata_api,ydata_api,zdata, + # kxords_api,kyords_api,kzords,xknots_api,yknots_api,zknots,emcoef integer nxcoef,nycoef,nzcoef real vlogw,w,xuse,yuse,zuse @@ -204,18 +204,18 @@ real function emissbs (vte,vnr,vnt) c vnr = ng/ne density ratio c vnt = ne*tau [sec/m**3]. - xuse=min(max(xdata(1),log10(vte)),xdata(nxdata)) - yuse=min(max(ydata(1),log10(vnr)),ydata(nydata)) + xuse=min(max(xdata_api(1),log10(vte)),xdata_api(nxdata_api)) + yuse=min(max(ydata_api(1),log10(vnr)),ydata_api(nydata_api)) zuse=min(max(zdata(1),log10(vnt)),zdata(nzdata)) - nxcoef=nxdata - nycoef=nydata + nxcoef=nxdata_api + nycoef=nydata_api nzcoef=nzdata icont = 0 - vlogw=B3VAL (xuse, yuse, zuse, 0, 0, 0, xknots, yknots, zknots, - & nxcoef, nycoef, nzcoef, kxords, kyords, kzords, - & emcoef, ldf, mdf, icont, iworki, work2, iflagi) + vlogw=B3VAL (xuse, yuse, zuse, 0, 0, 0, xknots_api, yknots_api, zknots, + & nxcoef, nycoef, nzcoef, kxords_api, kyords_api, kzords, + & emcoef, ldf_api, mdf, icont, iworki, work2, iflagi) w=10**vlogw emissbs=w @@ -227,8 +227,8 @@ real function emissbs (vte,vnr,vnt) real function z1avgbs (vte,vnr,vnt) implicit none real vnt,vnr,vte -Use(Imslwrk) # nxdata,nydata,nzdata,xdata,ydata,zdata, - # kxords,kyords,kzords,xknots,yknots,zknots,z1coef +Use(Imslwrk) # nxdata_api,nydata_api,nzdata,xdata_api,ydata_api,zdata, + # kxords_api,kyords_api,kzords,xknots_api,yknots_api,zknots,z1coef integer nxcoef,nycoef,nzcoef real vz1,xuse,yuse,zuse @@ -242,18 +242,18 @@ real function z1avgbs (vte,vnr,vnt) c vnr = ng/ne density ratio c vnt = ne*tau [sec/m**3]. - xuse=min(max(xdata(1),log10(vte)),xdata(nxdata)) - yuse=min(max(ydata(1),log10(vnr)),ydata(nydata)) + xuse=min(max(xdata_api(1),log10(vte)),xdata_api(nxdata_api)) + yuse=min(max(ydata_api(1),log10(vnr)),ydata_api(nydata_api)) zuse=min(max(zdata(1),log10(vnt)),zdata(nzdata)) - nxcoef=nxdata - nycoef=nydata + nxcoef=nxdata_api + nycoef=nydata_api nzcoef=nzdata icont = 0 - vz1=B3VAL (xuse, yuse, zuse, 0, 0, 0, xknots, yknots, zknots, - & nxcoef, nycoef, nzcoef, kxords, kyords, kzords, - & z1coef, ldf, mdf, icont, iworki, work2, iflagi) + vz1=B3VAL (xuse, yuse, zuse, 0, 0, 0, xknots_api, yknots_api, zknots, + & nxcoef, nycoef, nzcoef, kxords_api, kyords_api, kzords, + & z1coef, ldf_api, mdf, icont, iworki, work2, iflagi) z1avgbs=vz1 return @@ -264,8 +264,8 @@ real function z1avgbs (vte,vnr,vnt) real function z2avgbs (vte,vnr,vnt) implicit none real vnt,vnr,vte -Use(Imslwrk) # nxdata,nydata,nzdata,xdata,ydata,zdata, - # kxords,kyords,kzords,xknots,yknots,zknots,z2coef +Use(Imslwrk) # nxdata_api,nydata_api,nzdata,xdata_api,ydata_api,zdata, + # kxords_api,kyords_api,kzords,xknots_api,yknots_api,zknots,z2coef integer nxcoef,nycoef,nzcoef real vz2,xuse,yuse,zuse @@ -279,18 +279,18 @@ real function z2avgbs (vte,vnr,vnt) c vnr = ng/ne density ratio c vnt = ne*tau [sec/m**3]. - xuse=min(max(xdata(1),log10(vte)),xdata(nxdata)) - yuse=min(max(ydata(1),log10(vnr)),ydata(nydata)) + xuse=min(max(xdata_api(1),log10(vte)),xdata_api(nxdata_api)) + yuse=min(max(ydata_api(1),log10(vnr)),ydata_api(nydata_api)) zuse=min(max(zdata(1),log10(vnt)),zdata(nzdata)) - nxcoef=nxdata - nycoef=nydata + nxcoef=nxdata_api + nycoef=nydata_api nzcoef=nzdata icont = 0 - vz2=B3VAL (xuse, yuse, zuse, 0, 0, 0, xknots, yknots, zknots, - & nxcoef, nycoef, nzcoef, kxords, kyords, kzords, - & z2coef, ldf, mdf, icont, iworki, work2, iflagi) + vz2=B3VAL (xuse, yuse, zuse, 0, 0, 0, xknots_api, yknots_api, zknots, + & nxcoef, nycoef, nzcoef, kxords_api, kyords_api, kzords, + & z2coef, ldf_api, mdf, icont, iworki, work2, iflagi) z2avgbs=vz2 return diff --git a/api/fmombal.m b/api/fmombal.m index 2970b841..9fefab1d 100755 --- a/api/fmombal.m +++ b/api/fmombal.m @@ -2,14 +2,14 @@ c!include "../sptodp.h" c!include "api.h" ! SEE FMOMBAL SUBROUTINE FOR INPUT PARAMETER SPECIFICATIONS - subroutine coulfric(amu,denz2,dloglam,mntau,zi,capm,capn, + subroutine coulfric(amu,denz2,dloglam,mntau,zi_api,capm,capn, > ela,elab,tempa) implicit real (a-h,o-z), integer (i-n) Use(Reduced_ion_constants) real capm(KXA,miso,KXA,miso),denz2(miso), > capn(KXA,KXA,miso,miso),mntau(miso,miso), > amu(miso),tempa(miso), - > ela(*),elab(*),zi(miso,nzch) + > ela(*),elab(*),zi_api(miso,nzch) c c COMPUTES COULOMB COLLISION FRICTION TERMS c ELA AND ELAB, WHERE THE FRICTION FORCE IS @@ -206,24 +206,24 @@ integer natomic(misotope) amat(i)=zero enddo -c**** Setup zi, mass-density, Z*e*density arrays - call setden(amu,den,denmass,denz,denz2,zi) +c**** Setup zi_api, mass-density, Z*e*density arrays + call setden(amu,den,denmass,denz,denz2,zi_api) c**** Setup force arrays (gradp, gradt stuff) call setforce(den,denz,denmass,epar,gradp,gradt, > tempa,uneut,qneut,umass,fmomenta,nuion) c***** Compute Coulomb friction coefficients la, lab - call coulfric(amu,denz2,dloglam,mntau,zi,capm,capn, + call coulfric(amu,denz2,dloglam,mntau,zi_api,capm,capn, > ela,elab,tempa) c***** Compute response matrices for inversion of momentum equations c***** in charge-state space - call zrespond(den,denmass,zi,ela,nuion,nurec,uresp, + call zrespond(den,denmass,zi_api,ela,nuion,nurec,uresp, > usave,caplams,fmomenta,ldir) c***** Solve average-ion equations in mass-isotope space - call mrespond(elab,zi,den,denmass,uresp,sbar, + call mrespond(elab,zi_api,den,denmass,uresp,sbar, > KXA*(miso+1),amat,nurec,umass,ldir) c***** Compute individual ion flows (for all charge states) @@ -240,7 +240,7 @@ call mzrespond(elab,uresp,sbar,caplam,caplams,usol, c COMPUTE FRICTION FORCES c call getfrict(friction,fricc,caplam,denmass,ela,nuion, - . nurec,usol,zi) + . nurec,usol,zi_api) c ... End timing impurity calculations. if (istimingon .eq. 1 .and. misotope .gt. 1) @@ -251,26 +251,26 @@ call getfrict(friction,fricc,caplam,denmass,ela,nuion, c---- End of subroutine fmombal ---------------------------------------- c----------------------------------------------------------------------- subroutine getfrict(friction,fricc,caplam,denmass,ela, - > nuion,nurec,usol,zi) + > nuion,nurec,usol,zi_api) implicit real (a-h,o-z), integer (i-n) Use(Reduced_ion_constants) real friction(miso,nzch), usol(KXA,nzch,miso), . fricc(miso,nzch,5) real nuion(miso,0:nzch), nurec(miso,nzch) real denmass(miso,0:nzch), caplam(KXA,miso), - > ela(3,3,miso), zi(miso,nzch) + > ela(3,3,miso), zi_api(miso,nzch) m = 1 !ONLY FRICTION FORCE do misa = 1,miso nzmax = natom(misa) do nz = 1,nzmax friction(misa,nz) = - > zi(misa,nz) * ( ela(m,1,misa)*usol(1,nz,misa) + > zi_api(misa,nz) * ( ela(m,1,misa)*usol(1,nz,misa) > + ela(m,2,misa)*usol(2,nz,misa) > + ela(m,3,misa)*usol(3,nz,misa) + caplam(m,misa) ) - fricc(misa,nz,1) = zi(misa,nz)*ela(m,1,misa)*usol(1,nz,misa) - fricc(misa,nz,2) = zi(misa,nz)*ela(m,2,misa)*usol(2,nz,misa) - fricc(misa,nz,3) = zi(misa,nz)*ela(m,3,misa)*usol(3,nz,misa) - fricc(misa,nz,4) = zi(misa,nz)*caplam(m,misa) + fricc(misa,nz,1) = zi_api(misa,nz)*ela(m,1,misa)*usol(1,nz,misa) + fricc(misa,nz,2) = zi_api(misa,nz)*ela(m,2,misa)*usol(2,nz,misa) + fricc(misa,nz,3) = zi_api(misa,nz)*ela(m,3,misa)*usol(3,nz,misa) + fricc(misa,nz,4) = zi_api(misa,nz)*caplam(m,misa) friction(misa,nz) = friction(misa,nz) > - denmass(misa,nz)*usol(m,nz,misa)* > (nuion(misa,nz) + nurec(misa,nz))*al32(m) @@ -316,7 +316,7 @@ implicit real (a-h,o-z), integer (i-n) one = 1.0e0 c***** Calculate pi0 for greatest accuracy. pi0=acos(-one) -c***** Set zero for initializing arrays. +c***** Set zero for initializi_aping arrays. zero=0.e0 c***** Integer constants ilam1 = 1 @@ -352,14 +352,14 @@ integer natomic(misotope) end c---- End of subroutine initmombal ------------------------------------- c----------------------------------------------------------------------- - subroutine mrespond(elab,zi,den,denmass,uresp,sbar, + subroutine mrespond(elab,zi_api,den,denmass,uresp,sbar, > nmat,amat,nurec,umass,ldir) implicit real (a-h,o-z), integer (i-n) Use(Reduced_ion_constants) integer nmat real elab(KXA,miso,KXA*miso),uresp(KXA,nzch,miso,NBA), > sbar(KXA,miso+1), sbar0(KXA*(MXMISO1)), ubar0(KXA*(MXMISO1)), - > zi(miso,*),amat(nmat,nmat),denmass(miso,0:nzch) + > zi_api(miso,*),amat(nmat,nmat),denmass(miso,0:nzch) real den(miso,0:nzch), nurec(miso,nzch) real usum(MXNZCH), bmat(KXA*(MXMISO1)*KXA*(MXMISO1)) integer nrow(KXA*(MXMISO1)), z1 @@ -427,9 +427,9 @@ c AND U(1,1,alpha) MUST BE EXPRESSED IN TERMS OF AI, UBAR do 10 m = 1,KXA if( (misa.ne.mise) .or. (m.ne.mflow) )then sbar(m,misa) = - > sdot(nz,zi(misa,1),miso,uresp(m,z1,misa,iforc),KXA) + > sdot(nz,zi_api(misa,1),miso,uresp(m,z1,misa,iforc),KXA) amat(m+i,1+j) = - > -sdot(nz,zi(misa,1),miso,uresp(m,z1,misa,iacci),KXA) + > -sdot(nz,zi_api(misa,1),miso,uresp(m,z1,misa,iacci),KXA) endif do 10 mpmisb = 1,KXA*miso do nzb = 1,nz @@ -438,7 +438,7 @@ c AND U(1,1,alpha) MUST BE EXPRESSED IN TERMS OF AI, UBAR > + uresp(m,nzb,misa,ilam3)*elab(ilam3,misa,mpmisb) enddo if( (misa.ne.mise) .or. (m.ne.mflow) ) - > amat(m+i,mpmisb) = sdot(nz,zi(misa,z1),miso,usum,1) + > amat(m+i,mpmisb) = sdot(nz,zi_api(misa,z1),miso,usum,1) c c PERFORM SUMS OVER MISA (ISOTOPE MASS A) c @@ -673,7 +673,7 @@ real amu(miso), tempa(miso) end c---- End of subroutine neomn ------------------------------------------ c----------------------------------------------------------------------- - subroutine printit(amu,den,epar,zi,ela, + subroutine printit(amu,den,epar,zi_api,ela, > caplam,sbar,usol,fmom,denmass,nuion,nurec, > gradp,gradt,parcurrent,xirl,tempa,uneut,umass) implicit real (a-h,o-z), integer (i-n) @@ -683,7 +683,7 @@ real den(miso,0:nzch), caplam(KXA,miso) real denmass(miso,0:nzch) real nuion(miso,0:nzch), nurec(miso,nzch), > gradp(miso,*), gradt(miso,*), ela(3,3,miso), - > zi(miso,nzch), fmom(KXA,nzch,miso), amu(miso), + > zi_api(miso,nzch), fmom(KXA,nzch,miso), amu(miso), > tempa(miso), uneut(miso) c***** Print out average-ion quantities @@ -710,7 +710,7 @@ real nuion(miso,0:nzch), nurec(miso,nzch), nzmax = natom(misa) do 55 nz = 1,nzmax force = fmom(m,nz,misa) - friction = zi(misa,nz) * ( ela(m,1,misa)*usol(1,nz,misa) + friction = zi_api(misa,nz) * ( ela(m,1,misa)*usol(1,nz,misa) > + ela(m,2,misa)*usol(2,nz,misa) > + ela(m,3,misa)*usol(3,nz,misa) + caplam(m,misa) ) friction = friction - denmass(misa,nz)*usol(m,nz,misa)* @@ -724,7 +724,7 @@ real nuion(miso,0:nzch), nurec(miso,nzch), sumflow = sumflow + denmass(misa,nz)*usol(1,nz,misa) sumaflow= sumaflow+ denmass(misa,nz)*abs(usol(1,nz,misa)) write(*,550)nz,usol(m,nz,misa),force,friction, - > den(misa,nz),zi(misa,nz),nuion(misa,nz),nurec(misa,nz) + > den(misa,nz),zi_api(misa,nz),nuion(misa,nz),nurec(misa,nz) else if( m.eq.2 )then pres = den(misa,nz)*tempa(misa) write(*,550)nz,-2.5*usol(m,nz,misa)*pres, @@ -756,12 +756,12 @@ real nuion(miso,0:nzch), nurec(miso,nzch), end c---- End of subroutine printit ---------------------------------------- c----------------------------------------------------------------------- - subroutine setden(amu,den,denmass,denz,denz2,zi) + subroutine setden(amu,den,denmass,denz,denz2,zi_api) implicit real (a-h,o-z), integer (i-n) Use(Reduced_ion_constants) real amu(miso), denz2(miso) real den(miso,0:nzch), denmass(miso,0:nzch), denz(miso,*) - real zi(miso,*) + real zi_api(miso,*) c***** Sum over isotopes for electron density-isotope 1. do misa = 2,miso do nza = 1,natom(misa) @@ -776,8 +776,8 @@ real zi(miso,*) denz2(misa)=zero do nza = 1,natom(misa) denmass(misa,nza) = bmass*den(misa,nza) - zi(misa,nza) = den(misa,nza)*real(nza)**2 - denz2(misa)=denz2(misa) + zi(misa,nza) + zi_api(misa,nza) = den(misa,nza)*real(nza)**2 + denz2(misa)=denz2(misa) + zi_api(misa,nza) totmass = totmass + denmass(misa,nza) enddo enddo @@ -785,16 +785,16 @@ real zi(miso,*) do misa = 2,miso denmass(misa,0) = amu(misa)*promas*den(misa,0) enddo -c***** Provide pedestal for zi to avoid matrix inversion singularity - ziped = 1.E-4 +c***** Provide pedestal for zi_api to avoid matrix inversion singularity + zi_apiped = 1.E-4 do misa=1,miso zisum = zero do nza=1,natom(misa) - zi(misa,nza)=zi(misa,nza)/denz2(misa) + ziped - zisum = zisum + zi(misa,nza) + zi_api(misa,nza)=zi_api(misa,nza)/denz2(misa) + zi_apiped + zisum = zisum + zi_api(misa,nza) enddo - do nza = 1,natom(misa) !zi normed to 1 - zi(misa,nza) = zi(misa,nza)/zisum + do nza = 1,natom(misa) !zi_api normed to 1 + zi_api(misa,nza) = zi_api(misa,nza)/zisum enddo enddo return @@ -1045,14 +1045,14 @@ call scopy(k2,a(nrow(l)+k+1),1,a(nrow(l)+1),1) end c---- End of subroutine uinvm2 ----------------------------------------- c----------------------------------------------------------------------- - subroutine zrespond(den,denmass,zi,ela,nuion,nurec, + subroutine zrespond(den,denmass,zi_api,ela,nuion,nurec, > uresp,usave,caplams,fmom,ldir) implicit real (a-h,o-z), integer (i-n) Use(Reduced_ion_constants) parameter( NSIZE = 2*KXA*KXA*MXNZCH ) real zmat(KNX,MXMISO), source(KXA*MXNZCH,MXMISO) real amat(KNX), asource(KXA*MXNZCH*NBA) - real den(miso,0:nzch), zi(miso,nzch), ela(KXA,KXA,miso) + real den(miso,0:nzch), zi_api(miso,nzch), ela(KXA,KXA,miso) real nuion(miso,0:nzch), nurec(miso,nzch) real uresp(KXA,nzch,miso,NBA), usave(KXA,nzch,miso) real fmom(KXA,nzch,miso),denmass(miso,0:nzch) @@ -1126,7 +1126,7 @@ c SAME (SEE SUBROUTINE ZSOURCE FOR DETAILS) asource(i)=zero enddo - call zsource(asource,zi,den,fmom(1,1,misa),denmass,misa,nz) + call zsource(asource,zi_api,den,fmom(1,1,misa),denmass,misa,nz) nforc = nkz*(iforc-1) nacci = nkz*(iacci-1) if( ldir.le.1 )then @@ -1156,11 +1156,11 @@ call zsource(asource,zi,den,fmom(1,1,misa),denmass,misa,nz) c c LOWER DIAGONAL ELEMENTS (IONIZATION FROM LOWER CHARGE STATE) c - ozi = al32(m)/zi(misa,jz) + ozi_api = al32(m)/zi_api(misa,jz) if (jz.gt.1)then jp = jz-1 if( (m.eq.mp) )then - amat(index) = ozi*denmass(misa,jp)*nuion(misa,jp) + amat(index) = ozi_api*denmass(misa,jp)*nuion(misa,jp) if( ldir.le.1 ) > asource(is) = asource(is) - amat(index)*usave(mp,jp,misa) endif @@ -1172,7 +1172,7 @@ c LOWER DIAGONAL ELEMENTS (IONIZATION FROM LOWER CHARGE STATE) amat(index) = ela(m,mp,misa) if( (m.eq.mp) ) > amat(index) = amat(index) - denmass(misa,jz)* - > ( nuion(misa,jz) + nurec(misa,jz) )*ozi + > ( nuion(misa,jz) + nurec(misa,jz) )*ozi_api if( ldir.le.1 ) > asource(is) = asource(is) - amat(index)*usave(mp,jz,misa) c @@ -1182,7 +1182,7 @@ c UPPER DIAGONAL ELEMENTS (RECOMBINATION FROM HIGHER CHARGE STATE) index = index + KXA jp = jz+1 if( (m.eq.mp) )then - amat(index) = ozi*denmass(misa,jp)*nurec(misa,jp) + amat(index) = ozi_api*denmass(misa,jp)*nurec(misa,jp) if( ldir.le.1 ) > asource(is) = asource(is) - amat(index)*usave(mp,jp,misa) endif @@ -1227,20 +1227,20 @@ call scopy(nkz,xsol(noff),1,uresp(1,1,misa,ntype),1) end c---- End of subroutine zrespond --------------------------------------- c----------------------------------------------------------------------- - subroutine zsource(source,zi,den,fmom,denmass,misa,nz) + subroutine zsource(source,zi_api,den,fmom,denmass,misa,nz) implicit real (a-h,o-z), integer (i-n) Use(Reduced_ion_constants) - real zi(miso,*), den(miso,0:nz), source(KXA,nz,NBA), fmom(KXA,*) + real zi_api(miso,*), den(miso,0:nz), source(KXA,nz,NBA), fmom(KXA,*) real denmass(miso,0:nz) do 20 jz= 1,nz do 20 m = 1,KXA if( m.eq.1 )then source(m,jz,ilam1) = one - source(m,jz,iacci) = denmass(misa,jz) * anorm / zi(misa,jz) + source(m,jz,iacci) = denmass(misa,jz) * anorm / zi_api(misa,jz) else if( (m.eq.ilam2) .or. (m.eq.ilam3) )then source(m,jz,m) = one endif - source(m,jz,iforc) = fmom(m,jz)/zi(misa,jz) + source(m,jz,iforc) = fmom(m,jz)/zi_api(misa,jz) 20 continue return end diff --git a/bbb/ext_neutrals.m b/bbb/ext_neutrals.m index cb92f1a8..9d477114 100755 --- a/bbb/ext_neutrals.m +++ b/bbb/ext_neutrals.m @@ -173,10 +173,10 @@ c call gchange("MCN_sources",0) #wsor, sni, ... see c ... Write uedge plasma data if(isechocmdon) then - cmd="call writemcnfile("""//trim(bkufile)//""","""//trim(runid)//""")" + cmd="call writemcnfile("""//trim(bkufile)//""","""//trim(runid_ext)//""")" print *,trim(cmd) else - call writemcnfile(bkufile,runid) + call writemcnfile(bkufile,runid_ext) end if c ... Process plasma data into degas2 background file format diff --git a/bbb/oderhs.m b/bbb/oderhs.m index 2d1d38bb..2ed62f83 100755 --- a/bbb/oderhs.m +++ b/bbb/oderhs.m @@ -8554,8 +8554,6 @@ cc real(Size4) ranf c ... Pause from BASIS if a ctrl_c is typed call ruthere -c ... Count Jacobian evaluations, both for total and for this case - ijactot = ijactot + 1 #note: ijactot set 0 in exmain if icntnunk=0 ijac(ig) = ijac(ig) + 1 if (svrpkg.eq.'nksol') write(*,*) ' Updating Jacobian, npe = ', @@ -8709,6 +8707,9 @@ ccc call pandf1 (xc, yc, iv, neq, t, yl, wk) c sparse row format. call csrcsc (neq, 1, 1, rcsc, icsc, jcsc, jac, ja, ia) +c ... Count Jacobian evaluations, both for total and for this case + ijactot = ijactot + 1 #note: ijactot set 0 in exmain if icntnunk=0 + c ... Accumulate cpu time spent here. if (istimingon .eq. 1) ttjstor = ttjstor + gettime(sec4) - tsjstor return diff --git a/bbb/odesetup.m b/bbb/odesetup.m index ad2b6cb3..df759c64 100755 --- a/bbb/odesetup.m +++ b/bbb/odesetup.m @@ -58,8 +58,6 @@ cc Use(Rccoef) * -- local variables integer lda, lenk, ngspon, nispon, nuspon, ntgspon, ifld, isor, id character*60 runid - integer iprt_tfcx_warn - data iprt_tfcx_warn/1/ #Former Aux module variables integer igsp @@ -204,13 +202,6 @@ call xerrab("") endif enddo -c ... Check attempt to use deprecated variables fchemywi and fchemywo - if (fchemywi .ne. 1. .or. fchemywo .ne. 1.) then - call remark('**Input error: change fchemywi --> fchemygwi(igsp) - . and fchemywo --> fchemygwo(igsp)') - call xerrab("") - endif - c ... Check consistency of cngmom; should be zero if inertial neutrals if (isupgon(1).eq. 1 .and. cngmom(1).ne.0) then call remark('*** WARNING, likely Error: cngmom=1, isupgon=1') @@ -264,11 +255,6 @@ call remark('*** ERROR: use cgengpl, not cgpl for neut eng loss') call remark('*** Warning: yinc=2 recommended when isphion=1') endif -c ... Check if tfcx or tfcy set in input file - if ((tfcx>1.e-10 .or. tfcy>1.e-10) .and. iprt_tfcx_warn==1) then - call remark('*** WARNING: tfcx,y not active; see tgas instead') - iprt_tfcx_warn = 0 - endif c ... Check if isnfmiy=1 when geometry is snowflake > SF15 if (geometry=="snowflake45" .or. geometry=="snowflake75" .or. . geometry=="snowflake105" .or. geometry=="snowflake135" .or. @@ -6596,7 +6582,7 @@ call gchange("Interp",0) call comp_vertex_vals # gen plasma/neut values at rm,zm(,,4) endif endif - if (iprint .ge. 1) write(6,*) "Interpolants created; mype =", mype + write(6,*) "Interpolants created; mype =", mype endif 100 continue @@ -7392,7 +7378,7 @@ call gchange("Convdiffeqn",0) delt = courant*min(delx/vrzmax, 0.5*delx**2/drzmax) write(*,*) 'delt = ',delt dtdx = delt/delx - delto = tend/float(ntim) + delto = tendoned/float(ntim) ito = 1 timo(1) = 0. @@ -7428,7 +7414,7 @@ cc dens(1) = (4.*dens(2) - dens(3))/3. # Update boundary vals gampzt(ix,ito) = gampz(ix) enddo endif - if (tim > tend) exit + if (tim > tendoned) exit enddo return end diff --git a/flx/flx.v b/flx/flx.v index 1ac47a7a..a430a74f 100755 --- a/flx/flx.v +++ b/flx/flx.v @@ -16,25 +16,25 @@ nsearch integer /NSEARCH/ # maximum number search regions for contouring package ***** Flxin: -istchkon integer /0/ +istchkon integer /0/ +gridgen # switch for imposing limits on the polar angle about (rmagx,zmagx) # of (r,z) points on flux contours. # =0 no limits # =1 limits are defined by theta_split, thetax, dtheta_exclude, # dtheta_overlap_sol and dtheta_overlap_pf -isthmmxn integer /1/ +isthmmxn integer /1/ +gridgen # switch that controls expressions for thetamin and thetamax when # istchkon=1 for 'dnbot' configurations: # =0 uses same theta_split for in/outboard halves # =1 uses pi/twopi in place of theta_split for in/outboard halves -dtheta_exclude(1:2) real [radians] /2*1.5/ +dtheta_exclude(1:2) real [radians] /2*1.5/ +gridgen # angular width of region where SOL flux contours are excluded # index 1 refers to inboard flux contours; 2 refers to outboard contours -dtheta_overlap_sol(1:2) real [radians] /2*0.5/ +dtheta_overlap_sol(1:2) real [radians] /2*0.5/ +gridgen # angular width over which SOL flux contours can overlap # with flux contours in the adjacent region. # index 1 refers to inboard flux contours; 2 refers to outboard contours -dtheta_overlap_pf(1:2) real [radians] /2*0.25/ +dtheta_overlap_pf(1:2) real [radians] /2*0.25/ +gridgen # angular width over which p.f. flux contours can overlap # with flux contours in the adjacent region. thetax real [radians] @@ -45,13 +45,13 @@ thetamin(1:2) real [radians] thetamax(1:2) real [radians] # computed maximum poloidal angle of flux contour data # index 1 refers to inboard flux contours; 2 refers to outboard contours -theta1fac real /1.0/ +theta1fac real /1.0/ +gridgen # multiplicative factor for theta1 boundary of inner-half pf region -theta2fac real /0.0/ +theta2fac real /0.0/ +gridgen # multiplicative factor for theta2 boundary of outer-half pf region -ymax1fac real /1.0/ +ymax1fac real /1.0/ +gridgen # multiplies zmagx to set upper bound on search region 1 -ymax2fac real /3.0/ +ymax2fac real /3.0/ +gridgen # multiplies zseps to set upper bound on search region 2 imagx integer # horizontal (R) index of the refined-EFIT cell containing the magnetic axis @@ -65,99 +65,99 @@ icutoff1 integer # horizontal (R) index of the refined-EFIT cell containing xcutoff1 jcutoff1 integer # vertical (Z) index of the refined-EFIT cell containing ycutoff1 -slpyt real /1.0/ +slpyt real /1.0/ +gridgen # scaling factor for radial mesh near double-null separatrices - # - mesh size is scaled by slpyt near inner separatrix # - mesh size is scaled by 1/slpyt near outer separatrix -slp2fac real /1.0/ +slp2fac real /1.0/ +gridgen # scale factor for radial mesh in asymmetric double-null configurations # (see function rho3dn for details) # = 1. for uniform mesh between separatrices # < 1. for finer mesh at innermost separatrix # > 1. for coarser mesh at innermost separatrix -slp3fac real /1.0/ +slp3fac real /1.0/ +gridgen # scale factor for radial mesh in asymmetric double-null configurations # (see function rho3dn for details) # = 1. for uniform mesh between separatrices # < 1. for finer mesh at outermost separatrix # > 1. for coarser mesh at outermost separatrix -psifac real /1.0005/ +psifac real /1.0005/ +gridgen # multiplicative factor in normalized separatrix flux # to ensure that contour is just outside the x-point. psi0sep1 real # normalized flux value at lower x-point psi0sep2 real # normalized flux value at upper x-point -psi0max_inner real /1.07/ +regrid +psi0max_inner real /1.07/ +regrid +gridgen # normalized flux value at wall on inboard side of magnetic axis -psi0max_outer real /1.07/ +regrid +psi0max_outer real /1.07/ +regrid +gridgen # normalized flux value at wall on outboard side of magnetic axis -psi0min2_upper real /0.98/ +regrid +psi0min2_upper real /0.98/ +regrid +gridgen # normalized flux value at private flux boundary near upper x-point -psi0min2_lower real /0.98/ +regrid +psi0min2_lower real /0.98/ +regrid +gridgen # normalized flux value at private flux boundary near lower x-point -psi0min1 real /0.98/ +regrid +psi0min1 real /0.98/ +regrid +gridgen # normalized flux value at innermost core flux surface -psi0min2 real /0.98/ +regrid +psi0min2 real /0.98/ +regrid +gridgen # normalized flux value at innermost private flux surface -psi0sep real /1.00001/ +regrid +psi0sep real /1.00001/ +regrid +gridgen # normalized flux value at separatrix flux surface (just slightly outside) -psi0max real /1.07/ +regrid +psi0max real /1.07/ +regrid +gridgen # normalized flux value at wall on outboard side of magnetic axis -psi0lim real +psi0lim real +gridgen # normalized flux value where core radial mesh can be concentrated -sfaclim real /1.0/ +sfaclim real /1.0/ +gridgen # multiplicative factor for slope of psi0(iy) distribution at psi0lim # sfaclim = 1 --> nearly uniform # sfaclim < 1 --> more concentrated near psi0lim -alfcy_inner real /.0001/ +regrid +alfcy_inner real /.0001/ +regrid +gridgen # exponential factor for distribution of SOL flux contours (inboard) # alfcy << 1 --> uniform # alfcy >> 1 --> concentrated near separatrix -alfcy real /.0001/ +regrid +alfcy real /.0001/ +regrid +gridgen # exponential factor for distribution of SOL flux contours (outboard) # alfcy << 1 --> uniform # alfcy >> 1 --> concentrated near separatrix -xoverlap(2) real /5.0,4.0/ +xoverlap(2) real /5.0,4.0/ +gridgen +gridgen # overlap parameters for inboard and outboard flux contours # (in units of dxefit) rho(0:nym) _real # temporary array for normalized flux values on B2/UEDGE flux surfaces tflx(0:nym) _real # temporary array for indicies of B2/UEDGE flux surfaces -psitop(1:jdim) _real +psitop(1:jdim) _real +gridgen # array of un-normalized pol. mag. flux values across the midplane -psibot(1:jdim) _real +psibot(1:jdim) _real +gridgen # array of un-normalized pol. mag. flux values across both divertor plates -iseqdskr integer /0/ +iseqdskr integer /0/ +gridgen # switch (=1) for reflecting eqdsk data about the midplane; # used with geometry='uppersn' and for upper half-mesh of geometry='dnull' -kymesh integer /1/ +kymesh integer /1/ +gridgen # option flag for setting flux contour values that define radial mesh # =0 user sets arrays psitop (midplane) and psibot (div plates) # =1 psitop and psibot are analytic forms (see rho1); adjust via alfcy # =2 psitop and psibot are analytic forms (see rho5); uniform core -xcutoff1 real /0./ [m] +xcutoff1 real /0./ [m] +gridgen # inboard cutoff radius for flux contours -ycutoff1 real /0./ [m] +ycutoff1 real /0./ [m] +gridgen # lower vertical cutoff for flux contours -mdsefit integer /0/ +mdsefit integer /0/ +gridgen # Set this flag to 1 if the data is read in from Mdsplus # Prevents inflx from calling readefit. ***** Workdn: # some working arrays for full double-null configurations -psi0_mp_inner(0:nym) _real +psi0_mp_inner(0:nym) _real +gridgen # normalized flux values at radial mesh points of inner midplane -psi0_mp_outer(0:nym) _real +psi0_mp_outer(0:nym) _real +gridgen # normalized flux values at radial mesh points of outer midplane -psi0_dp_lower_inner(0:nym) _real +psi0_dp_lower_inner(0:nym) _real +gridgen # normalized flux values at radial mesh points of lower inner divertor plate -psi0_dp_lower_outer(0:nym) _real +psi0_dp_lower_outer(0:nym) _real +gridgen # normalized flux values at radial mesh points of lower outer divertor plate -psi0_dp_upper_inner(0:nym) _real +psi0_dp_upper_inner(0:nym) _real +gridgen # normalized flux values at radial mesh points of upper inner divertor plate -psi0_dp_upper_outer(0:nym) _real +psi0_dp_upper_outer(0:nym) _real +gridgen # normalized flux values at radial mesh points of upper outer divertor plate ***** Inpf0: @@ -186,14 +186,14 @@ yminf0(nsearch) _real [m] # vertical cutoff for contours at upper edge of efit domain ymaxf0(nsearch) _real [m] # vertical cutoff for contours at lower edge of efit domain -istcvon integer /0/ +istcvon integer /0/ +gridgen # obsolete option istcvon=1; use altsearch=2 instead. -altsearch integer /0/ +altsearch integer /0/ +gridgen # flag to change search path for private flux surfaces # altsearch=0 --> search vertically up toward x-point # altsearch=1 --> search vertically down from x-point # altsearch=2 --> search diagonally down and in from x-point -isetpath integer /0/ +isetpath integer /0/ +gridgen # set path for finding flux surfaces in midplane # isetpath=0 --> search radially outward from left boundary # isetpath=1 --> search radially inward from magnetic axis @@ -231,7 +231,7 @@ ncmax1 integer # flux contouring variables on the refined grid plflux(jdim) _real # temporary array of poloidal flux values in a search region -mrfac integer /4/ +mrfac integer /4/ +gridgen # mesh refinement factor relative to EFIT nx4 integer # number of grid surfaces in x-direction on refined grid @@ -246,7 +246,7 @@ f(nx4,ny4) _real ijumpf(jdim) _integer # core/p.f. discontinuity on flux contour j occurs between # data points ijumpf(j) and ijumpf(j)+1 -dsjumpf real /0.1/ [m] +dsjumpf real /0.1/ [m] +gridgen # criterion for "jump" discontinuity in core/p.f. contours ilast(jdim) _integer # comment needed diff --git a/grd/grd.v b/grd/grd.v index 68d29808..5b566a34 100755 --- a/grd/grd.v +++ b/grd/grd.v @@ -39,47 +39,47 @@ integer niwdim integer ***** Analgrd: # common block with parameters for cylindrical & rectangular (IDEAL) grid -radm real [m] /-1.e-4/ #minimum "radius" of cylinder or slab -radx real [m] /0.04/ #maximum "radius" of cylinder or slab -rad0 real [m] /0./ #location of "radial" separ'x for cylinder or slab -rscalcore real [ ] /1./ #scale fac to change radial dimen. of core region -za0 real [m] /0./ #position of left-hand axial boundary position -zax real [m] /1./ #position of right-hand axial boundary position -zaxpt real [m] /.75/ #position of x-point; uniform grid to here -tiltang real [deg] /0./ #inc. tilt angle of the divertor plate from 90 deg. -ixsnog integer /1/ #poloidal ix where nonorthogonal mesh begins -zxpt_reset real [m] /0./ #reset x-pt2 here if >0 for single-region mesh -alfyt real /-2./ #exponent coeff. for y (radial) grid nonuiformity -tnoty real /0./ #shift if r(t)~del_r*tanh radial mesh profile -sratiopf real /0./ #ratio of effec. alfyt in priv. flux region; calc - #internally to give same expans. rate if =0. -alfxt real /4.0/ #exponent coeff. for x (axial) grid nonuiformity -isadjalfxt integer /0/ #=1 changes alfxt slightly for smooth dx at ixpt2 -tctr real /0./ #relative location of dx maximum;0, left; 1., right -bpolfix real [T] /.3/ #poloidal B-field for cartesian grid -btfix real [T] /5./ #total B-field for cartesian grid -rmajfix real [m] /1.0/ #reference major radius for cartesian/cylindrical grid -sigma_bpol real [] /0.0/ #exponent for poloidal field for cartesian/cylindrical grid -sigma_btor real [] /0.0/ #exponent for toroidal field for cartesian/cylindrical grid -isgdistort integer /0/ #switch to distort poloidal grid -agsindx real /0./ #max amplitude shift of poloidal cell face linear +radm real [m] /-1.e-4/ +gridgen #minimum "radius" of cylinder or slab +radx real [m] /0.04/ +gridgen #maximum "radius" of cylinder or slab +rad0 real [m] /0./ +gridgen #location of "radial" separ'x for cylinder or slab +rscalcore real [ ] /1./ +gridgen #scale fac to change radial dimen. of core region +za0 real [m] /0./ +gridgen #position of left-hand axial boundary position +zax real [m] /1./ +gridgen #position of right-hand axial boundary position +zaxpt real [m] /.75/ +gridgen #position of x-point; uniform grid to here +tiltang real [deg] /0./ +gridgen #inc. tilt angle of the divertor plate from 90 deg. +ixsnog integer /1/ +gridgen #poloidal ix where nonorthogonal mesh begins +zxpt_reset real [m] /0./ +gridgen #reset x-pt2 here if >0 for single-region mesh +alfyt real /-2./ +gridgen #exponent coeff. for y (radial) grid nonuiformity +tnoty real /0./ +gridgen #shift if r(t)~del_r*tanh radial mesh profile +sratiopf real /0./ +gridgen #ratio of effec. alfyt in priv. flux region; calc + +gridgen #internally to give same expans. rate if =0. +alfxt real /4.0/ +gridgen #exponent coeff. for x (axial) grid nonuiformity +isadjalfxt integer /0/ +gridgen #=1 changes alfxt slightly for smooth dx at ixpt2 +tctr real /0./ +gridgen #relative location of dx maximum;0, left; 1., right +bpolfix real [T] /.3/ +gridgen #poloidal B-field for cartesian grid +btfix real [T] /5./ +gridgen #total B-field for cartesian grid +rmajfix real [m] /1.0/ +gridgen #reference major radius for cartesian/cylindrical grid +sigma_bpol real [] /0.0/ +gridgen #exponent for poloidal field for cartesian/cylindrical grid +sigma_btor real [] /0.0/ +gridgen #exponent for toroidal field for cartesian/cylindrical grid +isgdistort integer /0/ +gridgen #switch to distort poloidal grid +agsindx real /0./ +gridgen #max amplitude shift of poloidal cell face linear #linear in radial index iy -agsrsp real /0./ #max amplitude shift of poloidal cell face linear +agsrsp real /0./ +gridgen #max amplitude shift of poloidal cell face linear #in real space - use either asgindx or asgrsp/not b -iynod integer /0/ #iy about which asgindx distortion is centered -rnod real /0./ #value of y about which asgrsp distortion centered -ixdstar integer /1/ #ix where poloidal distortion begins +iynod integer /0/ +gridgen #iy about which asgindx distortion is centered +rnod real /0./ +gridgen #value of y about which asgrsp distortion centered +ixdstar integer /1/ +gridgen #ix where poloidal distortion begins ***** Torannulus: -acore real [m] /0.5/ #minor radius of core-edge bdry for mhdgeo=2 -rm0 real [m] /2./ #major radius of annulus for mhdgeo=2 -edgewid real [m] /0.1/ #width of simulated edge region for mhdgeo=2 -dthlim real [rad] /1e-4/ #polodial angular width of limiter cell -bpol0 real [T] /0.2/ #poloidal B-field for mhdgeo=2 (annulus) -btor0 real [T] /2./ #toroidal B-field for mhdgeo=2 on axis -radf(0:nym+1,0:4) _real [m] # minor radius of cells for mhdgeo=2 -thpf(1:nxm,0:4) _real [m] # poloidal angle of cells for mhdgeo=2 -ibpmodel integer /0/ #=0, Bpol=bpol0; =1, Bpol=bpol0*rm0/R +acore real [m] /0.5/ +gridgen #minor radius of core-edge bdry for mhdgeo=2 +rm0 real [m] /2./ +gridgen #major radius of annulus for mhdgeo=2 +edgewid real [m] /0.1/ +gridgen #width of simulated edge region for mhdgeo=2 +dthlim real [rad] /1e-4/ +gridgen #polodial angular width of limiter cell +bpol0 real [T] /0.2/ +gridgen #poloidal B-field for mhdgeo=2 (annulus) +btor0 real [T] /2./ +gridgen #toroidal B-field for mhdgeo=2 on axis +radf(0:nym+1,0:4) _real [m] +gridgen # minor radius of cells for mhdgeo=2 +thpf(1:nxm,0:4) _real [m] +gridgen # poloidal angle of cells for mhdgeo=2 +ibpmodel integer /0/ +gridgen #=0, Bpol=bpol0; =1, Bpol=bpol0*rm0/R ***** Magmirror: zu(1:nxm,1:nym,0:4) _real [m] #axial postions of vertices & center (0) @@ -195,23 +195,23 @@ iwsla(niwdim) _integer ***** Inmesh: # input data from user -isspnew integer /0/ +isspnew integer /0/ +gridgen # flag for source of strike-point data # isspnew = 0 --> strike point data from eqdsk files # = 1 --> user-specified values (rstrike(1:2),zstrike(1:2)) -rstrike(1:2) real [m] +rstrike(1:2) real [m] +gridgen # radial position of inboard:outboard strike point -zstrike(1:2) real [m] +zstrike(1:2) real [m] +gridgen # vertical position of inboard:outboard strike point -istpnew integer /0/ +istpnew integer /0/ +gridgen # flag for source of first-seed-point data (top of mesh) # istpnew = 0 --> 1st seed point at midplane for "dnbot" # = 1 --> user-specified values (rtpnew(1:2),ztpnew(1:2)) -rtpnew(1:2) real [m] +rtpnew(1:2) real [m] +gridgen # radial position of first inboard:outboard seed point -ztpnew(1:2) real [m] +ztpnew(1:2) real [m] +gridgen # vertical position of first inboard:outboard seed point -istptest(1:2) integer /2*0/ +istptest(1:2) integer /2*0/ +gridgen # flag for algorithm that defines top-of-mesh on inboard, (1), # and outboard, (2), separatrix flux contours: # = 0 --> test on R only @@ -219,10 +219,10 @@ istptest(1:2) integer /2*0/ # = 2 --> test on Z only ilmax(1:2) integer # index of last angle surface ( at divertor plate ) in each region -seedxp(idim,noregs) _real +seedxp(idim,noregs) _real +gridgen # normalized distance along core segment of separatrix # measured from first seed point (=0.) to x-point (=100.) -seedxpxl(idim,noregs) _real +seedxpxl(idim,noregs) _real +gridgen # normalized distance along divertor leg of separatrix measured from # x-point (=0.) to last seed point (=100.) seed(idim,noregs) _real @@ -245,7 +245,7 @@ y0g(noregs) _real # vertical position of the first seed point on the separatrix ylast(noregs) _real # vertical position of the last seed point on the separatrix -isztest(1:2) integer /2*0/ +isztest(1:2) integer /2*0/ +gridgen # flag for algorithm that defines end-of-mesh on inboard, (1), # and outboard, (2), separatrix flux contours: # = 0 --> test on R only @@ -253,7 +253,7 @@ isztest(1:2) integer /2*0/ # = 2 --> test on Z only epslon_lim real /1.e-3/ # ratio of limiter guard-cell x-width to adjacent cell -dalpha real /5./ +dalpha real /5./ +gridgen # fuzziness or overlap (in degrees) of angle limits associ. with alpha1 ***** Transit: @@ -295,7 +295,7 @@ iendc(noregs) _integer ***** UEgrid: # input/output data for defining the grid in the UEDGE code -ixtop integer +ixtop integer +gridgen # ix index of top "angle" surface opposite the x-point ***** Mmod: @@ -352,11 +352,11 @@ rtop2(ntop2) _real [m] ztop2(ntop2) _real [m] # vertical position of data points on the outboard top-of-mesh # reference surface -istream integer /0/ +istream integer /0/ +gridgen # option parameter for defining fixed upstream reference surface # istream=0 midplane+cut(ismmon=1) or top-of-mesh(ismmon=2) # istream=1 user-defined upstream surface arrays -isupstreamx integer /0/ +isupstreamx integer /0/ +gridgen # option parameter for defining fixed upstream reference surface # isupstreamx=0 midplane+cut(ismmon=1) or top-of-mesh(ismmon=2) # isupstreamx=1 orthogonal surface (SOL and PF) through x-point @@ -392,23 +392,23 @@ rdnstream2(ndnstream2) _real [m] zdnstream2(ndnstream2) _real [m] # vertical position of data points on the outboard downstream # reference surface -iplate integer /0/ +iplate integer /0/ +gridgen # option parameter for defining divertor plates # iplate=0 orthogonal plates # iplate=1 user-defined divertor plates -nplate1 integer +nplate1 integer +gridgen # number of data points on the inboard divertor plate surface -rplate1(nplate1) _real [m] +rplate1(nplate1) _real [m] +gridgen # radial position of data points on the inboard divertor plate surface -zplate1(nplate1) _real [m] +zplate1(nplate1) _real [m] +gridgen # vertical position of data points on the inboard divertor plate # surface -nplate2 integer +nplate2 integer +gridgen # number of data points on the outboard divertor plate surface -rplate2(nplate2) _real [m] +rplate2(nplate2) _real [m] +gridgen # radial position of data points on the outboard divertor plate # surface -zplate2(nplate2) _real [m] +zplate2(nplate2) _real [m] +gridgen # vertical position of data points on the outboard divertor plate # surface ntop integer @@ -455,17 +455,17 @@ zplate0(nptmp) _real [m] # temporary array for subroutine meshmod dsnorm(idim) _real # temporary array for subroutine meshmod2 -wtold real /0.5/ +wtold real /0.5/ +gridgen # weight factor for spatial smoothing of angle-like mesh surfaces # wtold=1.0 ==> no smoothing # wtold=0.5 ==> (1,2,1) relative weighting of (j-1,j,j+1) # wtold=0.0 ==> (1,0,1) relative weighting of (j-1,j,j+1) -nsmooth integer /2/ +nsmooth integer /2/ +gridgen # number of times to apply the smoothing algorithm to each # angle-like surface after non-orthogonal plate construction -fuzzm real [m] /1.0e-08/ +fuzzm real [m] /1.0e-08/ +gridgen # a measure of the "fuzziness" in meshpoint coordinates -delmax real [m] /1.0e-08/ +delmax real [m] /1.0e-08/ +gridgen # estimated maximum deviation of mesh points from exact flux # surfaces (used only in subroutine smooth) wtmesh1 real /0.5/ @@ -475,7 +475,7 @@ wtmesh1 real /0.5/ wtm1(idim) _real # temporary array of weight factors for merging limiter mesh # with original (unmodified) mesh -dmix0 real /0./ +dmix0 real /0./ +gridgen # normalized poloidal mixing length for combining mesh0 with mesh12 # =0. --> abrupt change from orthogonal mesh to mesh12 at upstream # position @@ -509,20 +509,20 @@ dsmesh3(idim) _real dsmeshff(idim) _real # temporary array for subroutine meshff # distance from top-of-mesh to modified meshpoints (flamefront) -cwtffu real /1./ +cwtffu real /1./ +gridgen # exponent for weight function variation between x-point and flamefront # wt(ix)=wtff*fac(ix)**cwtff where fac(ix) varies linearly from # zero at xpoint to unity at flamefront. -cwtffd real /1./ +cwtffd real /1./ +gridgen # exponent for weight function variation between flamefront and plate # wt(ix)=wtff*fac(ix)**cwtff where fac(ix) varies linearly from # zero at plate to unity at flamefront. -isxtform integer /1/ +isxtform integer /1/ +gridgen # flag for choosing various forms of x(ix) near flamefront # =1 --> specify slope factor at flamefront only (default) # =2 --> specify slope factor at flamefront and upstream # =3 --> specify slope factor at flamefront, upstream and downstream -iswtform integer /0/ +iswtform integer /0/ +gridgen # flag for choosing various forms of weight factors near flamefront # =0 --> wt(ix)=wtff for all ix (default) # =1 --> wt(ix)~wtff*abs(ix-ixlimit)**cwtff between flamefront and @@ -543,7 +543,7 @@ slpxffd1 real # slope factor at downstream limit (plate) of flamefront mesh # slpxuu < 1 makes mesh finer # slpxuu > 1 makes mesh coarser -nxdff1 integer /0/ +nxdff1 integer /0/ +gridgen # number of cells between flame front and divertor plate # on inner leg (region 1) wtff2 real @@ -568,26 +568,26 @@ nxdff2 integer /0/ ***** Refinex: # data used for mesh refinement near the x-point -isrefxptn integer /1/ +isrefxptn integer /1/ +gridgen # flag for choosing x-point mesh refinement algorithm # =0 old interpolation method # =1 new flux-surface-based method -nxmod integer /2/ +nxmod integer /2/ +gridgen # number of 'upstream' poloidal cells (per quadrant) in the # original mesh that we modify by calling subroutine refinex -alfxptl real /1./ # use as alfxpt for cells below(l) the x-pt; +alfxptl real /1./ +gridgen # use as alfxpt for cells below(l) the x-pt; # frac=(i/(nxxpt+nxmod))**alfxpt for extra x-pt grid spacing below x-pt -alfxpt2l real /1./ # use as alfxpt2 for cells below(l) the x-pt; +alfxpt2l real /1./ +gridgen # use as alfxpt2 for cells below(l) the x-pt; # frac2=(i/(nxxpt+nxmod-1))**alfxpt2 for mixing fixed lngth & # flux-surface length in adding extra x-pt cells below x-pt -alfxptu real /1./ # use as alfxpt for cells above(u) the x-pt; +alfxptu real /1./ +gridgen # use as alfxpt for cells above(u) the x-pt; # frac=(i/(nxxpt+nxmod))**alfxpt for extra x-pt grid spacing above x-pt -alfxpt2u real /1./ # use as alfxpt2 for cells above(u) the x-pt; +alfxpt2u real /1./ +gridgen # use as alfxpt2 for cells above(u) the x-pt; # frac2=(i/(nxxpt+nxmod-1))**alfxpt2 for mixing fixed lngth & # flux-surface length in adding extra x-pt cells above x-pt -alfxpt real /1./ # work var for alfxptl,u for below(l)/above(u) x-pt +alfxpt real /1./ +gridgen # work var for alfxptl,u for below(l)/above(u) x-pt # frac=(i/(nxxpt+nxmod))**alfxpt for setting extra x-pt grid spacing -alfxpt2 real /1./ # work var for alfxptsl,u for below(l)/above(u) x-pt +alfxpt2 real /1./ +gridgen # work var for alfxptsl,u for below(l)/above(u) x-pt # frac2=(i/(nxxpt+nxmod-1))**alfxpt2 for mixing fixed lngth & # flux-surface length in adding extra x-pt cells rsu(0:nym+2) _real [m] @@ -617,7 +617,7 @@ rmm(0:nym,0:nxm) _real [m] zmm(0:nym,0:nxm) _real [m] # working array that contains z-coord's of angle-like surfaces # in refined region near x-point -nsmoothx integer /8/ +nsmoothx integer /8/ +gridgen # number of times to apply the smoothing algorithm to each # angle-like surface after mesh refinement near the x-point @@ -629,21 +629,21 @@ xdat(ndat) _real [m] # data values for xfcn tdat(ndat) _real # data values for t -kxmesh integer /1/ +kxmesh integer /1/ +gridgen # switch parameter for choosing model that defines x-mesh : # kxmesh=0 use old model (manual def. of seed points) # kxmesh=1 use linear*rational form for x(t) # kxmesh=2 use linear*exponential form for x(t) # kxmesh=3 use spline form for x(t) # kxmesh=4 use exponential+spline form for x(t) -slpxt real /1.0/ +slpxt real /1.0/ +gridgen # slope enhancement factor for x(t) near the "top" of the core # plasma -alfx(2) real /2*0.1/ +alfx(2) real /2*0.1/ +gridgen # log( dx(n+1)/dx(n) ) for 'gas' cells near the divertor plates -dxgas(2) real [m] +dxgas(2) real [m] +gridgen # poloidal size of first 'gas' cell at inboard and outboard plates -nxgas(2) integer +nxgas(2) integer +gridgen # number of poloidal 'gas' cells at inboard and outboard plates dt1 real # normalized-index distance (t) from inboard plate to 'extra' data @@ -655,19 +655,19 @@ dt2 real # point dx2 real [m] # physical distance (x) from outboard plate to 'extra' data point -ileft integer /0/ +ileft integer /0/ +gridgen # type of end condition at the left endpoint # = 1 first derivative specified by dleft # = 2 second derivative specified by dleft -dleft real /0.0/ [m] +dleft real /0.0/ [m] +gridgen # derivative at left endpoint -iright integer /0/ +iright integer /0/ +gridgen # type of end condition at the right endpoint # = 1 first derivative specified by dright # = 2 second derivative specified by dright -dright real /0.0/ [m] +dright real /0.0/ [m] +gridgen # derivative at right endpoint -kord integer /4/ # order of spline (=4 for cubic interpolation) +kord integer /4/ +gridgen # order of spline (=4 for cubic interpolation) ndatp2 integer /NDATP2/ # NDATP2 = NDAT+2 kntopt integer /1/ # knot selection option flag tknt(ndatp2+4) _real # knot locations @@ -697,15 +697,15 @@ iysptrxu integer # radial index of cells just inside the separatrix ***** Expseed: # variables to control exponential poloidal mesh in sub exponseed -fraclplt(2) real /.6,.6/ #frac of divertor leg with near-plate spacing -alfxdiv(2) real /.18,.18/ #exponential factor for cell-size expansion -alfxcore(2) real /.4,.4/ #exponential factor for cell-size expansion -shift_seed_leg(2) real /0.,0./ #shift in seed away from Xpt toward plate -shift_seed_core(2) real /1.,1./ #shift in seed away from Xpt toward top -nxlplt(2) integer /12,12/ #num poloidal cells in leg-plate region -nxlxpt(2) integer /4,4/ #num pol cells in leg-xpt region; +fraclplt(2) real /.6,.6/ +gridgen #frac of divertor leg with near-plate spacing +alfxdiv(2) real /.18,.18/ +gridgen #exponential factor for cell-size expansion +alfxcore(2) real /.4,.4/ +gridgen #exponential factor for cell-size expansion +shift_seed_leg(2) real /0.,0./ +gridgen #shift in seed away from Xpt toward plate +shift_seed_core(2) real /1.,1./ +gridgen #shift in seed away from Xpt toward top +nxlplt(2) integer /12,12/ +gridgen #num poloidal cells in leg-plate region +nxlxpt(2) integer /4,4/ +gridgen #num pol cells in leg-xpt region; #note nxlplt+nxlxpt = nxleg (must check) -fcorenunif real /0.8/ #frac pol core mesh with expon mesh +fcorenunif real /0.8/ +gridgen #frac pol core mesh with expon mesh ***** Xfcn: # contains various x-mesh functional forms diff --git a/pyexamples/box2/box2_in.py b/pyexamples/box2/box2_in.py index 83bd75b8..8eaa6b4b 100644 --- a/pyexamples/box2/box2_in.py +++ b/pyexamples/box2/box2_in.py @@ -4,6 +4,8 @@ # This is a Python version of the box case from Andreas Holm ########################################################################### +from uedge import bbb, com, grd + #-Geometry bbb.mhdgeo=-1 #-set cartesian geometry @@ -96,7 +98,7 @@ #-Neutral gas propeties -bbb.tfcx=5.;bbb.tfcy=5. #Franck-Condon temperatures +#bbb.tfcx=5.;bbb.tfcy=5. #Franck-Condon temperatures bbb.cngfx=1.;bbb.cngfy=1. #turn-on grad(T_g) flux if =1 bbb.cngflox=1.;bbb.cngfloy=0. #turn-on drift with ions if =1 bbb.cngmom = 1. #ion-gas momentum transfer diff --git a/pyexamples/box2/plotcontour.py b/pyexamples/box2/plotcontour.py index fd273279..559998a5 100644 --- a/pyexamples/box2/plotcontour.py +++ b/pyexamples/box2/plotcontour.py @@ -1,10 +1,12 @@ import matplotlib +import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec +from uedge import bbb, com gs=gridspec.GridSpec(2, 2) plt.figure(10) plt.subplot(gs[0,0]) -CS = plt.contour(com.zm[:,:,0], com.rm[:,:,0], bbb.te/ev) +CS = plt.contour(com.zm[:,:,0], com.rm[:,:,0], bbb.te/bbb.ev) plt.clabel(CS, inline=1, fontsize=10) params = {'mathtext.default': 'regular' } plt.rcParams.update(params) @@ -14,7 +16,7 @@ plt.subplot(gs[0,1]) -CS = plt.contour(com.zm[:,:,0], com.rm[:,:,0], bbb.ti/ev) +CS = plt.contour(com.zm[:,:,0], com.rm[:,:,0], bbb.ti/bbb.ev) plt.clabel(CS, inline=1, fontsize=10) plt.title('T$\mathregular{_i}$ [ev]') plt.grid(True) diff --git a/pyexamples/box2/runcase.py b/pyexamples/box2/runcase.py index d0a39448..4d77e082 100644 --- a/pyexamples/box2/runcase.py +++ b/pyexamples/box2/runcase.py @@ -9,18 +9,21 @@ import matplotlib.gridspec as gridspec import numpy as np import os +from uedge.uedgeplots import * +from uedge.rundt import rundt #-import some utilities for using OS ###execfile(os.path.join(os.environ['HOME'], 'utils/python/osfun.py')) +def paws(): + programPause = input("Press the key to continue...") ##-how to do this better? #-in .bashrc: "export PYLIB=/home/umansky1/PyUEDGE/uedge/pylib" -execfile(os.environ['PYLIB']+"/plotmesh.py") -execfile(os.environ['PYLIB']+"/plotcontour.py") -execfile(os.environ['PYLIB']+"/plotvar.py") -execfile(os.environ['PYLIB']+"/paws.py") -execfile(os.environ['PYLIB']+"/osfun.py") +#execfile(os.environ['PYLIB']+"/plotmesh.py") +#execfile(os.environ['PYLIB']+"/plotcontour.py") +#execfile(os.environ['PYLIB']+"/plotvar.py") +#execfile(os.environ['PYLIB']+"/osfun.py") #execfile("../../plotmesh.py") #execfile("../../pylib/plotvar.py") #execfile("../../pylib/plotr.py") @@ -32,7 +35,7 @@ #-read UEDGE settings -execfile("box2_in.py") +from box2_in import * #-do a quick preliminary run to set all internals @@ -41,7 +44,7 @@ #-show grid plotmesh() -wait = raw_input("PAUSING, PRESS ENTER TO CONTINUE...") +wait = input("PAUSING, PRESS ENTER TO CONTINUE...") #-this should be done in uefacets @@ -66,13 +69,13 @@ bbb.isbcwdt=1 bbb.dtreal = 1e-14; bbb.itermx=30; bbb.exmain() bbb.t_stop=1e0 - bbb.rundt() + rundt() bbb.dtreal=1e20; bbb.isbcwdt=0; bbb.exmain() hdf5_save('mycase.h5') ###execfile('plotcontour.py') - plotcontour() + from plotcontour import * paws() diff --git a/pyexamples/d3dHsm/runcase.py b/pyexamples/d3dHsm/runcase.py index f1356b92..e8fe478f 100644 --- a/pyexamples/d3dHsm/runcase.py +++ b/pyexamples/d3dHsm/runcase.py @@ -9,23 +9,21 @@ import matplotlib.gridspec as gridspec import numpy as np import os +from uedge.uedgeplots import * +from uedge.rundt import rundt + #-import some utilities for using OS ###execfile(os.path.join(os.environ['HOME'], 'utils/python/osfun.py')) #-in .bashrc: "export PYLIB=/home/umansky1/PyUEDGE/uedge/pylib" -execfile(os.environ['PYLIB']+"/plotmesh.py") -execfile(os.environ['PYLIB']+"/plotcontour.py") -execfile(os.environ['PYLIB']+"/plotvar.py") -execfile(os.environ['PYLIB']+"/paws.py") -execfile(os.environ['PYLIB']+"/osfun.py") plt.ion() #-read UEDGE settings -execfile("rd_d3dHsm_in.py") +from rd_d3dHsm_in import * #-do a quick preliminary run to set all internals @@ -33,8 +31,8 @@ #-show grid -plotmesh(iso=1) -wait = raw_input("PAUSING, PRESS ENTER TO CONTINUE...") +plotmesh() +wait = input("PAUSING, PRESS ENTER TO CONTINUE...") #-run to steady state @@ -42,12 +40,13 @@ bbb.isbcwdt=1 bbb.dtreal = 1e-14; bbb.itermx=30; bbb.exmain() bbb.t_stop=1e0 -bbb.rundt() +print(rundt) +rundt.rundt() bbb.dtreal=1e20; bbb.isbcwdt=0; bbb.exmain() #-show some results -plotvar(bbb.te/bbb.ev) +plotmeshval(bbb.te/bbb.ev) #-export the solution in hdf5 file hdf5_save('mycase.h5') diff --git a/pyscripts/double.py b/pyscripts/double.py deleted file mode 100755 index d9c2be74..00000000 --- a/pyscripts/double.py +++ /dev/null @@ -1,24 +0,0 @@ -from uedge import * - -def uedouble(): - com.nxleg=2*com.nxleg - com.nxcore=2*com.nxcore - com.nycore=2*com.nycore - com.nysol=2*com.nysol - if com.nxomit > 0: - if com.geometry=="dnbot": - com.nxomit = com.nxleg[0,0]+com.nxcore[0,0] + 2*com.nxxptxx + 1 - else: - com.nxomit=2*(com.nxomit-2*com.nxxptxx) + 2*com.nxxpt # assumes com.nxomit removes 1/2 SOL - if com.nyomitmx == 1: com.nysol = 1 - if grd.kxmesh == 4: - dxgasold=grd.dxgas - alfxold=grd.alfx - grd.alfx=alfxold/2. - grd.dxgas=dxgasold*(exp(grd.alfx)-1)/(exp(alfxold)-1) - grd.nxgas=2*grd.nxgas - bbb.restart=1 - bbb.newgeo=1 - bbb.gengrid=1 - bbb.isnintp=1 - grd.ixdstar = com.nxcore[0,1]+1 diff --git a/pyscripts/osfun.py b/pyscripts/osfun.py deleted file mode 100644 index 0264e1f0..00000000 --- a/pyscripts/osfun.py +++ /dev/null @@ -1,17 +0,0 @@ -def date(): - os.system("date") - -def ls(opts=""): - os.system("ls " + opts) - -def more(fname): - os.system("more " + fname) - -def pwd(): - os.system("pwd") - -def cp(opts=""): - os.system("cp " + opts) - -def mv(opts=""): - os.system("mv " + opts) diff --git a/pyscripts/paws.py b/pyscripts/paws.py deleted file mode 100644 index 72837492..00000000 --- a/pyscripts/paws.py +++ /dev/null @@ -1,2 +0,0 @@ -def paws(): - programPause = raw_input("Press the key to continue...") diff --git a/pyscripts/rdcontdt.py b/pyscripts/rdcontdt.py deleted file mode 100755 index 79773ad3..00000000 --- a/pyscripts/rdcontdt.py +++ /dev/null @@ -1,251 +0,0 @@ -# This file runs a time-dependent case using dtreal. First, obtain a converged -# solution for a (usually small) dtreal; xuedge must report iterm=1 at the end. -# Then adjust control parameters in rdinitdt; read this file, which reads rdinitdt. -# If a mistake is made, to restart this file without a Jacobian evaluation, -# be sure to reset iterm=1 (=> last step was successful) - -# IMPORT UEDGE (assuming starting from ipython before any imports) -from .uedge import * -from .ruthere import * -from .uexec import * -from numpy import zeros - -# IMPORT HDF5 routines for saving solutions below -from .hdf5 import * - -# INITIALIZE PARAMS -- SHOULD BE DONE IN MASTER SCRIPT OR TERMINAL SESSION -# BEFORE INVOKING THIS SCRIPT -uexec("uedge.rdinitdt",returns=globals()) -no = 0;yes = 1 -echo = no - -# Set precisions of floating point output -###import print_options -###print_options.set_float_precision(4) - -# Check if successful time-step exists (bbb.iterm=1) -if (bbb.iterm == 1): - print("Initial successful time-step exists") - bbb.dtreal = bbb.dtreal*bbb.mult_dt #compensates dtreal divided by mult_dt below -else: - print("*---------------------------------------------------------*") - print("Need to take initial step with Jacobian; trying to do here") - print("*---------------------------------------------------------*") - bbb.icntnunk = 0 - bbb.exmain() - ruthere() - bbb.dtreal = bbb.dtreal*bbb.mult_dt #compensates dtreal divided by mult_dt below - -if (bbb.iterm != 1): - print("*--------------------------------------------------------------*") - print("Error: converge an initial time-step first; then retry rdcontdt") - print("*--------------------------------------------------------------*") - exit() - -nx=com.nx;ny=com.ny;nisp=com.nisp;ngsp=com.ngsp;numvar=bbb.numvar -isteon=bbb.isteon -if (i_stor==0): - ni_stor = zeros((bbb.n_stor,nx+1+1,ny+1+1,nisp),"d") # set time storage arrays - up_stor = zeros((bbb.n_stor,nx+1+1,ny+1+1,nisp),"d") - te_stor = zeros((bbb.n_stor,nx+1+1,ny+1+1),"d") - ti_stor = zeros((bbb.n_stor,nx+1+1,ny+1+1),"d") - ng_stor = zeros((bbb.n_stor,nx+1+1,ny+1+1,ngsp),"d") - tg_stor = zeros((bbb.n_stor,nx+1+1,ny+1+1,ngsp),"d") - phi_stor = zeros((bbb.n_stor,nx+1+1,ny+1+1),"d") - tim_stor = zeros((bbb.n_stor),"d") - dtreal_stor = zeros((bbb.n_stor),"d") - nfe_stor = zeros((bbb.n_stor),"l") - dt_stor = (bbb.tstor_e - bbb.tstor_s)/(bbb.n_stor - 1) - -i_stor = max(i_stor,1) # set counter for storage arrays -bbb.dt_tot = max(bbb.dt_tot,0.) -nfe_tot = max(nfe_tot,0) -deldt_0 = bbb.deldt -isdtsf_sav = bbb.isdtsfscal - -if (bbb.ipt==1 and bbb.isteon==1): # set ipt to te(nx,iysptrx+1) if no user value - ipt = bbb.idxte[nx-1,com.iysptrx] #note: ipt is local, bbb.ipt global - -bbb.irev = -1 # forces second branch of irev in ii1 loop below -if (bbb.iterm == 1): # successful initial run with dtreal - bbb.dtreal = bbb.dtreal/bbb.mult_dt # gives same dtreal after irev loop -else: # unsuccessful initial run; reduce dtreal - bbb.dtreal = bbb.dtreal/(3*bbb.mult_dt) # causes dt=dt/mult_dt after irev loop - -if (bbb.initjac == 0): bbb.newgeo=0 -dtreal_sav = bbb.dtreal -bbb.itermx = bbb.itermxrdc -bbb.dtreal = bbb.dtreal/bbb.mult_dt #adjust for mult. to follow; mult_dt in rdinitdt -bbb.dtphi = bbb.rdtphidtr*bbb.dtreal -neq=bbb.neq -svrpkg=bbb.svrpkg.tostring().strip() -# -bbb.ylodt = bbb.yl -bbb.pandf1 (-1, -1, 0, bbb.neq, 1., bbb.yl, bbb.yldot) -fnrm_old = sqrt(sum((bbb.yldot[0:neq]*bbb.sfscal[0:neq])**2)) -if (bbb.initjac == 1): fnrm_old=1.e20 -print("initial fnrm =",fnrm_old) - -for ii1 in range( 1, bbb.ii1max+1): - if (bbb.ismfnkauto==1): bbb.mfnksol = 3 - # adjust the time-step - if (bbb.irev == 0): - # Only used after a dt reduc. success. completes loop ii2 for fixed dt - bbb.dtreal = min(3*bbb.dtreal,bbb.t_stop) #first move forward after reduction - bbb.dtphi = bbb.rdtphidtr*bbb.dtreal - if (bbb.ismfnkauto==1 and bbb.dtreal > bbb.dtmfnk3): bbb.mfnksol = -3 - bbb.deldt = 3*bbb.deldt - else: - # either increase or decrease dtreal; depends on mult_dt - bbb.dtreal = min(bbb.mult_dt*bbb.dtreal,bbb.t_stop) - bbb.dtphi = bbb.rdtphidtr*bbb.dtreal - if (bbb.ismfnkauto==1 and bbb.dtreal > bbb.dtmfnk3): bbb.mfnksol = -3 - bbb.deldt = bbb.mult_dt*bbb.deldt - - bbb.dtreal = min(bbb.dtreal,bbb.dt_max) - bbb.dtphi = bbb.rdtphidtr*bbb.dtreal - if (bbb.ismfnkauto==1 and bbb.dtreal > bbb.dtmfnk3): bbb.mfnksol = -3 - bbb.deldt = min(bbb.deldt,deldt_0) - bbb.deldt = max(bbb.deldt,bbb.deldt_min) - nsteps_nk=1 - print('--------------------------------------------------------------------') - print('--------------------------------------------------------------------') - print(' ') - print('*** Number time-step changes = ',ii1,' New time-step = ', bbb.dtreal) - print('--------------------------------------------------------------------') - - bbb.itermx = bbb.itermxrdc - if (ii1>1 or bbb.initjac==1): # first time calc Jac if initjac=1 - if (bbb.irev == 1): # decrease in bbb.dtreal - if (bbb.numrev < bbb.numrevjmax and \ - bbb.numrfcum < bbb.numrevjmax+bbb.numfwdjmax): #dont recom bbb.jac - bbb.icntnunk = 1 - bbb.numrfcum = bbb.numrfcum + 1 - else: # force bbb.jac calc, reset numrev - bbb.icntnunk = 0 - bbb.numrev = -1 # yields api.zero in next statement - bbb.numrfcum = 0 - bbb.numrev = bbb.numrev + 1 - bbb.numfwd = 0 - else: # increase in bbb.dtreal - if (bbb.numfwd < bbb.numfwdjmax and \ - bbb.numrfcum < bbb.numrevjmax+bbb.numfwdjmax): #dont recomp bbb.jac - bbb.icntnunk = 1 - bbb.numrfcum = bbb.numrfcum + 1 - else: - bbb.icntnunk = 0 #recompute jacobian for increase dt - bbb.numfwd = -1 - bbb.numrfcum = 0 - bbb.numfwd = bbb.numfwd + 1 - bbb.numrev = 0 #bbb.restart counter for dt reversals - bbb.isdtsfscal = isdtsf_sav - bbb.ftol = min(bbb.ftol_dt, 0.01*fnrm_old) - bbb.ftol = max(bbb.ftol, bbb.ftol_min) - exmain() # take a single step at the present bbb.dtreal - ruthere() - if (bbb.iterm == 1): - bbb.dt_tot = bbb.dt_tot + bbb.dtreal - nfe_tot = nfe_tot + bbb.nfe[0,0] - bbb.ylodt = bbb.yl - bbb.pandf1 (-1, -1, 0, bbb.neq, 1., bbb.yl, bbb.yldot) - fnrm_old = sqrt(sum((bbb.yldot[0:neq-1]*bbb.sfscal[0:neq-1])**2)) - if (bbb.dt_tot>=0.9999999*bbb.t_stop or fnrm_old= t_stop **') - print('*****************************************************') - break - - bbb.icntnunk = 1 - bbb.isdtsfscal = 0 - for ii2 in range( 1, bbb.ii2max+1): #take ii2max steps at the present time-step - if (bbb.iterm == 1): - bbb.itermx = bbb.itermxrdc - bbb.ftol = min(bbb.ftol_dt, 0.01*fnrm_old) - bbb.ftol = max(bbb.ftol, bbb.ftol_min) - bbb.exmain() - ruthere() - if (bbb.iterm == 1): - bbb.ylodt = bbb.yl - bbb.pandf1 (-1, -1, 0, bbb.neq, 1., bbb.yl, bbb.yldot) - fnrm_old = sqrt(sum((bbb.yldot[0:neq-1]*bbb.sfscal[0:neq-1])**2)) - print("Total time = ",bbb.dt_tot,"; Timestep = ",bbb.dtreal) - print("variable index ipt = ",ipt, " bbb.yl[ipt] = ",bbb.yl[ipt]) - dtreal_sav = bbb.dtreal - bbb.dt_tot = bbb.dt_tot + bbb.dtreal - nfe_tot = nfe_tot + bbb.nfe[0,0] - hdf5_save(savefn) - if (bbb.dt_tot>=0.999999999999*bbb.t_stop or fnrm_old= t_stop **') - print('*****************************************************') - break - print(" ") -## Store variables if a storage time has been crossed - if (bbb.dt_tot >= dt_stor*i_stor and i_stor<=bbb.n_stor): - i_stor1 = i_stor-1 - ni_stor[i_stor1,:,:,:] = ni - up_stor[i_stor1,:,:,:] = up - te_stor[i_stor1,:,:] = te - ti_stor1[i_stor1,:,:] = ti - ng_stor[i_stor1,:,:,:] = ng - phi_stor1[i_stor1,:,:] = phi - tim_stor[i_stor1] = bbb.dt_tot - nfe_stor[i_stor1] = nfe_tot - dtreal_stor[i_stor1] = bbb.dtreal - i_stor = i_stor + 1 - ## End of storage section - - if (bbb.dt_tot>=bbb.t_stop or fnrm_old ydmax): - itrouble=ii - print("** Fortran index of trouble making equation is:") - print(itrouble+1) - break - print("** Number of variables is:") - print("numvar = ", numvar) - print(" ") - iv_t = (itrouble).__mod__(numvar) + 1 - print("** Troublemaker equation is:") - print("iv_t = ",iv_t) - print(" ") - print("** Troublemaker cell (ix,iy) is:") - print(bbb.igyl[itrouble,]) - print(" ") - print("** Timestep for troublemaker equation:") - print(bbb.dtuse[itrouble]) - print(" ") - print("** yl for troublemaker equation:") - print(bbb.yl[itrouble]) - print(" ") - echo=oldecho - ######## end of idtroub script ############################## - - if (bbb.dtreal < bbb.dt_kill): - print(' ') - print('*************************************') - print('** FAILURE: time-step < dt_kill **') - print('*************************************') - break - bbb.irev = 1 - print('*** Converg. fails for bbb.dtreal; reduce time-step by 3, try again') - print('----------------------------------------------------------------- ') - bbb.dtreal = bbb.dtreal/(3*bbb.mult_dt) - bbb.dtphi = bbb.rdtphidtr*bbb.dtreal - if (bbb.ismfnkauto==1 and bbb.dtreal > bbb.dtmfnk3): bbb.mfnksol = -3 - bbb.deldt = bbb.deldt/(3*bbb.mult_dt) - bbb.iterm = 1 -echo = yes diff --git a/pyscripts/rdinitdt.py b/pyscripts/rdinitdt.py deleted file mode 100755 index db60867c..00000000 --- a/pyscripts/rdinitdt.py +++ /dev/null @@ -1,37 +0,0 @@ -# Setup file to run time-dependently using dtreal -# Change dtreal for starting dt and savefname to change pfb file name -# Once variables are set, read rdrundt to execute a time-dependent run -from uedge import * - -i_stor = 0 -nfe_tot = 0 -savefn = "savedt.hdf5" # name of hdf5 savefile written every timestep - -bbb.rdtphidtr = 1e20 # ratio dtphi/dtreal -bbb.ismfnkauto = 1 # if =1, mfnksol=3 for dtreal 0 -bbb.initjac = 0 # if=1, calc initial Jac upon reading rdcontdt -bbb.numrevjmax = 2 # number of dt reductions before Jac recalculated -bbb.numfwdjmax = 1 # number of dt increases before Jac recalculated -###bbb.ismmaxuc = 1 # =1 for intern calc mmaxu; =0,set mmaxu & dont chng -bbb.irev = -1 # flag to allow reduced dt advance after cutback -bbb.rlx = 0.9 # max. change in variable at each linear iteration -bbb.itermx = 7 # max. number of linear iterations allowed -bbb.tstor_s = 1e-5 # beginning time for storing solution -bbb.tstor_e = 1e-3 # ending time for storing solution -bbb.n_stor = 0 # number of linearly spaced storage points -bbb.ipt = 1 # index of variable; value printed at step - # if ipt not reset from unity, ipt=idxte(nx,iysptrx+1) - diff --git a/pyscripts/rundt.py b/pyscripts/rundt.py index 14e60652..4c01f9c9 100644 --- a/pyscripts/rundt.py +++ b/pyscripts/rundt.py @@ -16,35 +16,28 @@ class UeRun(): ''' Class containing information on run ''' - def __init__(self, n_stor = False): - from time import time + def __init__(self, *args, n_stor = False, **kwargs): from numpy import array from uedge import bbb, com # TODO: Add restore/recover from timeslice # TODO: Add plot timeslice directly # NOTE: No -> Utilize direct I/O from file instead - self.tstart = time() - self.numvar = bbb.numvar - self.nx = com.nx - self.ny = com.ny - self.ixpt1 = com.ixpt1[0] - self.ixpt2 = com.ixpt2[0] - self.iysptrx = com.iysptrx self.equationkey = array([b'te', b'ti', b'phi', b'up', b'ni', b'ng', b'tg']) - self.classvars = ['slice_ni', 'slice_ng', 'slice_up', 'slice_te', - 'slice_ti', 'slice_tg', 'slice_phi', 'slice_dttot', 'time', - 'fnorm', 'nfe', 'dt_tot', 'dtreal', 'ii1', 'ii2', 'ii1fail', - 'ii2fail', 'dtrealfail', 'itrouble', 'troubleeq', 'troubleindex', - 'ylfail', 'isteon', 'istion', 'isupon', 'isphion', 'isupgon', - 'isngon', 'istgon', 'ishymol', 'nisp', 'ngsp', 'nhsp', 'nhgsp', - 'nzsp', 'b0', 'ncore', 'pcoree', 'pcorei', 'internaleq', - 'internalspecies', 'yldotsfscalfail'] - - # Intiialize all variables to empty lists in class - for var in self.classvars: - self.__setattr__(var, []) + self.setupvars = {} + for var in [ 'numvar', 'isteon', 'istion', 'isupon', 'isphion', + 'isupgon', 'isngon', 'istgon', 'ishymol', 'nisp', 'ngsp', 'nhsp', + 'nhgsp', 'nzsp', 'b0', 'ncore', 'pcoree', 'pcorei', 'nx', 'ny', + 'iysptrx', 'ixpt1', 'ixpt2', 'neq' + ]: + try: + self.setupvars[var] = getattr(bbb, var) + except: + self.setupvars[var] = getattr(com, var) + self.setupvars['ixpt1'] = self.setupvars['ixpt1'][0] + self.setupvars['ixpt2'] = self.setupvars['ixpt2'][0] + super().__init__(*args, **kwargs) def itroub(self): ''' Function that displays information on the problematic equation ''' @@ -58,88 +51,79 @@ def itroub(self): 'Ion momentum', 'Ion density', 'Gas density', 'Gas temperature'] # Find the fortran index of the troublemaking equation - self.neq = bbb.neq - self.itrouble.append(deepcopy(argmax(abs(bbb.yldot*\ + self.classvars['itrouble'].append(deepcopy(argmax(abs(bbb.yldot*\ bbb.sfscal)[:bbb.neq])+1)) print("** Fortran index of trouble making equation is:\n{}".format(\ - self.itrouble[-1])) + self.classvars['itrouble'][-1])) # Print equation information print("** Number of equations solved per cell:\n numvar = {}\n"\ - .format(self.numvar)) + .format(bbb.numvar)) - self.troubleeq.append(mod(self.itrouble[-1]-1, bbb.numvar)+1) + self.classvars['troubleeq'].append(mod(self.classvars['itrouble'][-1]-1, bbb.numvar)+1) species = '' - self.internaleq.append([abs(x - self.itrouble[-1]).min() for x in \ + self.classvars['internaleq'].append([abs(x - self.classvars['itrouble'][-1]).min() for x in \ self.equations].index(0)) - if self.equations[self.internaleq[-1]].ndim == 3: - self.internalspecies.append( where(\ - self.equations[self.internaleq[-1]] == self.itrouble[-1])\ + if self.equations[self.classvars['internaleq'][-1]].ndim == 3: + self.classvars['internalspecies'].append( where(\ + self.equations[self.classvars['internaleq'][-1]] == self.classvars['itrouble'][-1])\ [-1][0] + 1) - species = ' of species {}'.format(self.internalspecies[-1]) + species = ' of species {}'.format(self.classvars['internalspecies'][-1]) else: - self.internalspecies.append(0) + self.classvars['internalspecies'].append(0) print('** Troublemaker equation is:\n{} equation{}: iv_t={}\n'\ - .format(equationsdescription[self.internaleq[-1]], species, - self.troubleeq[-1])) + .format(equationsdescription[self.classvars['internaleq'][-1]], species, + self.classvars['troubleeq'][-1])) # Display additional information about troublemaker cell - self.troubleindex.append(deepcopy(bbb.igyl[self.itrouble[-1]-1,])) - self.dtrealfail.append(deepcopy(bbb.dtreal)) - self.ylfail.append(deepcopy(bbb.yl[self.itrouble[-1]-1])) - self.yldotsfscalfail.append(deepcopy((bbb.yldot*bbb.sfscal)\ - [self.itrouble[-1]-1])) + self.classvars['troubleindex'].append(deepcopy(bbb.igyl[self.classvars['itrouble'][-1]-1,])) + self.classvars['dtrealfail'].append(deepcopy(bbb.dtreal)) + self.classvars['ylfail'].append(deepcopy(bbb.yl[self.classvars['itrouble'][-1]-1])) + self.classvars['yldotsfscalfail'].append(deepcopy((bbb.yldot*bbb.sfscal)\ + [self.classvars['itrouble'][-1]-1])) print('** Troublemaker cell (ix,iy) is:\n' + \ - '{}\n'.format(self.troubleindex[-1])) + '{}\n'.format(self.classvars['troubleindex'][-1])) print('** Timestep for troublemaker equation:\n' + \ - '{:.4e}\n'.format(self.dtrealfail[-1])) + '{:.4e}\n'.format(self.classvars['dtrealfail'][-1])) print('** yl for troublemaker equation:\n' + \ - '{:.4e}\n'.format(self.ylfail[-1])) - print('** yl*sfscal for troublemaker equation:\n' + \ - '{:.4e}\n'.format(self.yldotsfscalfail[-1])) + '{:.4e}\n'.format(self.classvars['ylfail'][-1])) + print('** yldot*sfscal for troublemaker equation:\n' + \ + '{:.4e}\n'.format(self.classvars['yldotsfscalfail'][-1])) - def savesuccess(self, ii1, ii2, savedir, savename, fnrm=None): + def savesuccess(self, savename, fnrm=None): from time import time from uedge import bbb from copy import deepcopy - self.time.append(time()) - if fnrm is None: - bbb.pandf1 (-1, -1, 0, bbb.neq, 1., bbb.yl, bbb.yldot) - self.fnorm.append(deepcopy((sum((bbb.yldot[:bbb.neq]*\ - bbb.sfscal[:bbb.neq])**2))**0.5)) - else: - self.fnorm.append(fnrm) - self.nfe.append(deepcopy(bbb.nfe)) - self.dt_tot.append(deepcopy(bbb.dt_tot)) - self.dtreal.append(deepcopy(bbb.dtreal)) - self.ii1.append(ii1) - self.ii2.append(ii2) - self.neq = bbb.neq try: - self.save('{}_UeCase.hdf5'.format(savefname.split('.')[0])) + self.classvars['time'].append(time()) except: pass - - self.save_intermediate(savedir, savename) + if 'fnorm' in self.classvars.keys(): + if fnrm is None: + bbb.pandf1 (-1, -1, 0, bbb.neq, 1., bbb.yl, bbb.yldot) + self.classvars['fnorm'].append(deepcopy((sum((bbb.yldot[:bbb.neq]*\ + bbb.sfscal[:bbb.neq])**2))**0.5)) + else: + self.classvars['fnorm'].append(fnrm) + for var in ['nfe', 'dt_tot', 'dtreal']: + if var in self.classvars: + self.classvars[var].append(deepcopy(getattr(bbb, var))) + self.save_intermediate(savename) def store_timeslice(self): from copy import deepcopy from uedge import bbb - self.slice_ni.append(deepcopy(bbb.ni)) - self.slice_ng.append(deepcopy(bbb.ng)) - self.slice_up.append(deepcopy(bbb.up)) - self.slice_te.append(deepcopy(bbb.te)) - self.slice_ti.append(deepcopy(bbb.ti)) - self.slice_tg.append(deepcopy(bbb.tg)) - self.slice_phi.append(deepcopy(bbb.phi)) - self.slice_dttot.append(deepcopy(bbb.dt_tot)) + for var in ['ni', 'ng', 'up', 'te', 'ti', 'tg', 'phi', 'dt_tot']: + self.classvars['slice_{}'.format(var)].append(deepcopy(getattr(bbb, var))) - def save_intermediate(self, savedir, savename): + + def save_intermediate(self, savename): from uedge.hdf5 import hdf5_save + from os.path import exists from uedge import bbb, com from h5py import File @@ -149,56 +133,39 @@ def save_intermediate(self, savedir, savename): for var in [ 'nisp', 'ngsp', 'nhsp', 'nhgsp', 'nzsp']: self.__setattr__(var, com.__getattribute__(var)) - try: - hdf5_save('{}/{}_last_ii2.hdf5'.format(savedir,savename)) - except: - print('Folder {} not found, saving output to cwd...'\ - .format(savedir)) - hdf5_save('{}_last_ii2.hdf5'.format(savename)) - # Try to store ready-made Case-file, if possible - try: - self.save('{}/{}_last_ii2_Case.hdf5'.format(savedir,savename)) - except: - pass + + if not exists('/'.join(savename.split('/')[:-1])): + print('Folder {} not found, saving output to cwd...'\ + .format(savename.split('/')[0])) + hdf5_save(savename) + + try: - file = File('{}/{}_last_ii2.hdf5'.format(savedir, savename), 'r+') + self.save(savename)#.replace('.hdf5', '_UeCase.hdf5')) except: - file = File('{}_last_ii2.hdf5'.format(savename), 'r+') - - file.require_group('convergence') - group = file['convergence'] - group.create_dataset('t_start', data=self.tstart) - group.create_dataset('numvar', data=self.numvar) - group.create_dataset('neq', data=self.neq) - group.create_dataset('nx', data=self.nx) - group.create_dataset('ny', data=self.ny) - group.create_dataset('ixpt1', data=self.ixpt1) - group.create_dataset('ixpt2', data=self.ixpt2) - group.create_dataset('iysptrx', data=self.iysptrx) - group.create_dataset('equationkey', data=self.equationkey) - group.create_dataset('itermx', data=self.itermx) - group.create_dataset('incpset', data=self.incpset) - group.create_dataset('ii1max', data=self.ii1max) - group.create_dataset('ii2max', data=self.ii1max) - group.create_dataset('numrevjmax', data=self.numrevjmax) - group.create_dataset('numfwdjmax', data=self.numfwdjmax) - group.create_dataset('numtotjmax', data=self.numtotjmax) - group.create_dataset('rdtphidtr', data=self.rdtphidtr) - group.create_dataset('deldt_min', data=self.deldt_min) - group.create_dataset('rlx', data=self.rlx) - - for var in self.classvars: - group.create_dataset(var, data=self.__getattribute__(var)) - - file.close() - - def convergenceanalysis(savefname, savedir='../solutions', fig=None, + hdf5_save(savename) + + with File(savename, 'r+') as file: + file.require_group('convergence') + group = file['convergence'] + group.create_dataset('t_start', data=self.tstart) + group.create_dataset('equationkey', data=self.equationkey) + for var, value in self.setupvars.items(): + group.create_dataset(var, data=value) + for var, value in self.classsetup.items(): + group.create_dataset(var, data=value) + for var in self.classvars: + group.create_dataset(var, data=self.classvars[var]) + print('Intermediate solution written to {}'.format(savename)) + + def convergenceanalysis(self, savefname, fig=None, xaxis = 'exmain', logx = False, color='k', label=None, ylim = (None, None)): from h5py import File from matplotlib.pyplot import subplots + from os.path import exists from datetime import timedelta from matplotlib.ticker import FuncFormatter from numpy import cumsum, ones @@ -210,49 +177,49 @@ def convergenceanalysis(savefname, savedir='../solutions', fig=None, print('Three subplots required for plots! Aborting...') return f = fig - try: - file = File('{}/{}'.format(savedir, savefname), 'r') - except: - print('File {}/{} not found. Aborting!'.format(savedir, - savefname)) + + if not exists(savefname): + print('File {} not found. Aborting!'.format(savefname)) return - data = file['convergence'] - try: + + with File(savefname, 'r') as file: data = file['convergence'] - except: - print('Convergence data not found in {}/{}. Aborting!'.format(\ - savedir, savefname)) - return - - if xaxis == 'exmain': - xlabel = 'Exmain calls' - xones = ones(data['ii2'][()].shape) - x = cumsum(xones) - elif xaxis == 'nfe': - xlabel = 'nfe internal iterations' - x = cumsum(data['nfe'][()][:, 0, 0]) - elif xaxis == 'time': - xlabel = 'Total wall-clock time [HH:MM]' - x = [timedelta(t - data['t_start'][()]) for t in data['time'][()]] - x = data['time'][()] - data['t_start'][()] + try: + data = file['convergence'] + except: + print('Convergence data not found in {}. Aborting!'.format(\ + savefname)) + return - if logx is True: - ax[0].loglog(x, data['fnorm'][()], '-', color=color, label=label) - ax[1].loglog(data['dt_tot'][()], data['fnorm'][()], '-', - color=color, label=label) - ax[2].loglog(x, data['dtreal'][()], '-', color=color, label=label) - else: - ax[0].semilogy(x, data['fnorm'][()], '-', color=color, label=label) - ax[1].semilogy(data['dt_tot'][()], data['fnorm'][()], '-', - color=color, label=label) - ax[2].semilogy(x, data['dtreal'][()], '-', color=color, - label=label) + if xaxis == 'exmain': + xlabel = 'Exmain calls' + xones = ones(data['ii2'][()].shape) + x = cumsum(xones) + elif xaxis == 'nfe': + xlabel = 'nfe internal iterations' + x = cumsum(data['nfe'][()][:, 0, 0]) + elif xaxis == 'time': + xlabel = 'Total wall-clock time [HH:MM]' + x = [timedelta(t - data['t_start'][()]) for t in data['time'][()]] + x = data['time'][()] - data['t_start'][()] + + if logx is True: + ax[0].loglog(x, data['fnorm'][()], '-', color=color, label=label) + ax[1].loglog(data['dt_tot'][()], data['fnorm'][()], '-', + color=color, label=label) + ax[2].loglog(x, data['dtreal'][()], '-', color=color, label=label) + else: + ax[0].semilogy(x, data['fnorm'][()], '-', color=color, label=label) + ax[1].semilogy(data['dt_tot'][()], data['fnorm'][()], '-', + color=color, label=label) + ax[2].semilogy(x, data['dtreal'][()], '-', color=color, + label=label) + ax[1].set_title('Total exmain evaluations: {}'.format\ + (len(data['dtreal'][()]))) ax[0].set_xlabel(xlabel) ax[1].set_xlabel('Accumulated plasma simualtion time [s]') ax[2].set_xlabel(xlabel) - ax[1].set_title('Total exmain evaluations: {}'.format\ - (len(data['dtreal'][()]))) ax[0].set_ylabel('Initial fnorm') ax[1].set_ylabel('Initial fnorm') ax[2].set_ylabel('Time-step (dtreal) [s]') @@ -269,106 +236,169 @@ def convergenceanalysis(savefname, savedir='../solutions', fig=None, return f - def failureanalysis(savefname, savedir='../solutions', equation=None): + def failureanalysis(self, savefname, equation=None, N=slice(None), geometry=False): from h5py import File + from os.path import exists from matplotlib.pyplot import subplots from numpy import histogram, zeros from matplotlib.collections import PolyCollection - f, ax = subplots(2,1, figsize=(10, 7)) - try: - file = File('{}/{}'.format(savedir, savefname), 'r') - except: - print('File {}/{} not found. Aborting!'.format(savedir, - savefname)) - try: - data = file['convergence'] - except: - print('Convergence data not found in {}/{}. Aborting!'.format(\ - savedir, savefname)) + # TODO: add option for plotting offending cell on geometric grid + if geometry is True: + f, ax = subplots(1,2, figsize=(12, 8), + gridspec_kw={'width_ratios': [1, 1.2]}) + f.subplots_adjust(top=0.98) + else: + f, ax = subplots(2,1, figsize=(10, 7)) + + if not exists(savefname): + print('File {} not found. Aborting!'.format(savefname)) return + + with File(savefname, 'r') as file: + try: + data = file['convergence'] + except: + print('Convergence data not found in {}. Aborting!'.format(\ + savefname)) + return - if equation is not None: - iequation = [x.decode('UTF-8') for x in data['equationkey']]\ - .index(equation) - # Bin the equation errors - nspecies = 1/(data['nisp'][()] + 1) - nbins = 7*data['nisp'][()] - counts, bins = histogram((data['internaleq'][()]+\ - data['internalspecies']*nspecies)-0.5, bins=nbins, - range=(-0.5,6.5)) - h, e = histogram(data['internaleq'][()] - 0.5, bins=7, - range=(-0.5,6.5)) - ax[0].bar([x for x in range(7)], h, width=1, edgecolor='k', - color=(0, 87/255, 183/255)) - ax[0].hist(bins[3*data['nisp'][()]:-1], bins[3*data['nisp'][()]:], - weights=counts[3*data['nisp'][()]:], edgecolor='k', - color=(255/255, 215/255, 0)) - ax[0].set_xticks(range(7)) - ax[0].set_xticklabels([x.decode('UTF-8') for x in \ - data['equationkey'][()]]) - ax[0].grid(linestyle=':', linewidth=0.5, axis='y') - ax[0].set_xlim((-0.5,6.5)) - ax[0].set_ylabel('Counts') - for i in range(7): - ax[0].axvline(i-0.5, linewidth=1, color='k') - - # Visualize error locations - nx = data['nx'][()] - ny = data['ny'][()] - ixpt1 = data['ixpt1'][()] - ixpt2 = data['ixpt2'][()] - iysptrx = data['iysptrx'][()] - frequency = zeros((nx+2, ny+2)) - - cells = [] - for i in range(nx+2): - for j in range(ny+2): - cells.append([[i-.5, j-.5], [i+.5, j-.5], - [i+.5, j+.5], [i-.5, j+.5]]) - - polys = PolyCollection(cells, edgecolors='k', linewidth=0.5, - linestyle=':') - - for i in range(len(data['itrouble'])): - coord = data['troubleindex'][()][i] - if equation is None: - frequency[coord[0], coord[1]] += 1 - elif iequation == data['internaleq'][()][i]: - frequency[coord[0], coord[1]] += 1 - - polys.set_cmap('binary') + if equation is not None: + iequation = [x.decode('UTF-8') for x in data['equationkey']]\ + .index(equation) + # Bin the equation errors + nspecies = 1/(data['nisp'][()] + 1) + nbins = 7*data['nisp'][()] + counts, bins = histogram((data['internaleq'][N]+\ + data['internalspecies'][N]*nspecies)-0.5, bins=nbins, + range=(-0.5,6.5)) + h, e = histogram(data['internaleq'][N] - 0.5, bins=7, + range=(-0.5,6.5)) + ax[0].bar([x for x in range(7)], h, width=1, edgecolor='k', + color=(0, 87/255, 183/255)) + ax[0].hist(bins[3*data['nisp'][()]:-1], bins[3*data['nisp'][()]:], + weights=counts[3*data['nisp'][()]:], edgecolor='k', + color=(255/255, 215/255, 0)) + ax[0].set_xticks(range(7)) + ax[0].set_xticklabels([x.decode('UTF-8') for x in \ + data['equationkey'][()]]) + ax[0].grid(linestyle=':', linewidth=0.5, axis='y') + ax[0].set_xlim((-0.5,6.5)) + ax[0].set_ylabel('Counts') + for i in range(7): + ax[0].axvline(i-0.5, linewidth=1, color='k') + + # Visualize error locations + nx = data['nx'][()] + ny = data['ny'][()] + ixpt1 = data['ixpt1'][()] + ixpt2 = data['ixpt2'][()] + iysptrx = data['iysptrx'][()] + frequency = zeros((nx+2, ny+2)) + + cells = [] + if geometry is True: + try: + rm = file['grid/com/rm'][()] + zm = file['grid/com/zm'][()] + except: + raise KeyError('Grid data not found in {}'.format(\ + savefname)) + for i in range(nx+2): + for j in range(ny+2): + cell = [] + for k in [1, 2, 4, 3]: + cell.append((rm[i, j, k], zm[i, j, k])) + cells.append(cell) + else: + for i in range(nx+2): + for j in range(ny+2): + cells.append([[i-.5, j-.5], [i+.5, j-.5], + [i+.5, j+.5], [i-.5, j+.5]]) + + polys = PolyCollection(cells, edgecolors='k', + linewidth=0.5-0.45*geometry, linestyle=':') + + for i in range(len(data['itrouble'][N])): + coord = data['troubleindex'][N][i] + if equation is None: + frequency[coord[0], coord[1]] += 1 + elif iequation == data['internaleq'][N][i]: + frequency[coord[0], coord[1]] += 1 + + if geometry is False: + polys.set_cmap('binary') + else: + polys.set_cmap('binary') polys.set_array(frequency.reshape(((nx+2)*(ny+2),))) cbar = f.colorbar(polys, ax=ax[1]) cbar.ax.set_ylabel('N trouble'+' for {}'.format(equation)*\ (equation is not None), va='bottom', labelpad=20) - ax[1].plot([.5, nx+.5, nx+.5, .5, .5], [.5, .5, ny+.5, ny+.5, .5], - 'k-', linewidth=1) ax[1].set_xlabel('Poloidal index') ax[1].set_ylabel('Radial index') ax[1].add_collection(polys) - ax[1].plot([.5, nx+.5],[iysptrx+.5, iysptrx+.5], 'k-', - linewidth=1) - ax[1].plot([ixpt1+.5, ixpt1+.5], [.5, iysptrx+.5], 'k-', - linewidth=1) - ax[1].plot([ixpt2+.5, ixpt2+.5], [.5, iysptrx+.5], 'k-', - linewidth=1) - - file.close() + if geometry is False: + ax[1].plot([.5, nx+.5, nx+.5, .5, .5], [.5, .5, ny+.5, ny+.5, .5], + 'k-', linewidth=1) + ax[1].plot([.5, nx+.5],[iysptrx+.5, iysptrx+.5], 'k-', + linewidth=1) + ax[1].plot([ixpt1+.5, ixpt1+.5], [.5, iysptrx+.5], 'k-', + linewidth=1) + ax[1].plot([ixpt2+.5, ixpt2+.5], [.5, iysptrx+.5], 'k-', + linewidth=1) + else: + ax[1].set_aspect('equal') + ax[1].set_ylim((0.95*zm.min(), 1.05*zm.max())) + ax[1].set_xlim((0.95*rm.min(), 1.05*rm.max())) + return f + def restorevalues(self): + ''' Restores the original UEDGE values ''' + from uedge import bbb + for key, value in self.orig.items(): + bbb.__setattr__(key, value) + + def exmain_isaborted(self): + ''' Checks if abort is requested ''' + from uedge import bbb + bbb.exmain() + # Abort flag set, abort case + if bbb.exmain_aborted == 1: + # Reset flag + bbb.exmain_aborted == 0 + # Restore parameters modified by script + try: + self.restorevalues(self) + except: + self.restorevalues() + return True + + def message(self, string, separator='-', pad='', seppad = '', + nseparator=1): + ''' Prints formatted self.message to stdout ''' + # TODO: add formatting for len>75 strings + msg = pad.strip() + ' ' + string.strip() + ' ' + pad.strip() + for i in range(nseparator): + print(seppad + separator*(len(msg)-2*len(seppad)) + seppad) + print(msg) + print(seppad + separator*(len(msg)-2*len(seppad)) + seppad) + + + def converge(self, dtreal=2e-9, ii1max=5000, ii2max=5, itermx=7, ftol=1e-5, dt_kill=1e-14, t_stop=100, dt_max=100, ftol_min = 1e-9, incpset=7, n_stor=0, storedist='lin', numrevjmax=2, numfwdjmax=1, numtotjmax=0, tstor=(1e-3, 4e-2), ismfnkauto=True, dtmfnk3=5e-4, mult_dt=3.4, - reset=True, initjac=False, rdtphidtr=1e20, deldt_min=0.04, rlx=0.9, + reset=True, initjac=True, rdtphidtr=1e20, deldt_min=0.04, rlx=0.9, - tsnapshot=None, savedir='../solutions', ii2increase=1.5): + tsnapshot=None, savedir='../solutions', ii2increase=0, savefname=None, + message=None): ''' Converges the case by increasing dt @@ -440,8 +470,11 @@ def converge(self, dtreal=2e-9, ii1max=5000, ii2max=5, itermx=7, ftol=1e-5, from numpy import linspace, logspace, log10, append from copy import deepcopy from uedge import bbb + from time import time from os.path import exists + self.tstart = time() + # TODO: add kwarg savename # Check if requested save-directory exists: if not, write to cwd if not exists(savedir): @@ -449,68 +482,52 @@ def converge(self, dtreal=2e-9, ii1max=5000, ii2max=5, itermx=7, ftol=1e-5, savedir)) savedir = '.' + if savefname is None: + autoname = True + try: + savefname = self.casename + except: + savefname = 'dtrun' + else: + autoname = False + + savefname = '{}/{}'.format(savedir, savefname) + + if autoname is True: + if exists('{}_last_ii2.hdf5'.format(savefname)): + i = 2 + while exists('{}_{}_last_ii2.hdf5'.format(savefname, i)): + i += 1 + savefname = '{}_{}'.format(savefname, i) + savefname = '{}_last_ii2.hdf5'.format(savefname) + + self.classvars = {} + for var in ['slice_ni', 'slice_ng', + 'slice_up', 'slice_te', 'slice_ti', 'slice_tg', 'slice_phi', + 'slice_dt_tot', 'time', 'fnorm', 'nfe', 'dt_tot', 'dtreal', 'ii1', + 'ii2', 'ii1fail', 'ii2fail', 'dtrealfail', 'itrouble', 'troubleeq', + 'troubleindex', 'ylfail', 'internaleq', 'internalspecies', + 'yldotsfscalfail']: + self.classvars[var] = [] self.orig = {} - self.orig['itermx'] = deepcopy(bbb.itermx) - self.orig['dtreal'] = deepcopy(bbb.dtreal) - self.orig['icntnunk'] = deepcopy(bbb.icntnunk) - self.orig['ftol'] = deepcopy(bbb.ftol) - self.orig['mfnksol'] = deepcopy(bbb.mfnksol) - self.orig['rlx'] = deepcopy(bbb.rlx) - self.orig['deldt'] = deepcopy(bbb.deldt) - self.orig['isdtsfscal'] = deepcopy(bbb.isdtsfscal) - self.orig['incpset'] = deepcopy(bbb.incpset) - + for var in ['itermx', 'dtreal', 'icntnunk', 'ftol', 'mfnksol', 'rlx', + 'deldt', 'isdtsfscal', 'incpset']: + self.orig[var] = deepcopy(getattr(bbb, var)) if numtotjmax == 0: numtotjmax = numrevjmax + numfwdjmax - - - self.itermx = itermx - self.incpset = incpset - self.ii1max = ii1max - self.ii2max = ii2max - self.numrevjmax = numrevjmax - self.numfwdjmax = numfwdjmax - self.numtotjmax = numtotjmax - self.rdtphidtr = rdtphidtr - self.deldt_min = deldt_min - self.rlx = rlx - + self.classsetup = {} + for var in ['itermx', 'incpset', 'ii1max', 'ii2max', 'numrevjmax', + 'numfwdjmax', 'numtotjmax', 'rdtphidtr', 'deldt_min', 'rlx']: + self.classsetup[var] = locals()[var] # TODO: Add variable to control reduciton factor? # TODO: Should dtreal = min(x, t_stop) actually be t_stop or dt_max? - def restorevalues(self): - ''' Restores the original UEDGE values ''' - for key, value in self.orig.items(): - bbb.__setattr__(key, value) - - def message(string, separator='-', pad='', seppad = '', - nseparator=1): - ''' Prints formatted message to stdout ''' - # TODO: add formatting for len>75 strings - message = pad.strip() + ' ' + string.strip() + ' ' + pad.strip() - for i in range(nseparator): - print(seppad + separator*(len(message)-2*len(seppad)) + seppad) - print(message) - print(seppad + separator*(len(message)-2*len(seppad)) + seppad) - def scale_timestep(scaling): ''' Increases/decreases time-step ''' bbb.dtreal *= scaling - def exmain_isaborted(self): - ''' Checks if abort is requested ''' - from uedge import bbb - bbb.exmain() - # Abort flag set, abort case - if bbb.exmain_aborted == 1: - # Reset flag - bbb.exmain_aborted == 0 - # Restore parameters modified by script - restorevalues(self) - return True - - def issuccess(self, t_stop, ftol_min): + def issuccess(self, t_stop, ftol_min, savefname): ''' Checks if case is converged ''' from datetime import timedelta from time import time @@ -520,25 +537,32 @@ def issuccess(self, t_stop, ftol_min): bbb.pandf1 (-1, -1, 0, bbb.neq, 1., bbb.yl, bbb.yldot) self.fnrm_old = sum((bbb.yldot[:bbb.neq-1]*\ bbb.sfscal[:bbb.neq-1])**2)**0.5 - self.savesuccess(ii1, ii2, savedir, bbb.label[0].strip(\ - ).decode('UTF-8'), self.fnrm_old) + self.savesuccess(savefname, + #bbb.label[0].strip(\).decode('UTF-8'), + self.fnrm_old) if (bbb.dt_tot>=t_stop or self.fnrm_old= t_stop'*(bbb.dt_tot >= t_stop), pad='**', separator='*') print('Total runtime: {}'.format(timedelta( seconds=round(time()-self.tstart)))) - restorevalues(self) + try: + self.restorevalues(self) + except: + self.restorevalues() return True def isfail(dt_kill): ''' Checks whether to abandon case ''' if (bbb.dtreal < dt_kill): - message('FAILURE: time-step < dt_kill', pad='**', + self.message('FAILURE: time-step < dt_kill', pad='**', separator='*') - restorevalues(self) + try: + self.restorevalues(self) + except: + self.restorevalues() return True def setmfnksol(ismfnkauto, dtmfnk3): @@ -572,22 +596,29 @@ def calc_fnrm(): bbb.rlx = rlx bbb.dtreal = dtreal bbb.ftol = ftol + ii1 = 0 + ii2 = 0 if (bbb.iterm == 1) and (bbb.ijactot > 0): - message('Initial successful time-step exists', separator='') + self.message('Initial successful time-step exists', separator='') else: - message('Need to take initial step with Jacobian; ' + \ + self.message('Need to take initial step with Jacobian; ' + \ 'trying to do here', seppad='*') # Ensure time-step is taken bbb.icntnunk = 0 + if message is not None: + print(message) # Take timestep and see if abort requested - if exmain_isaborted(self): + if self.exmain_isaborted(): return # Increase time # Verify time-step was successful if (bbb.iterm != 1): - restorevalues(self) - message('Error: converge an initial time-step first; then ' + \ - 're-execute command', seppad='*') + try: + self.restorevalues(self) + except: + self.restorevalues() + self.message('Error: converge an initial time-step first; ' + 'then re-execute command', seppad='*') return bbb.incpset = incpset bbb.itermx = itermx @@ -628,7 +659,7 @@ def calc_fnrm(): bbb.dtreal = min([bbb.dtreal, dt_max]) bbb.dtphi = rdtphidtr*bbb.dtreal bbb.deldt = min([bbb.deldt, deldt_0, deldt_min]) - message('Number of time-step changes = ''{} New time-step: {:.2E}\n'\ + self.message('Number of time-step changes = ''{} New time-step: {:.2E}\n'\ .format((ii1+1), bbb.dtreal), pad='***', nseparator=1) # Enter for every loop except first, unless intijac == True @@ -661,9 +692,14 @@ def calc_fnrm(): # Dynamically decrease ftol as the initial ftol decreases bbb.ftol = max(min(ftol, 0.01*self.fnrm_old),ftol_min) # Take timestep and see if abort requested - if exmain_isaborted(self): + if message is not None: + print(message) + if self.exmain_isaborted(): return - if issuccess(self, t_stop, ftol_min): + if bbb.iterm == 1: + self.classvars['ii1'].append(ii1) + self.classvars['ii2'].append(ii2) + if issuccess(self, t_stop, ftol_min, savefname): return bbb.icntnunk = 2 bbb.isdtsfscal = 0 @@ -676,14 +712,19 @@ def calc_fnrm(): bbb.ftol = max(min(ftol, 0.01*self.fnrm_old),ftol_min) # Take timestep and see if abort requested - message("Inner iteration #{}".format(ii2+1), nseparator=0, + self.message("Inner iteration #{}".format(ii2+1), nseparator=0, separator='') - if exmain_isaborted(self): + if message is not None: + print(message) + if self.exmain_isaborted(): return - if issuccess(self, t_stop, ftol_min): + if bbb.iterm == 1: + self.classvars['ii1'].append(ii1) + self.classvars['ii2'].append(ii2) + if issuccess(self, t_stop, ftol_min, savefname): return - message("Total time = {:.4E}; Timestep = {:.4E}".format(\ + self.message("Total time = {:.4E}; Timestep = {:.4E}".format(\ bbb.dt_tot-bbb.dtreal,bbb.dtreal), nseparator=0, separator='') # TODO: replace with more useful information @@ -706,19 +747,528 @@ def calc_fnrm(): self.itroub() ''' ISFAIL ''' if isfail(dt_kill): - self.save_intermediate(savedir, bbb.label[0].strip()\ - .decode('UTF-8')) + self.save_intermediate(savefname) +# bbb.label[0].strip().decode('UTF-8')) break irev = 1 - message('Converg. fails for bbb.dtreal; reduce time-step by '+\ + self.message('Converg. fails for bbb.dtreal; reduce time-step by '+\ '3, try again', pad = '***', nseparator=0) scale_timestep(1/(3*mult_dt)) bbb.dtphi = rdtphidtr*bbb.dtreal bbb.deldt *= 1/(3*mult_dt) setmfnksol(ismfnkauto, dtmfnk3) # bbb.iterm = 1 -# bbb.iterm = -1 # Ensure subsequent repetitions work as intended + # ii1max iterations completed without convergence + # Indicate solver failure regardless of whether last solve was completer + # or not. + bbb.iterm = -1 # Ensure subsequent repetitions work as intended + # Return False to signal fail + return False + + + def continuation_solve(self, + var, + target=None, + savedir=None, + index=None, + dt=0.02, + dtreal=1e-5, + ftol=1e-5, + iicont_max=7, + iicont_fail_max=3, + cutoff=1e6, + itermx=5, + incpset=8, + initres=1000, + deltafac=2, + dtdeltafac=2, + staticiter_max=3, + dtlim=1e4, + ii1max=150, + saveres=7, + **kwargs): + """ Solves var up to target using continuation method + var - string for variable, or nested dict of variables and targets + If nested dict, top level is variable name and subdict contains + keys 'target' and optionally 'index'. These values override the + kwargs of the solver. If index is not set, defaults to None + target - target value for variable + index - tuple containing index of var if var is array. + In case of extended array slicing, e.g. setting ranges of + var, index must be defined as tuples of indices and/or + slice objects. E.g. to set var[1:4,:,0], one would use + index = (slice(1,4), slice(None), 0). If ranges of var is + set, ensure target has the matching, appropriate dimensions + or is a float + dt - time-step used when solving using continuation + staticiter_max - maximum consequtive static iterations before callong + time-depenent run + dtreal - time-step used when resorting to time-dependent solve + ftol - ftol to be attained + iicont_max - iterations in continuation loop + iicont_fail_max - max fails in continuation loop before going time-dependent + cutoff - cutoff value for required interations at current delta to attain target + itermx - max number of iterations to take within time-step + incpset - max iterations before re-calculating jacobian + initres - inital resolution, i.e. delta = target-var/initres + deltafac - factor by which delta is increased after iicont_max increases + dtdeltafac - factor by which delta is increased by for time-dependent run + dtlim - Cutoff on (target-var)/delta for which a time-dependent run is triggered + ii1max - maximum number of dt iterations to take + saveres - how many successful iterations are allowed between saves are dumped + kwargs passed to rundt + """ + from uedge import bbb + from time import time + from os import makedirs + from copy import deepcopy + from shutil import rmtree + from os.path import exists + from numpy import array + from Forthon import package, packageobject + # TODO: icntnunk=0 on fail only? + self.tstart = time() + # ==== HELPERS ==== + def changevar(): + """ Changes var by delta """ + setvar(min(self.lastsuccess + self.delta, 1)) + + def setvar(value): + """ Sets var to value """ + for key, subdict in self.var.items(): + newvar = subdict['origvar'] + value * subdict['deltavar'] + if subdict['index'] is None: + setattr(subdict['pkgobj'], key, newvar) + else: + getattr(subdict['pkgobj'], key)[subdict['index']] = newvar + + def getvar(key, subdict): + """ Returns current value of var """ + # Variable is array: only set requested index + if 'index' not in subdict: + return getattr(subdict['pkgobj'], key) + else: + try: + return getattr(subdict['pkgobj'], key)[subdict['index']] + except: + return getattr(subdict['pkgobj'], key) + + def issuccess(): + """ Checks whether case is converged at target value of var """ + from time import time + from datetime import timedelta + self.lastsuccess += min(self.delta, 1-self.lastsuccess) + if self.lastsuccess >= 1: + print('\n===== TARGET VALUE ACHIEVED: REDUCE FNORM ====') + staticiter() + bbb.dtreal = 1e20 + + print('++++++++++++++++++++++++++++++++++++++++') + print('+++++ CONTINUATION SOLVE SUCCEEDED +++++') + print('++++++++++++++++++++++++++++++++++++++++') + print('Total runtime: {}'.format(timedelta( + seconds=round(time()-self.tstart)))) + self.savesuccess('SUCCESS_{}.hdf5'.format(self.savedir)) + self.restorevalues() + return True + else: + if self.isave >= self.saveres: + self.isave = 0 + self.savesuccess(self.savefname.format('{:.3f}'.format(self.lastsuccess).replace('.','p'))) + else: + self.isave += 1 + + def printinfo(printdelta=True): + """ Print current solve status """ + # TODO: add printing of indices? + ljust = 42 + print(' -Variable(s) being solved:') + for key, _ in self.var.items(): + print(' '.ljust(ljust-3), '-',key) + print('{}{:.3f}%'.format(' -Progress'.ljust(ljust), + self.lastsuccess*100 + )) + if printdelta is not False: + print('{}{:.3f}%'.format(' -Advancing by'.ljust(ljust), + self.delta*100 + )) + print('{}{}'.format(' -Steps to target at current delta:'.ljust(ljust), + (int((1 - self.lastsuccess)/self.delta)))) + print('') + + def isdeltatoosmall(): + """ Returns True if delta is below cutoff """ + if self.delta < 1/self.cutoff: + return True + return False + + def dostaticiter(): + """ Returns True if delta is above dt-run threshold """ + if self.delta > 1/self.dtlim: + return True + return False + + + def staticiter(): + """ Tries to reduce initial fnrm by iterating at constant var """ + from uedge import bbb + from copy import deepcopy + # TODO: always dump save on successful steady-state solution? + printinfo(False) + setvar(self.lastsuccess) + dtreal_orig = deepcopy(bbb.dtreal) + ftol_orig = deepcopy(bbb.ftol) + bbb.dtreal = dtreal_orig/100 + while bbb.dtreal < 1e5: + ftol_old = deepcopy(bbb.ftol) + # Dynamically decreasing fnorm + bbb.pandf1 (-1, -1, 0, bbb.neq, 1., bbb.yl, bbb.yldot) + fnorm_old = (sum((bbb.yldot[:bbb.neq]*\ + bbb.sfscal[:bbb.neq])**2))**0.5 + if bbb.dtreal > 1: + bbb.dtreal = 1e20 + print('\n===== STATIC ITERATION AT DTREAL={:.2e} ====='.format(bbb.dtreal)) + bbb.ftol = max(min(bbb.ftol, 0.01*fnorm_old), 1e-9) + # Take a converging step + if self.exmain_isaborted(): + setvar(self.lastsuccess) + bbb.dtreal = dtreal_orig + bbb.ftol = ftol_orig + return None + # Temporary fix to avoid Segfault when using icntnunk=1! + bbb.ijactot = 2 + bbb.icntnunk = 1 + bbb.ftol = ftol_old + bbb.dtreal = bbb.dtreal*5 + # Assert steady-state + if bbb.iterm != 1: + bbb.dtreal = dtreal_orig + setvar(self.lastsuccess) + bbb.ftol = ftol_orig + print('\n===== STATIC FNRM REDUCTION FAILED =====\n') + return False + self.savesuccess(self.savefname.format('{:.3f}_staticiter'.format(\ + self.lastsuccess).replace('.','p') + )) + print('===== CONVERGED AT STEADY STATE: RETURNING TO MAIN LOOP =====') + bbb.dtreal = dtreal_orig + bbb.ftol = ftol_orig + return True + + def dtsolve(dtdeltafac): + """ Perform a time-dependent solve """ + from uedge import bbb + printinfo(False) + bbb.icntnunk = 0 + # Make backup of original values before entering dt run + # TODO: compress the storing/setting/resetting of vars + tstart_cont = deepcopy(self.tstart) + incpset_cont = deepcopy(bbb.incpset) + itermx_cont = deepcopy(bbb.itermx) + dtreal_cont = deepcopy(bbb.dtreal) + ftol_cont = deepcopy(bbb.ftol) + orig_cont = deepcopy(self.orig) + classvars_cont = deepcopy(self.classvars) + classsetup_cont = deepcopy(self.classsetup) + bbb.incpset = 5 + bbb.itermx = 30 + bbb.dtreal = 1e-5 + bbb.ftol = 1e-8 + abort = False + # Ensure a first time-step can be taken + dtdelta = self.lastsuccess + dtdeltafac/100 + while bbb.iterm != 1: + dtdelta = self.lastsuccess + dtdeltafac/100 + setvar(dtdelta) + if self.exmain_isaborted(): + setvar(self.lastsuccess) + return False + # If the iteration does not succeed, reduce the time-step + if bbb.iterm != 1: + # Advancing less than .5% dt is probably an indication of a failure + if dtdeltafac < 0.5: + abort = True + break + dtdeltafac /=2 + if abort is False: + self.converge(dtreal=dtreal, savedir='.', + savefname=self.savefname.format('{:.3f}_dtrun'.format(\ + dtdelta).replace('.','p')).replace('.hdf5',''), + message='Solving for delta={:.3f}%'.format(dtdelta*100), + ii1max=ii1max, **kwargs) + if bbb.iterm == 1: + self.lastsuccess = dtdelta + # Ensure original values are still being kept track of + self.orig = orig_cont + bbb.incpset = incpset_cont + bbb.itermx = itermx_cont + bbb.dtreal = dtreal_cont + bbb.ftol = ftol_cont + orig_cont = deepcopy(self.orig) + self.classvars = classvars_cont + self.classsetup = classsetup_cont + self.tstart = tstart_cont + bbb.dtreal = dt + if bbb.iterm != 1: + print('=====================================') + print('===== CONTINUATION SOLVE FAILED =====') + print('=====================================') + setvar(self.lastsuccess) + print('Progress upon abortion: {:.3f}%'.format( + self.lastsuccess*100) + ) + return False + else: + return True + + if isinstance(var, str): + if target is None: + raise KeyError("Keyword argument 'target' must be set" + " when evaluating a single variable!" + ) + elif isinstance(var, dict): + if savedir is None: + raise KeyError("Must define save directory name when " + "evalauting multiple variables!" + ) + self.savedir = savedir + if self.savedir is None: + self.savedir = bbb.label[0].decode('UTF-8').strip() + autosavedir = True + else: + autosavedir = False + if len(self.savedir) == 0: + self.savedir = '{}_solve'.format(var) + autosavedir = True + if autosavedir: + if exists(self.savedir): + i=2 + while exists(self.savedir+'_{}'.format(i)): + i += 1 + print('Run-folder {} exists: saving to {}_{}'.format(\ + self.savedir, self.savedir, i) + ) + self.savedir = '{}_{}'.format(self.savedir, i) + else: + print('Saving intermediate saves to {}'.format(self.savedir)) + else: + if exists(self.savedir): + rmtree(self.savedir) + # Finalize savename w/ placeholder + self.savefname = '{}/progress{{}}.hdf5'.format(self.savedir) + makedirs(self.savedir) + # Additional variables to be stored in save files + self.classvars = {} + for dictkeys in ['delta', 'progress', 'delta_fail', 'progress_fail']: + self.classvars[dictkeys] = [] + # Find package variable + # TODO: Use dict to set one or multiple variables instead of string? + + # Set up variables that need to be accessed by helpers + self.dtlim = dtlim + self.cutoff = cutoff + # Set up control flags + start = True + nstaticiter = 0 + iicont = 0 + iicont_fail = 0 + dtiter = 0 + self.isave = saveres + self.saveres = saveres + self.delta = 1/initres + self.lastsuccess = 0 + + # Set up dictionary with variables and targets + self.var = {} + # Construct the variable dictionary based on the user input + if isinstance(var, str): + self.var = { var: { + 'index': index, + 'target': target + }} + else: + self.var = var + for key, subdict in self.var.items(): + if 'target' not in subdict.keys(): + raise KeyError("Target not set for variable {}!".format(key)) + elif 'index' not in subdict.keys(): + self.var[key]['index'] = None + # Complement the user input with gradients for each variable + for key, subdict in self.var.items(): + if isinstance(subdict['target'], list): + subdict['target'] = array(subdict['target']) + for pkg in package(): + pkgobj = packageobject(pkg) + if key in pkgobj.varlist(): + subdict['pkg'] = pkg + subdict['pkgobj'] = pkgobj + if subdict['index'] is not None: + try: + getattr(subdict['pkgobj'], key)[subdict['index']] + except: + raise ValueError('{} does not accommodate requested indices!'.format(key)) + subdict['origvar'] = deepcopy(getvar(key, subdict)) + subdict['deltavar'] = deepcopy(subdict['target'] - subdict['origvar']) + try: + dist = abs(subdict['deltavar']).max() + except: + dist = abs(subdict['deltavar']) + if dist < 1e-10: + raise ValueError('Target equals current value for {}, aborting!'.format(key)) + + self.classsetup = {} + for key, subdict in self.var.items(): + self.setupvars[key] = getattr(subdict['pkgobj'], key) + self.classsetup['initial_{}'.format(key)] = getvar(key, subdict) + self.classsetup['delta_{}'.format(key)] = subdict['deltavar'] + self.classsetup['target_{}'.format(key)] = subdict['target'] + if isinstance(subdict['index'], (tuple, slice)): + self.classsetup['index_{}'.format(key)] = str(subdict['index']).encode("ascii", "ignore") + elif subdict['index'] is not None: + self.classsetup['index_{}'.format(key)] = subdict['index'] + else: + self.classsetup['index_{}'.format(key)] = False + + # Record original solver settings + self.orig = {} + for ovar in ['itermx', 'dtreal', 'icntnunk', 'ftol', 'incpset', + 'ismmaxuc', 'mmaxu' + ]: + self.orig[ovar] = deepcopy(getattr(bbb, ovar)) + # Take the initial time-step + changevar() + # Set remaining solver settings + bbb.dtreal = dt + bbb.ftol = ftol + bbb.incpset = incpset + bbb.itermx = itermx +# bbb.ismmaxuc = 0 +# bbb.mmaxu = 70 + if (bbb.iterm == 1) and (bbb.ijactot > 0): + self.message('Initial successful time-step exists', separator='') + else: + self.message('Need to take initial step with Jacobian; ' + \ + 'trying to do here', seppad='*') + printinfo() + # Ensure preconditioner is calculated + bbb.icntnunk = 0 + # Take timestep and see if abort requested + if self.exmain_isaborted(): + setvar(self.lastsuccess) + return + # Increase time + # Verify time-step was successful + if (bbb.iterm != 1): + self.restorevalues() + setvar(self.lastsuccess) + self.message('Error: converge an initial time-step first; then ' + \ + 're-execute command', seppad='*') + return + # Start outer loop + while True: + # Ensure we don't exceed the target value + # Re-eval preconditioner: omit first step + if start is False: + bbb.icntnunk = 0 + else: + bbb.icntnunk = 1 + start = False + # Reset counters + iicont = 0 + iicont_fail = 0 + # Start inner loop + while iicont < iicont_max: + print('===== MAIN LOOP {}/{} ====='.format(iicont+1, iicont_max)) + printinfo() + # TODO: how to get intial fnrm rather than last fnorm? + fnorm_old = (sum((bbb.yldot[:bbb.neq]*\ + bbb.sfscal[:bbb.neq])**2))**0.5 + bbb.ftol = min(max(fnorm_old/100, 1e-6), ftol) + # Take step and check whether abort is requested + if self.exmain_isaborted(): + setvar(self.lastsuccess) + return + bbb.ftol = ftol + bbb.icntnunk = 1 + # Failure to take time-step + if bbb.iterm != 1: + # Check whether an excorbant number of iterations are necessary to converge: fail if yes + self.classvars['delta_fail'].append(self.delta) + self.classvars['progress_fail'].append(self.lastsuccess) + if isdeltatoosmall(): + print('=====================================') + print('===== CONTINUATION SOLVE FAILED =====') + print('=====================================') + print('Last successful step for {}: {:.4e}'.format(\ + self.var, self.lastsuccess) + ) + return + # Start trying to reset convergence + else: + # Reset success counter + iicont = 0 + iicont_fail += 1 + # Delta has been reduced too many times: do some moves to knock the case loose + if iicont_fail >= iicont_fail_max: + print('===== RECOVERY LOOP ENTERED =====') + # Iterate the case at the last successful + bbb.icntnunk = 0 + if dostaticiter(): + staticstatus = staticiter() + # Flag for abortion + if staticstatus is None: + return + elif staticstatus is False: + dtcall=True + else: + dtcall=False + nstaticiter += 1 + # Enter loop for time-dependent simulations + if (bbb.iterm != 1) or (dtcall is True) or (nstaticiter==staticiter_max): + print('===== ENTER TIME-DEPENDENT SOLVE =====') + if dtsolve(dtdeltafac) is False: + return + else: + # Force-break inner loop + iicont = 1e5 + dtcall = False + bbb.icntnunk = 0 + if issuccess(): + return + changevar() + # Initial fnrm successfully reduced, continue increasing var + else: + iicont_fail = 0 + changevar() + # Try to attain convergence again by decreasing change in var + else: + print('\n===== FAIL {}/{} ====='.format(iicont_fail, iicont_fail_max)) + print('===== DECREASE DELTA AND TRY AGAIN =====') + # Reset to last successful step + setvar(self.lastsuccess) + # Decrease change in variable + self.delta /= deltafac + # Set to new, decreased change + changevar() + bbb.icntnunk = 0 + # Success: update counter and last successful value + else: + self.classvars['delta'].append(self.delta) + self.classvars['progress'].append(self.lastsuccess) + if issuccess(): + return + # Check whether to increment time-step + if iicont == iicont_max-1: + print('\n===== INNER LOOP COMPLETED: ADVANCING DELTA =====') + nstaticiter = 0 + self.delta *= 1.1*deltafac + else: + print('\n===== SUCCESS: ADVANCING VARIABLE =====') + iicont += 1 + changevar() + print('EXITED LOOP: YOU SHOULD NOT BE HERE...') + return def rundt(**kwargs): runcase=UeRun() diff --git a/pyscripts/ruthere.py b/pyscripts/ruthere.py deleted file mode 100755 index 7b42c2d3..00000000 --- a/pyscripts/ruthere.py +++ /dev/null @@ -1,66 +0,0 @@ -from numpy import * -import time -import signal -import sys - -############################################################################# -# From Dave Grote: -# --- Setup signal handler to capture Control-C -# --- To use, first call arminterrupt(). Then at the place where the interrupt -# --- is allowed, call ruthere(). This will raise a KeyboardInterrupt if -# --- Control-C had been pressed. -# --- When a interrupt request is received, all this handler does is set a -# --- flag to that effect. Then, a subsequent call to ruthere will check -# --- that flag, and if set, raise an exception. This allows a graceful -# --- stop with the current time step completed. -# --- Set the following two in case ruthere is called before arminterrupt. -_defaultcontrolC = signal.getsignal(signal.SIGINT) -_controlCrecieved = False -savetracebacklimit = 0 - -def _handlecontrolC(signum, frame): - global _controlCrecieved - _controlCrecieved = True - -def ruthere(reset=True): - """ -Checks if an interrupt was requested (usually control-C). If so, then raise -an exception. If reset is True, restore the original interrupt handler so that the -calling code does not have to, and so that, if there is an exception, it gets -restored (since the calling code is not returned to). - """ - global _controlCrecieved - global _defaultcontrolC - global savetracebacklimit - if _controlCrecieved: - if reset: - signal.signal(signal.SIGINT, _defaultcontrolC) - _controlCrecieved = False - raise KeyboardInterrupt("Interrupt requested") - -def arminterrupt(): - global _controlCrecieved - global _defaultcontrolC - global savetracebacklimit - _controlCrecieved = False - _defaultcontrolC = signal.getsignal(signal.SIGINT) - try: - savetracebacklimit = sys.tracebacklimit - except: - savetracebacklimit = None - signal.signal(signal.SIGINT, _handlecontrolC) - sys.tracebacklimit = 0 - -def disarminterrupt(): - global _defaultcontrolC - global savetracebacklimit - signal.signal(signal.SIGINT, _defaultcontrolC) - sys.tracebacklimit = savetracebacklimit - - -#========================================================================= - -arminterrupt() - - - diff --git a/pyscripts/uedgeplots.py b/pyscripts/uedgeplots.py index 502dfb86..ce1aa8db 100755 --- a/pyscripts/uedgeplots.py +++ b/pyscripts/uedgeplots.py @@ -44,13 +44,16 @@ def plotmesh(ixmin=None, ixmax=None, iymin=None, iymax=None, The plot axis limits may be specified with r_rmin,r_max,z_min,z_max. """ - zrefl = com.zm - zlim = com.ylim - zreflbdry = com.zbdry - if str(com.geometry) == str([b'uppersn ']): - zrefl = 2.0 * com.zmid - com.zm - zlim = 2.0 * com.zmid - com.ylim - zreflbdry = 2.0 * com.zmid - com.zbdry + try: + zrefl = com.zm + zlim = com.ylim + zreflbdry = com.zbdry + if str(com.geometry) == str([b'uppersn ']): + zrefl = 2.0 * com.zmid - com.zm + zlim = 2.0 * com.zmid - com.ylim + zreflbdry = 2.0 * com.zmid - com.zbdry + except: + pass if ixmin == None: ixmin = com.nxomit @@ -112,13 +115,16 @@ def plotanymesh(verts, r_min=None, r_max=None, z_min=None, z_max=None, title=Non The plot axis limits may be specified with r_rmin,r_max,z_min,z_max. """ - zrefl = com.zm - zlim = com.ylim - zreflbdry = com.zbdry - if str(com.geometry) == str([b'uppersn ']): - zrefl = 2.0 * com.zmid - com.zm - zlim = 2.0 * com.zmid - com.ylim - zreflbdry = 2.0 * com.zmid - com.zbdry + try: + zrefl = com.zm + zlim = com.ylim + zreflbdry = com.zbdry + if str(com.geometry) == str([b'uppersn ']): + zrefl = 2.0 * com.zmid - com.zm + zlim = 2.0 * com.zmid - com.ylim + zreflbdry = 2.0 * com.zmid - com.zbdry + except: + pass if r_min == None: r_min = np.min(verts[:,:,:,0]) @@ -181,13 +187,16 @@ def plotmeshval(val, ixmin=None, ixmax=None, iymin=None, iymax=None, The plot axis limits may be specified with r_rmin,r_max,z_min,z_max. """ - zrefl = com.zm - zlim = com.ylim - zreflbdry = com.zbdry - if str(com.geometry) == str([b'uppersn ']): - zrefl = 2.0 * com.zmid - com.zm - zlim = 2.0 * com.zmid - com.ylim - zreflbdry = 2.0 * com.zmid - com.zbdry + try: + zrefl = com.zm + zlim = com.ylim + zreflbdry = com.zbdry + if str(com.geometry) == str([b'uppersn ']): + zrefl = 2.0 * com.zmid - com.zm + zlim = 2.0 * com.zmid - com.ylim + zreflbdry = 2.0 * com.zmid - com.zbdry + except: + pass if ixmin == None: ixmin = com.nxomit @@ -263,13 +272,16 @@ def plotanymeshval(verts,z, r_min=None, r_max=None, z_min=None, z_max=None, titl The plot axis limits may be specified with r_rmin,r_max,z_min,z_max. """ - zrefl = com.zm - zlim = com.ylim - zreflbdry = com.zbdry - if str(com.geometry) == str([b'uppersn ']): - zrefl = 2.0 * com.zmid - com.zm - zlim = 2.0 * com.zmid - com.ylim - zreflbdry = 2.0 * com.zmid - com.zbdry + try: + zrefl = com.zm + zlim = com.ylim + zreflbdry = com.zbdry + if str(com.geometry) == str([b'uppersn ']): + zrefl = 2.0 * com.zmid - com.zm + zlim = 2.0 * com.zmid - com.ylim + zreflbdry = 2.0 * com.zmid - com.zbdry + except: + pass if r_min == None: r_min = com.rm.min() @@ -333,13 +345,16 @@ def mkdensityfile(filename, ival, renmin=None, renmax=None, samples=[500, 500, 5 about 1cm r resolution and .6cm in z. """ - zrefl = com.zm - zlim = com.ylim - zreflbdry = com.zbdry - if str(com.geometry) == str([b'uppersn ']): - zrefl = 2.0 * com.zmid - com.zm - zlim = 2.0 * com.zmid - com.ylim - zreflbdry = 2.0 * com.zmid - com.zbdry + try: + zrefl = com.zm + zlim = com.ylim + zreflbdry = com.zbdry + if str(com.geometry) == str([b'uppersn ']): + zrefl = 2.0 * com.zmid - com.zm + zlim = 2.0 * com.zmid - com.ylim + zreflbdry = 2.0 * com.zmid - com.zbdry + except: + pass if renmin == None: renmin = np.min(ival) diff --git a/pyscripts/uexec.py b/pyscripts/uexec.py deleted file mode 100755 index e5e51522..00000000 --- a/pyscripts/uexec.py +++ /dev/null @@ -1,30 +0,0 @@ - - -import sys -try: - from importlib import reload,import_module -except: - from importlib import import_module - - -import builtins - - -def uexec(mname,returns=globals()): - if mname in sys.modules: - _m = reload(sys.modules[mname]) - else: - _m = import_module(mname) - - # is there an __all__? if so respect it - if "__all__" in _m.__dict__: - names = _m.__dict__["__all__"] - else: - # otherwise we import all names that don't begin with _ - names = [x for x in _m.__dict__ if not x.startswith("_")] - - # now drag them in - for k in names: - #print k,getattr(_m,k) - returns[k] = getattr(_m,k) - diff --git a/setup.py b/setup.py index d8cfcb8a..7a9ce7c4 100755 --- a/setup.py +++ b/setup.py @@ -187,7 +187,7 @@ def makeobjects(pkg): cmdclass={'build': uedgeBuild, 'clean': uedgeClean}, test_suite="pytests", - install_requires=['forthon', 'easygui'], + install_requires=['forthon'], # note that include_dirs may have to be expanded in the line above classifiers=['Programming Language :: Python', 'Programming Language :: Python :: 3'] diff --git a/wdf/wdf.m b/wdf/wdf.m index 7fc10e0b..0e98847e 100755 --- a/wdf/wdf.m +++ b/wdf/wdf.m @@ -5,7 +5,7 @@ Use(Xpoint_indices) # ixpt1,ixpt2 from com package Use(Share) # geometry from com package Use(Dimwdf) # nptsvb,nptshb,nptskb,nptsw,npsegxz -Use(Auxw) # ixpt1b,ixpt2b,ixtop1b,ixtop2b, +Use(Auxw) # ixpt1b_wdf,ixpt2b_wdf,ixtop1b,ixtop2b, # nohzsb,novzsb,nosegsxzb Use(Options) # iswdfon @@ -21,10 +21,10 @@ call xerrab("") c allocate arrays for DEGAS - ixpt1b = ixpt1(1) - ixtop1b = ixpt1b + (ixpt2(1)-ixpt1(1))/2 + ixpt1b_wdf = ixpt1(1) + ixtop1b = ixpt1b_wdf + (ixpt2(1)-ixpt1(1))/2 ixtop2b = ixtop1b + 1 - ixpt2b = ixpt2(1) + ixpt2b_wdf = ixpt2(1) if (geometry .eq. "dnbot") then ixtop1b=ixtop1b-1 # we omit guard cells at the midplane ixtop2b=ixtop2b+1 # of up/down symmetric double nulls @@ -151,7 +151,7 @@ subroutine rdgrd3 (iunit) Use(RZ_grid_info) # rm,zm from com package Use(Linkbbb) # nxbbb,nybbb from com package Use(Dimwdf) -Use(Auxw) # ixpt1b,ixpt2b,ixtop1b,ixtop2b +Use(Auxw) # ixpt1b_wdf,ixpt2b_wdf,ixtop1b,ixtop2b Use(Linkgrd) Use(Degas1) Use(Degas2) @@ -167,7 +167,7 @@ subroutine rdgrd3 (iunit) i = 1 gridx(j,i,1) = rm(ixtop1b,iy,4) - rgrid1w gridz(j,i,1) = zm(ixtop1b,iy,4) - do ix = ixtop1b, ixpt1b+1, -1 + do ix = ixtop1b, ixpt1b_wdf+1, -1 i = i + 1 gridx(j,i,1) = rm(ix,iy,3) - rgrid1w gridz(j,i,1) = zm(ix,iy,3) @@ -176,9 +176,9 @@ subroutine rdgrd3 (iunit) gridx(j,i,1) = gridx(j,i-1,1) gridz(j,i,1) = gridz(j,i-1,1) i = i + 1 - gridx(j,i,1) = rm(ixpt1b,iy,4) - rgrid1w - gridz(j,i,1) = zm(ixpt1b,iy,4) - do ix = ixpt1b, 1, -1 + gridx(j,i,1) = rm(ixpt1b_wdf,iy,4) - rgrid1w + gridz(j,i,1) = zm(ixpt1b_wdf,iy,4) + do ix = ixpt1b_wdf, 1, -1 i = i + 1 gridx(j,i,1) = rm(ix,iy,3) - rgrid1w gridz(j,i,1) = zm(ix,iy,3) @@ -198,7 +198,7 @@ subroutine rdgrd3 (iunit) i = 1 gridx(j,i,1) = rm(ixtop2b,iy,3) - rgrid1w gridz(j,i,1) = zm(ixtop2b,iy,3) - do ix = ixtop2b, ixpt2b + do ix = ixtop2b, ixpt2b_wdf i = i + 1 gridx(j,i,1) = rm(ix,iy,4) - rgrid1w gridz(j,i,1) = zm(ix,iy,4) @@ -207,9 +207,9 @@ subroutine rdgrd3 (iunit) gridx(j,i,1) = gridx(j,i-1,1) gridz(j,i,1) = gridz(j,i-1,1) i = i + 1 - gridx(j,i,1) = rm(ixpt2b+1,iy,3) - rgrid1w - gridz(j,i,1) = zm(ixpt2b+1,iy,3) - do ix = ixpt2b+1, nxbbb + gridx(j,i,1) = rm(ixpt2b_wdf+1,iy,3) - rgrid1w + gridz(j,i,1) = zm(ixpt2b_wdf+1,iy,3) + do ix = ixpt2b_wdf+1, nxbbb i = i + 1 gridx(j,i,1) = rm(ix,iy,4) - rgrid1w gridz(j,i,1) = zm(ix,iy,4) @@ -225,11 +225,11 @@ subroutine rdgrd3 (iunit) endwhile c Finally, the "magnetic axis" for core and private flux regions j = nybbb + 2 - do i = 1, ixtop1b-ixpt1b+2 + do i = 1, ixtop1b-ixpt1b_wdf+2 gridx(j,i,1) = rmagxw - rgrid1w gridz(j,i,1) = zmagxw enddo - do i = ixtop1b-ixpt1b+3, nptshb + do i = ixtop1b-ixpt1b_wdf+3, nptshb gridx(j,i,1) = 0.5*(gridx(j-1,nptshb,1)+gridx(j+1,nptshb,1)) gridz(j,i,1) = 0.5*(gridz(j-1,nptshb,1)+gridz(j+1,nptshb,1)) enddo @@ -277,13 +277,13 @@ subroutine rdgrd3 (iunit) xwall(itot,2)=gridx(j,i,1) zwall(itot,2)=gridz(j,i,1) j=j-1 # jump to inboard half of p.f. - do i=ixtop1b+3,ixtop1b-ixpt1b+3,-1 + do i=ixtop1b+3,ixtop1b-ixpt1b_wdf+3,-1 itot=itot+1 xwall(itot,2)=gridx(j,i,1) zwall(itot,2)=gridz(j,i,1) enddo j=j+2 # jump to outboard half of p.f. - do i=ixpt2b-ixtop2b+4,nxbbb-ixtop2b+4 + do i=ixpt2b_wdf-ixtop2b+4,nxbbb-ixtop2b+4 itot=itot+1 xwall(itot,2)=gridx(j,i,1) zwall(itot,2)=gridz(j,i,1) @@ -302,13 +302,13 @@ subroutine rdgrd3 (iunit) c Third, clockwise around the core itot=0 j=nybbb+1 - do i=ixtop1b-ixpt1b+2,1,-1 + do i=ixtop1b-ixpt1b_wdf+2,1,-1 itot=itot+1 xwall(itot,3)=gridx(j,i,1) zwall(itot,3)=gridz(j,i,1) enddo j=j+2 - do i=1,ixtop1b-ixpt1b+2 + do i=1,ixtop1b-ixpt1b_wdf+2 itot=itot+1 xwall(itot,3)=gridx(j,i,1) zwall(itot,3)=gridz(j,i,1) @@ -365,7 +365,7 @@ call xerrab("**** bbb-wdf file not found") Use(Xpoint_indices) # from com package Use(Linkbbb) Use(Dimwdf) -Use(Auxw) # ixpt1b,ixpt2b,ixtop1b,ixtop2b,nosegsxzb,novzsb +Use(Auxw) # ixpt1b_wdf,ixpt2b_wdf,ixtop1b,ixtop2b,nosegsxzb,novzsb Use(Linkgrd) Use(Degas1) Use(Degas2) @@ -376,16 +376,16 @@ call xerrab("**** bbb-wdf file not found") do ix=1,nxbbb do iy=1,nybbb c translate UEDGE cell index (ix,iy) to DEGAS cell index (iv,jh) - if ((1 .le. ix) .and. (ix .le. ixpt1b)) then + if ((1 .le. ix) .and. (ix .le. ixpt1b_wdf)) then jh = nybbb + 1 - iy iv = ixtop1b + 3 - ix - elseif ((ixpt1b+1 .le. ix) .and. (ix .le. ixtop1b)) then + elseif ((ixpt1b_wdf+1 .le. ix) .and. (ix .le. ixtop1b)) then jh = nybbb + 1 - iy iv = ixtop1b + 1 - ix - elseif ((ixtop2b .le. ix) .and. (ix .le. ixpt2b)) then + elseif ((ixtop2b .le. ix) .and. (ix .le. ixpt2b_wdf)) then jh = nybbb + 2 + iy iv = ix - ixtop2b + 1 - elseif ((ixpt2b+1 .le. ix) .and. (ix .le. nxbbb)) then + elseif ((ixpt2b_wdf+1 .le. ix) .and. (ix .le. nxbbb)) then jh = nybbb + 2 + iy iv = ix - ixtop2b + 3 endif @@ -410,12 +410,12 @@ c translate UEDGE plate index (nxm,iy) to DEGAS wall segment index (iseg) c translate UEDGE private flux wall index (ix) to DEGAS private region c wall segment index (iseg) on wall number 2 - do ix=1,ixpt1b + do ix=1,ixpt1b_wdf iseg = 1 + ix currxzt(iseg,1,2,1) = - abs(fngysibbb(ix))# NOTE negative for puffing enddo - do ix=ixpt2b+1,nxbbb - iseg = 1 + ixpt1b + 1 + (ix-ixpt2b) + do ix=ixpt2b_wdf+1,nxbbb + iseg = 1 + ixpt1b_wdf + 1 + (ix-ixpt2b_wdf) currxzt(iseg,1,2,1) = - abs(fngysibbb(ix)) enddo @@ -438,7 +438,7 @@ c wall segment index (iseg) on wall number 1 subroutine defaultz integer jb,jj,ib,ii Use(Dimwdf) -Use(Auxw) # ixpt1b +Use(Auxw) # ixpt1b_wdf Use(Degas1) Use(Degas2) Use(Eqdsk) @@ -468,7 +468,7 @@ c wall segment index (iseg) on wall number 1 novbs = ib novzs = novbs - 1 - ivnull = ixtop1b - ixpt1b + 2 # vertical index of non-physical zone + ivnull = ixtop1b - ixpt1b_wdf + 2 # vertical index of non-physical zone # between core and private flux c.....set kzone arrays diff --git a/wdf/wdf.v b/wdf/wdf.v index b4c1036a..b2c5417d 100755 --- a/wdf/wdf.v +++ b/wdf/wdf.v @@ -51,10 +51,10 @@ nixw integer ***** Auxw: # some parameters used to transform UEDGE indices to DEGAS indices -ixpt1b integer # ix index of pf cell just below x-point on inboard half +ixpt1b_wdf integer # ix index of pf cell just below x-point on inboard half ixtop1b integer # ix index of cell at top of inboard half mesh ixtop2b integer # ix index of cell at top of outboard half mesh -ixpt2b integer # ix index of core cell just above x-point on outboard half +ixpt2b_wdf integer # ix index of core cell just above x-point on outboard half nohzsb integer # number of DEGAS horizontal zones novzsb integer # number of DEGAS vertical zones nosegsxzb integer # number of DEGAS wall segments @@ -258,9 +258,9 @@ sndspd real # plasma flow mach number where sound speed = sqrt((te+ti)/mi) t0puff(npns) _real [eV] /.025/ # energy of puffed source of neutrals -te0 real [eV] +te0_wdf real [eV] # electron temperature in the case of a constant plasma -ti0(npis) _real [eV] +ti0_wdf(npis) _real [eV] # ion temperature in the case of a constant plasma wtmin0 real # minimum test flight weight