C -*- Mode: Fortran; -*- c----------------------------------------------------------------- c Ravi Samtaney c Copyright 2004 c Princeton Plasma Physics Laboratory c All Rights Reserved c----------------------------------------------------------------- c $Log: init.F,v $ c Revision 1.2 2005/10/08 17:26:31 samtaney c Same equilibrium as in ICNSP paper. This one shows the correct c radial transport of mass with pellet injection. c c Revision 1.1.1.1 2005/09/12 17:28:10 samtaney c Original source. c c Revision 1.2 2004/09/09 21:02:24 samtaney c Added initial equilibrium field - analytical soln to G-S as in Pellet case. c c Revision 1.1 2004/06/18 22:29:57 samtaney c Vanilla initial conditions - just for testing. c c----------------------------------------------------------------- c Subroutine to set up the initial conditions c--------------------------------------------------------------------------- c subroutine InitialConditions(ux,eqsrc) use mesh use mesh_common use properties use iounits use BoundaryPressTemp use MagAxisToroidalField use PelletParams use OptionFlags implicit none double precision:: ux(IXLO:IXHI,IYLO:IYHI,IZLO:IZHI,nvar) double precision:: eqsrc(1:nxlocal,1:nylocal,1:nzlocal,nvar) c c=======Declarations========= real*8 rho,press, Bi, Bj, Bk, u,v,w c integer:: i,j,k character*24 fdate real*8 kx,ky,lamda real*8 p1,rho1,u1,amach,gas1,gas2,amp double precision:: RMajor, a, R0 double precision:: eps, aspect, beta double precision:: alpha, pi, alpha0 double precision:: pio2ka double precision:: dx,xin,x,y double precision:: fun(ixlo:ixhi) double precision:: dfdxi(ixlo:ixhi) double precision:: d2fdxi2(ixlo:ixhi) double precision:: fmin, g0 double precision:: psi, dpsidx,dpsidy, pisq double precision:: d2psidx2,d2psidy2 double precision:: zeta,coszeta, sinzeta, g0 double precision:: tempAxis integer ig, jg c namelist/molwts/wmol namelist/gammas/gamma namelist/viscous/Reynolds,Lundquist, Prandtl c open(16,file='prop.inp',form='formatted') read(16,molwts) write(6,molwts) read(16,gammas) write(6,gammas) read(16,viscous) write(6,viscous) close(16) rgas=1.D0/wmol c pi=4.D0*datan(1.D0) c Scale factor alpha=1.0D0 c RMajor = radius of inner boundary RMajor=3.D0 RMajor=9.D0/(4.D0*datan(1.D0))-1.D0 R0=RMajor a=1.D0 aspect=1.D0 c For torus with aspect=1 RMajor=R0=3, a=1 c alpha0=0.2757368396764204D0 c For torus with aspect=1 with R0=9/pi-1, a=1 alpha0=0.4815088503327338D0 c Radial function parameters eps=a/(R0+a) c eps=0.D0 beta=0.25D0*pi*pi/aspect**2 Btaxis=1.D0*(RMajor+a) g0=Btaxis psi0=1.5D0*a/(alpha0*pi**2)* & (g0)/(RMajor+a) c g0=0.D0 c lamda=pi*pi*alpha0 pio2ka=0.5D0*pi/(aspect*a) pfac=0.5D0*alpha0*pi*pi/(a*a*(R0+a)*(R0+a)) c Reference quantities wmol=2.D0 RUniversal=8314.5D0 c rgas=RUniversal/wmol c This is a non-dim gas constant - Eos is p=rho T rgas=1.D0 c Permeability of free space in Henry/m or kg-m/Coulomb^2 mu0=4.D0*pi * 1.D-07 c Reference Magnetic field in Tesla (kg/s/coulomb) Bref=1.D0 c Bref=1.D0 c CDXU c Bref=0.375D0 c DIII-D c Bref=2.0D0 c ITER c Bref=5.0D0 c Bref=BToroidalMax c c Alternative pref is Bref^2/mu0 (in N/m^2) pref=Bref**2/mu0 c Reference no. of molecules /m^3 n0ref=1.5D19 c Mass of each proton in Kg Mi=1.6605 D-27 c Reference density rhoRef=n0ref*Mi*wmol c rhoRef=1.6605*wmol*1.D-08 rhoRef=Mi*wmol*n0Ref c Reference temperature in K TRef=pref/(rhoRef*RUniversal/wmol) c Convert K to eV K2eV=1.1604D04 c Reference temperature in eV TRefeV=Tref/K2eV c c Reference velocity is the alfven velocity defined c as Bref/sqrt(rhoRef*mu0) velRef=Bref/dsqrt(rhoRef*mu0) c c Reference length scale is 1m lref=1.D0 lref=0.26D0 c c Reference time scale tauRef tauRef=lref/velRef c c Initial radius of pellet in m rpinit=0.001D0 c rpinit=2.756654269586026D-04 c rpinit=0.00012D0 c rpinit=0.0D0 c Radius - nondimensional rpinit=rpinit/lref rp=rpinit c dx=2.D0*a/nx c c Up is the radial, axial and azimuthal velocities c up(1)=24000.D0/velRef up(1)=(LaunchType)*1000.D0/velRef up(2)=0.D0 up(3)=0.D0 c xp is the initial radial, axial and azimuthal positions c For square torus if(LaunchType.eq.1) then xpinit(1)=Rmajor+2*dx else if (LaunchType.eq.-1) then xpinit(1)=Rmajor+2.D0*a-2*dx endif xpinit(2)=0.D0 xpinit(3)=0.D0 xp=xpinit c For DIII-D type params c$$$ if(LaunchFlag.eq.1) then c$$$ xpinit(1)=1.3D0 c$$$ else if (LaunchFlag.eq.-1) then c$$$ xpinit(1)=2.2D0 c$$$ endif c$$$ xpinit(2)=0.D0 c$$$ xpinit(3)=0.D0 c$$$ xp=xpinit c call iglobal(iprocx,1,ig) call jglobal(iprocy,1,jg) c c allocate(Rtrajectory(nx+1)) if(up(1).gt.0.D0) then do i=1,nx+1,1 Rtrajectory(i)=RMajor+(ig+i-2)*dx enddo else do i=1,nx+1,1 Rtrajectory(i)=RMajor+2.D0*a - (ig+i-2)*dx enddo endif c c Ambient pressure c pambient=0.05D0 pambient=0.03D0/mu0*alpha*alpha c if(iprocx.eq.1.and.iprocy.eq.1.and.iprocz.eq.1) then write(ihis,*) 'Reference Quantities' write(ihis,*) 'Density (kg/m^3)', rhoRef write(ihis,*) 'Number density (/m^3)',n0Ref write(ihis,*) 'Pressure (Pa)', pRef write(ihis,*) 'Temperature (K)', TRef write(ihis,*) 'Temperature (eV)', TRefeV write(ihis,*) 'Magnetic field (T)', Bref write(ihis,*) 'Velocity (m/s)', velRef write(ihis,*) 'Energy (J/m^3)', rhoRef*velRef**2 write(ihis,*) 'Length (m)', lref write(ihis,*) 'time (s)', tauRef write(ihis,*) write(ihis,*) 'Ambient pressure (Pa)', pambient write(ihis,*) write(ihis,*) 'Plasma gamma', gamma write(ihis,*) 'Universal gas constant (J/kmol/K)', RUniversal write(ihis,*) 'Molecular weight', wmol write(ihis,*) 'Gas constant (nondim)', rgas write(ihis,*) 'EoS: p=rho T' write(ihis,*) 'Reynolds Number:', Reynolds write(ihis,*) 'Lundquist Number:', Lundquist write(ihis,*) 'Prandtl Number:', Prandtl write(ihis,*) write(ihis,*) 'Pellet Parameters' write(ihis,*) 'Initial radius (m)', rpinit*lref write(ihis,*) 'Initial radial position(m)', xpinit(1)*lref write(ihis,*) 'Initial radial velocity (m/s)', up(1)*velRef write(ihis,*) 'Initial radial velocity (nondim)', up(1) if(up(1).lt.0.D0) then write(ihis,*) 'LFS pellet launch' else write(ihis,*) 'HFS pellet launch' endif do i=1,nx+1,1 write(ihis,*) 'Rtrajectory=',i,Rtrajectory(i) enddo write(ihis,*) endif c dx=2.D0*a/nx write(6,*) 'lamda=',lamda write(6,*) 'Scaled alpha=',alpha write(6,*) 'Scaled beta=',beta write(6,*) 'eps=',eps write(6,*) 'Initial conditions',RMajor rho=1.D0 press=pambient/pRef u=0.D0 v=0.D0 w=0.D0 c Bi=0.00D0 Bj=0.0D0 Bk=0.D0 c call iglobal(iprocx,1,ig) call jglobal(iprocy,1,jg) do i=ixlo,ixhi,1 x = RMajor+(i +ig -1- 0.5D0)*dx xin=(x-R0-a)/a call yfun(fun(i),dfdxi(i),d2fdxi2(i),alpha0,beta, & lamda,eps,xin) write(6,*) 'INIT',i,xin,fun(i), dfdxi(i) enddo do k=izlo,izhi,1 do j=iylo,iyhi,1 c y = -a+(j + jg -1 - 0.5D0)*dx do i=ixlo,ixhi,1 c x=RMajor+(i+ig -1 -0.5D0)*dx y = zc(i,j) x=Rc(i,j) c write(6,*) 'INIT',i,j,Rc(i,j), zc(i,j), x,y zeta=pio2ka*y coszeta=dcos(zeta) sinzeta=dsin(zeta) psi=psi0*fun(i)*coszeta dpsidx=1.D0/a*dfdxi(i)*coszeta d2psidx2=1.D0/a/a*d2fdxi2(i)*coszeta dpsidy=-pio2ka*fun(i)*sinzeta d2psidy2=-pio2ka**2*fun(i)*coszeta press=pambient/pRef+pfac*psi*psi*alpha*alpha bi=-dpsidy/(x)*psi0*alpha bj=dpsidx/(x)*psi0*alpha c press=0.1D0 c bi=0.D0 c bj=0.D0 c bk=Btaxis/x bk=0.D0 if(i.lt.nx/2) then ux(i,j,k,1)=rho ux(i,j,k,2)=rho*u ux(i,j,k,3)=rho*v ux(i,j,k,4)=rho*w ux(i,j,k,5)=Bi ux(i,j,k,6)=Bj ux(i,j,k,7)=Bk ux(i,j,k,8)=press/(gamma-1.D0)+ & 0.5*rho*(u*u+v*v+w*w)+0.5D0*(bi*bi+bj*bj+bk*bk) else ux(i,j,k,1)=1.D0*rho ux(i,j,k,2)=rho*u ux(i,j,k,3)=rho*v ux(i,j,k,4)=rho*w ux(i,j,k,5)=Bi ux(i,j,k,6)=Bj ux(i,j,k,7)=Bk ux(i,j,k,8)=1.D0*press/(gamma-1.D0)+ & 0.5*rho*(u*u+v*v+w*w)+0.5D0*(bi*bi+bj*bj+bk*bk) endif write(6,*) 'INIT EQCHECK',i,j,(x-R0-a)/a, & psi0*d2psidx2+psi0*d2psidy2 & -psi0*dpsidx/x & +x*x*2.D0*pfac*psi enddo enddo enddo eqsrc=0.D0 c Order of calling computeRHS and EMCT is important. c computeRHS sets 5:7 terms of RHS - and these are c overwritten by EMCT. In the time marching also c these are important to be called in this order. if(eqsrcFlag.eq.0) then return else if(eqsrcFlag.eq.1) then call ComputeRHS(ux,eqsrc,0.D0) if(emCTFlag.eq.1) then call EMCT(ux,eqsrc,0.D0) endif else if(eqsrcFlag.eq.2) then call ComputeRHSInit(ux,eqsrc,0.D0) if(emCTFlag.eq.1) then call EMCT(ux,eqsrc,0.D0) endif endif c c Plasma temp on axis = p_axis/pref * rho (nondim) and c multiplied by TRefeV gives it in ev TempAxis=(pambient)/pref & *TRefeV c pbdry=pambient/pRef c Taking density as unity TempBdry=pbdry c c Temporary assignment of TempAblation tempAblation=0.D0 c if(iprocx.eq.1.and.iprocy.eq.1.and.iprocz.eq.1) then write(ihis,*) 'Plasma Beta (based on mag axis values)=', & 2.D0*mu0*(pambient)/ & (Btaxis/(rMajor+a))**2 write(ihis,*) 'Boundary Temperature=', & TempBdry endif c return end c----------------------------------------------------------------------- subroutine yfunOLD(y,dy,d2y,alpha,beta,lamda,eps,x) integer ncoef parameter (ncoef=100) c real*8 c(-5:ncoef) real*8 dcda(-5:ncoef) real*8 alpha,beta,lamda,eps real*8 y,dy,x,d2y integer n c c(-5)=0.D0 c(-4)=0.D0 c(-3)=0.D0 c(-2)=0.D0 c(-1)=0.D0 c(0)=0.D0 c c(1)=-1.D0 c(1)=-4.D0*datan(1.D0) dcda(0)=0.D0 dcda(1)=0.D0 c do n=2,ncoef,1 c(n)=-1.D0/(n*(n-1))*( & eps*(n-1)*(n-3)*c(n-1)+ & (lamda)*(c(n-2)+ & 3.D0*eps*c(n-3)+ & 3.D0*eps**2*c(n-4)+ & eps**3*c(n-5))- & beta*(c(n-2)+eps*c(n-3))) enddo c y=0.D0 c dy=ncoef*c(ncoef)*x c d2y=ncoef*(ncoef-1)*c(ncoef)*x dy=0.D0 d2y=0.D0 c do n=ncoef,1,-1 y=(y+(c(n)))*x enddo do n=ncoef,2,-1 dy=(dy+(n*c(n)))*x enddo do n=ncoef,3,-1 d2y=(d2y+(n*(n-1)*c(n)))*x enddo dy=(dy+(c(1))) d2y=(d2y+(2*c(2))) return end c----------------------------------------------------------------------- c----------------------------------------------------------------------- subroutine yfun(y,dy,d2y,alpha,beta,lamda,eps,x) integer ncoef parameter (ncoef=100) c real*8 c(-5:ncoef) real*8 dcda(-5:ncoef) real*8 alpha,beta,lamda,eps real*8 y,dy,x,d2y integer n c c(-5)=0.D0 c(-4)=0.D0 c(-3)=0.D0 c(-2)=0.D0 c(-1)=0.D0 c(0)=-1.D0 c c(1)=-4.D0*datan(1.D0) c These R0s are the centerline radius c For R0=9/pi c(1)=-0.6799476563639605D0 c For R0=7/pi c c(1)=-0.8700992160947173D0 c For R0=6/pi c c(1)=-1.011848779711143D0 c For R0=5/pi c c(1)= -1.210283580515592D0 dcda(0)=0.D0 dcda(1)=0.D0 c do n=2,ncoef,1 c(n)=-1.D0/(n*(n-1))*( & eps*(n-1)*(n-3)*c(n-1)+ & (lamda)*(c(n-2)+ & 3.D0*eps*c(n-3)+ & 3.D0*eps**2*c(n-4)+ & eps**3*c(n-5))- & beta*(c(n-2)+eps*c(n-3))) enddo c y=0.D0 c dy=ncoef*c(ncoef)*x c d2y=ncoef*(ncoef-1)*c(ncoef)*x dy=0.D0 d2y=0.D0 c do n=ncoef,1,-1 y=(y+(c(n)))*x enddo do n=ncoef,2,-1 dy=(dy+(n*c(n)))*x enddo do n=ncoef,3,-1 d2y=(d2y+(n*(n-1)*c(n)))*x enddo y=y+c(0) dy=(dy+(c(1))) d2y=(d2y+(2*c(2))) return end c-----------------------------------------------------------------------