program main !*****************************************************************************80 ! !! MAIN is the main program for CYSTAL_QED. ! ! Discussion: ! ! CRYSTAL_QED seeks parameters that minimize a crystallization cost functional. ! ! Discussion: ! ! The DQED package is used to solve the bounded and constrained least ! squares problem. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! implicit none ! ! Parameters that don't depend on anything else. ! integer, parameter :: liopt = 17 integer, parameter :: lropt = 1 integer, parameter :: mcon = 6 integer, parameter :: mequa = 1 integer, parameter :: mvars = 7 integer, parameter :: nt = 5 ! ! First order parameters. ! integer, parameter :: ldfj = mcon+mequa integer, parameter :: nall = mcon+2*mvars+nt+1 ! ! Second order parameters. ! integer, parameter :: liwork = 3*mcon+9*mvars+4*nt+nall+11 integer, parameter :: lwork = nall*nall+4*nall+mvars*nt+33*mvars & +mequa*nt+mvars*nt+14*nt+9*mcon+26+3*nall+2 double precision bl(mvars+mcon) double precision bu(mvars+mcon) double precision dnrm2 external dqed_evaluate double precision fj(ldfj,mvars+1) double precision fnorm external func double precision fx(mcon+mequa) integer i integer igo integer ind(mvars+mcon) integer iopt(liopt) integer iopti integer ivars integer iwork(liwork) integer nvars double precision ropt(lropt) real tarray(2) double precision temp real time1 real time2 double precision value double precision work(lwork) double precision xpar(mvars) call timestamp ( ) ! ! Read the system CPU clock at the starting time. ! call cpu_time ( time1 ) ! ! Say hello. ! call hello ( liwork, lwork ) ! ! Tell DQED the size of the work arrays. ! iwork(1) = lwork iwork(2) = liwork ! ! Set the initial parameter values. ! ivars = 0 ! ! These are Hui Zhang's original crucible shape parameters. ! if ( ivars == 0 ) then nvars = 7 xpar(1) = 0.015D+00 xpar(2) = 0.03D+00 xpar(3) = 0.046D+00 xpar(4) = 0.07D+00 xpar(5) = 0.10D+00 xpar(6) = 0.14D+00 xpar(7) = 0.18D+00 ! ! These crucible shape parameters represent the optimum solution ! when the nodes are monotonically constrained. ! else if ( ivars == 1 ) then nvars = 7 xpar(1) = 0.064D+00 xpar(2) = 0.064D+00 xpar(3) = 0.064D+00 xpar(4) = 0.25D+00 xpar(5) = 0.25D+00 xpar(6) = 0.25D+00 xpar(7) = 0.25D+00 ! ! These crucible shape parameters cause the grid routine to fail, ! by generating degenerate control volumes. ! else if ( ivars == 2 ) then nvars = 7 value = 0.30D+00 xpar(1) = 3.0D+00*value/5.0D+00 xpar(2) = value xpar(3) = value-(value-0.25D+00)/5.0D+00 xpar(4) = value-2.0D+00*(value-0.25D+00)/5.0D+00 xpar(5) = value-3.0D+00*(value-0.25D+00)/5.0D+00 xpar(6) = value-4.0D+00*(value-0.25D+00)/5.0D+00 xpar(7) = value-4.5D+00*(value-0.25D+00)/5.0D+00 end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Initial parameter values:' write ( *, '(a)' ) ' ' do i = 1, nvars write ( *, '(i6,g14.6)' ) i, xpar(i) end do ! ! Set the constraints on the variables. ! My first guess is to just restrain all the variables to be ! between 0 and 0.25. ! if ( ivars == 0 .or. ivars == 1 .or. ivars == 2 ) then ind(1) = 3 bl(1) = 0.0D+00 bu(1) = 0.25D+00 ind(2) = 3 bl(2) = 0.0D+00 bu(2) = 0.25D+00 ind(3) = 3 bl(3) = 0.0D+00 bu(3) = 0.25D+00 ind(4) = 3 bl(4) = 0.0D+00 bu(4) = 0.25D+00 ind(5) = 3 bl(5) = 0.0D+00 bu(5) = 0.25D+00 ind(6) = 3 bl(6) = 0.0D+00 bu(6) = 0.25D+00 ind(7) = 3 bl(7) = 0.0D+00 bu(7) = 0.25D+00 ! ! The following 6 constraints force monotonicity. ! ind(8) = 1 bl(8) = 0.0D+00 bu(8) = 0.0D+00 ind(9) = 1 bl(9) = 0.0D+00 bu(9) = 0.0D+00 ind(10) = 1 bl(10) = 0.0D+00 bu(10) = 0.0D+00 ind(11) = 1 bl(11) = 0.0D+00 bu(11) = 0.0D+00 ind(12) = 1 bl(12) = 0.0D+00 bu(12) = 0.0D+00 ind(13) = 1 bl(13) = 0.0D+00 bu(13) = 0.0D+00 end if ! ! Set IOPTI, the switch which chooses ! just one function call (for checking), ! just one jacobian call (for checking), ! or optimization. ! iopti = 0 write ( *, '(a)' ) ' ' if ( iopti == 0 ) then write ( *, '(a)' ) 'Just call the function evaluator once.' else if ( iopti == 1 ) then write ( *, '(a)' ) 'Just call the jacobian approximator once.' else write ( *, '(a)' ) 'Optimize the problem.' end if ! ! If IOPTI = 0, just call the function evaluator once. ! if ( iopti == 0 ) then call func(fx,iopt,mcon,mequa,nvars,ropt,xpar) fnorm = dnrm2 ( mequa, fx(mcon+1), 1 ) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) 'Two-norm of residual is ', fnorm ! ! If IOPTI = 1, just call the jacobian approximator once. ! else if ( iopti == 1 ) then call diffor(fj,func,fx,iopt,ldfj,mcon,mequa,nvars,ropt,xpar) ! ! If IOPTI = 2, then carry out an optimization. ! else igo = 0 ! ! IOPT(1) = 2 means change the value of ITMAX to IOPT(2). ! iopt(1) = 2 iopt(2) = 10 write ( *, '(a)' ) ' ' write ( *, * ) 'Maximum number of minimization steps is ',iopt(2) ! ! IOPT(3) = 99 means no more options. ! iopt(3) = 99 call dqed ( dqed_evaluate, mequa,nvars,mcon,ind,bl,bu,xpar,fj,ldfj, & fnorm,igo,iopt,ropt,iwork,work) write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'DQED return flag IGO = ', igo write ( *, '(a,i6)' ) 'DQED real work array requirement = ', iwork(1) write ( *, '(a,i6)' ) 'DQED integer work array requirement = ', iwork(2) write ( *, '(a,g14.6)' ) 'Two-norm of residual is ', fnorm write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DQED returned the parameter values:' write ( *, '(a)' ) ' ' do i = 1,nvars write ( *, '(i6,g14.6)' ) i,xpar(i) end do write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) 'with the function value ', fnorm call func(fx,iopt,mcon,mequa,nvars,ropt,xpar) fnorm = dnrm2(mequa,fx(mcon+1),1) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I check the output function value!' write ( *, '(a,g14.6)' ) 'I get a function value ', fnorm end if ! ! Having completed the desired task, read the system clock and stop. ! call cpu_time ( time2 ) write ( *, '(a)' ) 'CRYSTAL_QED:' temp = time2-time1 write ( *, '(a,g14.6,a)' ) ' Total elapsed CPU time = ',temp,' seconds,' temp = temp/60.0 write ( *, '(a,g14.6,a)' ) ' = ',temp,' minutes.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CRYSTAL_QED:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine dqed_evaluate ( xpar, fj, ldfj, igo, iopt, ropt ) !*****************************************************************************80 ! !! DQED_EVALUATE evaluates certain functions being treated by DQED. ! ! Discussion: ! ! DQED_EVALUATE also evaluates partial derivatives. ! ! The user has NVARS variables XPAR(I), and is trying to minimize ! the square root of the sum of the squares of MEQUA functions ! F(I)(XPAR), subject to MCON constraints which have the form ! ! BL(I) < = G(I)(XPAR) <= BU(I) ! ! where either the left or right bounding inequality may be dropped. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, double precision XPAR(*). ! XPAR is an array of length NVARS, containing the values of the ! independent variables at which the functions and partial ! derivatives should be evaluated. ! ! Output, double precision FJ(LDFJ,NVARS+1). ! ! If IGO is nonzero, then partial derivatives must ! be placed in the first NVARS columns of FJ, as ! follows: ! ! Rows I = 1 to MCON, and columns J = 1 to NVARS ! should contain the partials of the I-th constraint ! function G(I) with respect to the J-th variable. ! ! Rows I = MCON+1 to MCON+MEQUA, and columns J = 1 to NVARS ! should contain the partials of the (I-MCON)-th nonlinear ! function F(I-MCON) with respect to the J-th variable. ! ! Regardless of the value of IGO, column NVARS+1 must be ! assigned the values of the constraints and functions, ! as follows: ! ! Rows I = 1 to MCON, column J = NVARS+1, should contain ! the value of G(I). ! ! Rows I = MCON+1 to MCON+MEQUA, column J = NVARS+1, should ! contain the value of F(I-MCON). ! ! Input, integer LDFJ. ! LDFJ is the leading dimension of FJ, which must be at least ! MCON+MEQUA. ! ! Input/output, integer IGO. ! On input, IGO tells the user whether the partial derivatives ! are needed. ! ! 0, the partial derivatives are not needed. ! nonzero, the partial derivatives are needed. ! ! On output, the user may reset the input value of IGO if one ! of two situations is encountered: ! ! 99, the functions, constraints, or partial derivatives ! could not be evaluated at the given input point X. Request ! that DQED reject that point, and try a different one. ! ! Any other value, abort the run. ! ! Input, integer IOPT(*), the integer option array. ! ! Input, double precision ROPT(*), the real option array. ! implicit none ! integer ldfj integer, parameter :: mcon = 6 integer, parameter :: mequa = 1 integer, parameter :: nvars = 7 ! double precision fj(ldfj,nvars+1) external func double precision fx(mcon+mequa) integer igo integer iopt(*) double precision ropt(*) double precision xpar(nvars) ! if ( igo /= 0 ) then call diffor(fj,func,fx,iopt,ldfj,mcon,mequa,nvars,ropt,xpar) end if call func(fj(1,nvars+1),iopt,mcon,mequa,nvars,ropt,xpar) return end subroutine dqedhd ( x, fj, ldfj, igo, iopt, ropt ) !*****************************************************************************80 ! !! DQEDHD evaluates functions and derivatives for DQED. ! ! Discussion: ! ! For our purposes, this is a dummy routine, and is not called. ! The compiler is happier when it is here, though. ! ! The user problem has MCON constraint functions, ! MEQUA least squares equations, and involves NVARS ! unknown variables. ! ! When this subprogram is entered, the general (near) ! linear constraint partial derivatives, the derivatives ! for the least squares equations, and the associated ! function values are placed into the array FJ. ! all partials and functions are evaluated at the point ! in X. then the subprogram returns to the calling ! program unit. typically one could do the following ! steps: ! ! if ( igo /= 0 ) then ! place the partials of the i-th constraint function with respect to ! variable j in the array fj(i,j), i = 1,...,mcon, j=1,...,nvars. ! ! place the values of the i-th constraint equation into fj(i,nvars+1). ! ! if ( igo /= 0 ) then ! place the partials of the i-th least squares equation with respect ! to variable j in the array fj(i,j), i = 1,...,mequa, j = 1,...,nvars. ! ! place the value of the i-th least squares equation into fj(i,nvars+1). ! implicit none ! integer ldfj integer, parameter :: nvars = 8 ! double precision fj(ldfj,nvars+1) integer igo integer iopt(*) integer mcon integer mequa integer mode double precision ropt(*) double precision x(nvars) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DQEDHD - Fatal error!' write ( *, '(a)' ) ' This is a dummy routine.' stop end subroutine func ( fx, iopt, mcon, mequa, nvars, ropt, xpar ) !*****************************************************************************80 ! !! FUNC communicates between the optimizer and the constitutive routines. ! ! Discussion: ! ! FUNC receives a set of parameter values as input, solves the ! resulting constitutive equations, and evaluates the derived ! cost functional, whose value is returned to the calling routine. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IOPT(*), the integer option array. ! ! Input, integer MCON. ! MCON is the number of constraints imposed on the independent ! variables that compose a feasible solution. ! ! Input, integer MEQUA. ! MEQUA is the number of components of the nonlinear function. ! ! Input, integer NVARS. ! NVARS is the number of parameters or independent variables ! upon which the function depends. ! ! Input, double precision ROPT(*). ! ROPT is the real option array. ! ! Input, double precision XPAR(*). ! XPAR is an array of length NVARS, containing the values of the ! parameters. ! implicit none ! integer, parameter :: maxbot = 20 integer mcon integer mequa integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nk = 14 integer, parameter :: ns = 10 integer nvars ! double precision ae1(ni,nj) double precision ae2(ni,nj) double precision ak1(ni,nj) double precision ak2(ni,nj) double precision area(ni,nj) double precision areal double precision areas double precision areat double precision b double precision b1jbl(nj) double precision b2jbl(nj) double precision b3jbl(ni) double precision birad double precision bo double precision cappa double precision cd double precision ce1 double precision ce2 double precision cfo double precision cinc double precision cinco double precision cmu double precision cost double precision cvn double precision damax double precision delt double precision dtm double precision epsad double precision epsil double precision ewall double precision f(ni,nj,ns) double precision fcsl double precision fjeta(ni,nj) double precision fjksi(ni,nj) double precision fks double precision fksl double precision fma double precision fmax(ns) double precision fn(ni,nj,ns) double precision fnu double precision fo(ni,nj,ns) double precision fr double precision frsl double precision fx(mcon+mequa) double precision gam(ni,nj) double precision gamt(ni,nj) double precision grash double precision hamag double precision heta(ni,nj) double precision hf double precision hksi(ni,nj) integer i integer icost integer icrys integer inturb integer iopt(1) integer iplot integer ipref integer iprint integer iter integer izone integer j integer jcrys integer jpref integer k integer l0 integer l1 integer last integer lastt logical lblk(ns) logical lconv logical lortho logical lsolve(ns) integer m0 integer m1 integer mode integer nbot integer, save :: ncall = 0 integer ndt integer nf integer np integer npc integer nsolve(ns) integer ntimes(ns) double precision orth double precision pr double precision r(ni,nj) double precision ra double precision rdtm double precision re double precision recb double precision rect double precision relax(nk) double precision res(ns) double precision rho(ni,nj) double precision rhocon double precision ropt(*) double precision rpr double precision rueta(ni,nj) double precision ruksi(ni,nj) double precision sige double precision sigk double precision sigma double precision sigt double precision smax double precision smooth double precision ssum double precision stel double precision stes double precision t1 double precision t2 double precision t3 double precision tal double precision tanca double precision tanca2 double precision tas double precision tend double precision tf double precision tinit character ( len = 25 ) title(ns) double precision tnow double precision tw double precision vave double precision vol(ni,nj) double precision vort(ni,nj) double precision x(ni,nj) double precision xbot(maxbot) double precision xc(ni,nj) double precision xlen double precision xpar(nvars) double precision y(ni,nj) double precision ybot(maxbot) double precision yc(ni,nj) double precision ylen ! ncall = ncall+1 ! ! Initialize the data. ! call inidat(ae1,ae2,ak1,ak2,area,b,b1jbl,b2jbl,b3jbl, & birad,bo,cappa,cd,ce1, & ce2,cfo,cinc,cmu,cost,cvn,delt,dtm,epsad,epsil,ewall,f,fcsl, & fks,fksl,fma,fn,fnu,fo,fr,frsl,gam,gamt,grash,hamag,heta,hf, & hksi,icost,icrys,inturb,iplot,ipref,iprint,jcrys, & jpref,l0,l1,last,lastt,lblk,lortho,lsolve,m0,m1,maxbot,mode, & nbot,ndt,ni,nj,nk,np,npc,ns,nsolve,ntimes,orth, & pr,ra,rdtm, & re,recb,rect,relax,rho,rhocon,rpr,rueta,ruksi,sige,sigk, & sigma,sigt,smooth,stel,stes,tal,tanca,tanca2,tas,tend,tf, & tinit,title,tnow,tw,vave,vol,vort,x,xbot,xc,xlen,y,ybot,yc,ylen) ! ! TEMPORARY ! do i = 1,7 xbot(i+1) = xpar(i) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNC:' write ( *, '(a)' ) ' Input parameters:' do i = 1,7 write ( *, '(i6,g14.6)' ) i,xpar(i) end do ! ! Print out data. ! if ( ncall == 1 ) then call prdat(b,birad,bo,cappa,cfo,cvn,delt,dtm,fcsl,fksl,fma, & fnu,fr,frsl,grash,icost,icrys,inturb,iprint, & jcrys,l0,last,lastt,m0,mode,nbot,ns,nsolve,ntimes,orth, & pr,ra,rdtm,recb,rect,rhocon,smooth,stel,stes,tanca,tanca2, & tend,tf,tinit,title,tw,vave,xlen,ylen) end if ! ! Generate the initial grid XC, YC. ! call inigrd(icrys,jcrys,l0,m0,nbot,xbot,xc,xlen,ybot,yc,ylen) ! ! Adapt the grid. ! do k = 1, 10 call adapt(cvn,epsad,icrys,iprint,jcrys,l1,m1,nbot, & orth,smooth,xbot,xc,ybot,yc) end do ! ! Once XC and YC are determined, compute the control volume areas. ! call doarea(area,areal,areas,areat,icrys,jcrys,l0,m0,ni,nj,xc,yc) ! ! Each time iteration begins at this point. ! 10 continue ! ! Begin the iterative solution of the state equations in the ! solid zone. ! izone = 1 l1 = l0 m1 = jcrys ! ! Only the temperature (variable 5) needs to be solved for. ! lsolve(1) = .false. lsolve(2) = .false. lsolve(3) = .false. lsolve(4) = .false. lsolve(5) = .true. lsolve(6) = .false. lsolve(7) = .false. lsolve(8) = .false. lsolve(9) = .false. lsolve(10) = .false. do iter = 0,last lconv = .true. ! ! Set the X and Y coordinates of the primary nodes from XC, YC, ! the coordinates of the corner nodes. ! ! Note that this is only done for the left portion of the region! ! call setx ( l1, m1, ni, nj, x, xc, y, yc ) ! ! Compute various geometric quantities required for computing ! derivatives. ! call setgeo(ae1,ae2,ak1,ak2,heta,hksi,l1,m1,mode,ni,nj, & r,vol,x,xc,y,yc) ! ! Estimate pressure and momentum at control volume interfaces. ! if ( ndt == 0 ) then if ( iter == 0 ) then call setp(fjeta,fjksi,heta,hksi,l1,m1,ni,nj,f(1,1,3)) call setru(heta,hksi,l1,m1,ni,nj,rho,rueta,ruksi, & f(1,1,1),f(1,1,2),x,y) end if end if call setup(ae1,ae2,ak1,ak2,b1jbl,b2jbl,b3jbl,birad, & cappa,cd,ce1,ce2,cmu,epsil,ewall,f, & fcsl,fjeta,fjksi,fksl,fma,fmax,fn,fo,frsl,gamt,grash,hamag, & heta,hksi,inturb,icrys,iter,izone,jcrys,l1,lblk,lconv,lortho, & lsolve,m1,mode,nf,np,npc,nsolve,ntimes,pr,r,rdtm,re,recb, & rect,res,relax,rho,rhocon,rpr,rueta,ruksi,sige,sigk,sigt, & smax,ssum,stel,stes,tal,tas,tf,tw,vol,x,xc,y,yc) if ( iprint > 0 ) then call output(cfo,iter,izone,res,smax,ssum,f(1,1,5), & tnow,f(1,1,1),f(1,1,6)) end if if ( lconv .and. iter >= 5)go to 20 end do if ( iprint > 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNC - Warning!' write ( *, '(a)' ) ' The solid iteration has not converged' write ( *, '(a,i6,a)' ) ' after ',iter,' iterations.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Solid norm and max relative change:' write ( *, '(a)' ) ' ' do i = 1,ns if ( lsolve(i) ) then fmax(i) = damax ( l0*m0, f(1,1,i), 1 ) write(*,'(i3,2x,a25,2g14.6)')i,title(i),fmax(i),res(i) end if end do end if ! ! Zone 2 (LIQUID) calculation ! 20 continue do i = 1,icrys do j = 1,jcrys f(i,j,5) = fo(i,j,5) f(i,j,9) = fo(i,j,9) end do end do izone = 2 l1 = icrys m1 = m0 do i = 1,icrys ruksi(i,1) = 0.0D+00 ruksi(i,m0) = 0.0D+00 rueta(i,2) = 0.0D+00 rueta(i,m0) = 0.0D+00 end do do j = 1,m0 rueta(1,j) = 0.0D+00 rueta(icrys,j) = 0.0D+00 ruksi(2,j) = 0.0D+00 ruksi(icrys,j) = 0.0D+00 end do if ( inturb == 0 ) then lsolve(1) = .true. lsolve(2) = .true. lsolve(3) = .true. lsolve(4) = .true. lsolve(5) = .true. lsolve(6) = .false. lsolve(7) = .false. lsolve(8) = .false. lsolve(9) = .false. lsolve(10) = .false. else lsolve(1) = .true. lsolve(2) = .true. lsolve(3) = .true. lsolve(4) = .true. lsolve(5) = .true. lsolve(6) = .false. lsolve(7) = .true. lsolve(8) = .true. lsolve(9) = .false. lsolve(10) = .false. end if do iter = 1,last lconv = .true. ! ! Set the turbulent viscosity. ! if ( inturb == 0 ) then do i = 2,l1-1 do j = 2,m1-1 gamt(i,j) = 0.0D+00 end do end do else if ( inturb == 1 ) then do i = 2,l1-1 do j = 2,m1-1 gamt(i,j) = (1.0D+00-relax(12))*gamt(i,j) & +relax(12)*cmu*rho(i,j)*f(i,j,7)**2/f(i,j,8) end do end do end if ! ! Set the X, Y values. ! call setx(l1,m1,ni,nj,x,xc,y,yc) ! ! Compute various geometric quantities required for computing ! derivatives. ! call setgeo(ae1,ae2,ak1,ak2,heta,hksi,l1,m1,mode,ni,nj, & r,vol,x,xc,y,yc) ! ! Set the right hand side of certain flux boundary conditions. ! call setup(ae1,ae2,ak1,ak2,b1jbl,b2jbl,b3jbl,birad, & cappa,cd,ce1,ce2,cmu,epsil,ewall,f, & fcsl,fjeta,fjksi,fksl,fma,fmax,fn,fo,frsl,gamt,grash,hamag, & heta,hksi,inturb,icrys,iter,izone,jcrys,l1,lblk,lconv, & lortho, & lsolve,m1,mode,nf,np,npc,nsolve,ntimes,pr,r,rdtm,re,recb, & rect,res,relax,rho,rhocon,rpr,rueta,ruksi,sige,sigk,sigt, & smax,ssum,stel,stes,tal,tas,tf,tw,vol,x,xc,y,yc) ! ! Compute the stream function PSI. ! f(2,2,10) = 0.0D+00 do i = 2,icrys if ( i > 2 ) then f(i,2,10) = f(i-1,2,10)-rueta(i-1,2)*ae1(i-1,2) & -0.5D+00*(ruksi(i-1,1)+ruksi(i,1))*ae2(i-1,2) end if do j = 3,m1-1 t1 = (rueta(i,j-1)+rueta(i,j)) t2 = (rueta(i-1,j-1)+rueta(i-1,j)) t3 = 0.5D+00*(hksi(i,j-1)*t2+hksi(i-1,j-1)*t1)/(hksi(i,j-1)+hksi(i-1,j-1)) f(i,j,10) = f(i,j-1,10)+ruksi(i,j-1)*ak1(i,j-1)-t3*ak2(i,j-1) end do end do f(1,1,10) = f(2,2,10) do j = 2,m0 f(1,j,10) = f(2,j,10) end do do i = 2,l0 f(i,1,10) = f(i,2,10) end do ! ! Optional printout. ! if ( iprint > 0 ) then call output(cfo,iter,izone,res,smax,ssum,f(1,1,5),tnow,f(1,1,1),f(1,1,6)) end if if ( lconv)go to 30 end do if ( iprint > 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNC - Warning!' write ( *, '(a)' ) ' The liquid iteration has not converged' write ( *, '(a,i6,a)' ) ' after ',iter,' iterations.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Liquid norms and max relative change:' write ( *, '(a)' ) ' ' do i = 1,ns if ( lsolve(i) ) then fmax(i) = damax(l0*m0,f(1,1,i),1) write(*,'(i3,2x,a25,2g14.6)')i,title(i),fmax(i),res(i) end if end do end if ! ! We have computed the state variables for the current time. ! 30 continue ! ! Compute the vorticity VORT. ! call vortic ( heta, hksi, icrys, l0, m0, f(1,1,1), f(1,1,2), vort, xc, yc ) ! ! Evaluate the cost function integrand at the current time. ! cinco = cinc call setcst(area,cinc,gam,icost,icrys,jcrys,l0,m0,ni,nj,f(1,1,5), & tf,f(1,1,1),f(1,1,2),vave,vort,xc,yc) ! ! Update the total cost, by adding the estimated contribution ! from the current time interval. ! ! TEMPORARY ! if ( ndt > 0 ) then ! cost = cost+0.5D+00*delt*(cinco+cinc) cost = cost+0.5D+00*dtm*(cinco+cinc) end if ! ! If we've reached the end time, then write out final data, ! possibly save a restart file, and stop. ! if ( ndt >= lastt ) then if ( iprint > 0 ) then call pmod(ipref,jcrys,jpref,l1,m1,f(1,1,3)) end if cost = cost/(sqrt(areal)*(tend-tinit)) ! ! If requested, write out plot data. ! Sadly, it is necessary to call SETX to generate the X, Y arrays ! for the entire region. Previous calls only generate them for ! a portion of the region. ! if ( iplot == 1 ) then call setx(l0,m0,ni,nj,x,xc,y,yc) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'HEY, ABOUT TO CALL RSWRIT.' call rswrit(b1jbl,b2jbl,b3jbl,cost,f(1,1,9),gamt, & icrys,jcrys,l0,m0,nbot,f(1,1,3),f(1,1,4),f(1,1,10), & rueta,ruksi,f(1,1,5),f(1,1,8),f(1,1,7),tnow,f(1,1,1), & f(1,1,2),vort,f(1,1,6),x,xbot,xc,y,ybot,yc) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'HEY, ABOUT TO CALL WRTEC.' call wrtec(gamt,l1,m1,f(1,1,3),f(1,1,10),f(1,1,5),f(1,1,8), & f(1,1,7),f(1,1,1),f(1,1,2),f(1,1,6),x,y) end if write ( *, '(a,g14.6)' ) 'FUNC: COST = ',cost do i = 1,mcon fx(i) = xpar(i+1)-xpar(i) end do fx(mcon+1) = cost return end if ! ! Update the number of steps, and the current time. ! ndt = ndt+1 tnow = tinit+ndt*dtm ! ! Save a copy of the current data. ! do i = 1,l0 do j = 1,m0 do k = 1,ns fo(i,j,k) = f(i,j,k) end do end do end do call movgrd(b1jbl,b2jbl,bo,delt,fksl,fr,frsl,icrys, & iprint,jcrys,l0,m0,f(1,1,3),pr,re,stel,tanca,tanca2,xc, & xlen,yc) l1 = l0 m1 = m0 call adapt(cvn,epsad,icrys,iprint,jcrys,l1,m1,nbot, & orth,smooth,xbot,xc,ybot,yc) ! ! Once XC and YC are determined, compute the control volume areas. ! call doarea(area,areal,areas,areat,icrys,jcrys,l0,m0,ni,nj,xc,yc) go to 10 end subroutine adapt ( cvn, epsad, icrys, iprint, jcrys, l1, m1, nbot, & orth, smooth, xbot, xc, ybot, yc ) !*****************************************************************************80 ! !! ADAPT executes MAGG, the multizone adaptive grid generation. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none integer nbot integer, parameter :: ni = 64 integer, parameter :: nj = 64 double precision a1 double precision a11 double precision a12 double precision a2 double precision a21 double precision a22 double precision a3 double precision b1 double precision b2 double precision b3 double precision c1 double precision c2 double precision c3 double precision cvn double precision det double precision dot double precision dx double precision dxc(ni,nj) double precision dxmax double precision dy double precision dyc(ni,nj) double precision epsad double precision gn double precision gt double precision gx double precision gy integer i integer icrys integer iprint integer j integer jcrys integer l1 integer m1 integer nin double precision orth double precision pb(ni,nj) double precision pc1 double precision pc2 double precision pd(ni,nj) double precision ratio double precision res1 double precision res2 double precision rmax double precision s1 double precision s2 double precision scale double precision smooth double precision ssmot double precision ssweg double precision vmag2 double precision wn double precision wt double precision wx double precision wy double precision xbot(nbot) double precision xc(ni,nj) double precision xdisp double precision xij double precision xjac double precision xn double precision xnn double precision xpp double precision xt double precision xtn double precision xtt double precision ybot(nbot) double precision yc(ni,nj) double precision ydisp double precision yij double precision yn double precision ynn double precision ypp double precision yt double precision ytn double precision ytt ! dxmax = 0.0D+00 do i = 1,l1 do j = 1,m1 dxc(i,j) = 0.0D+00 dyc(i,j) = 0.0D+00 end do end do do i = 2,l1 do j = 2,m1 if ( i < icrys ) then pc1 = 0.95D+00*((i-(2+icrys)/2)/10.0D+00)**2+0.05D+00 else pc1 = 0.8D+00*((i-l1)/12.0)**2+0.2D+00 end if if ( j <= jcrys ) then pc2 = 0.9D+00*((j-(2+jcrys)/2)/10.0D+00)**2+0.05D+00 else pc2 = 0.9D+00*((j-(m1+jcrys)/2)/10.0D+00)**2+0.05D+00 end if pb(i,j) = pc1*pc2 pd(i,j) = pc1*pc2 end do end do ! ! This DO loop iteration is an iterative solution of a linear system. ! do nin = 1,100 ssmot = 0.0D+00 ssweg = 0.0D+00 rmax = 0.0D+00 if ( iprint > 0 ) then if ( nin == 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ADAPT: Iter DXMAX XC(20,20) YC(20,20)' // & ' XC(3,M1) YC(3,M1) RMAX' end if if ( mod(nin,50) == 0 ) then write(*,'(6x,i5,6g10.3)') & nin,dxmax,xc(20,20),yc(20,20),xc(3,m1),yc(3,m1),rmax end if end if do j = 3,m1-1 do i = 3,l1-1 xt = 0.5D+00*(xc(i+1,j)-xc(i-1,j)) xn = 0.5D+00*(xc(i,j+1)-xc(i,j-1)) yt = 0.5D+00*(yc(i+1,j)-yc(i-1,j)) yn = 0.5D+00*(yc(i,j+1)-yc(i,j-1)) xjac = xt*yn-xn*yt ! ! Refer to table 2.1, page 14, of Hui Zhang's thesis, for ! a description of these coefficients. ! a1 = -2.0D+00*smooth*(xt*yt+xn*yn)*(xn*xn+yn*yn)/xjac**3 a2 = 4.0D+00*smooth*(xt*yt+xn*yn)*(xt*xn+yt*yn)/xjac**3 a3 = -2.0D+00*smooth*(xt*yt+xn*yn)*(xt*xt+yt*yt)/xjac**3 b1 = 2.0D+00*smooth*(yt*yt+yn*yn)*(xn*xn+yn*yn)/xjac**3 b2 = -4.0D+00*smooth*(yt*yt+yn*yn)*(xt*xn+yt*yn)/xjac**3 b3 = 2.0D+00*smooth*(yt*yt+yn*yn)*(xt*xt+yt*yt)/xjac**3 c1 = 2.0D+00*smooth*(xt*xt+xn*xn)*(xn*xn+yn*yn)/xjac**3 c2 = -4.0D+00*smooth*(xt*xt+xn*xn)*(xt*xn+yt*yn)/xjac**3 c3 = 2.0D+00*smooth*(xt*xt+xn*xn)*(xt*xt+yt*yt)/xjac**3 a1 = a1+orth*pd(i,j)*xn*yn a2 = a2+orth*pd(i,j)*(xt*yn+xn*yt) a3 = a3+orth*pd(i,j)*xt*yt b1 = b1+orth*pd(i,j)*xn*xn b2 = b2+2.0D+00*orth*pd(i,j)*(2.0D+00*xt*xn+yt*yn) b3 = b3+orth*pd(i,j)*xt*xt c1 = c1+orth*pd(i,j)*yn*yn c2 = c2+2.0D+00*orth*pd(i,j)*(xt*xn+2.0D+00*yt*yn) c3 = c3+orth*pd(i,j)*yt*yt a1 = a1-2.0D+00*cvn*pb(i,j)*xn*yn a2 = a2+2.0D+00*cvn*pb(i,j)*(xt*yn+xn*yt) a3 = a3-2.0D+00*cvn*pb(i,j)*xt*yt b1 = b1+2.0D+00*cvn*pb(i,j)*yn*yn b2 = b2-4.0D+00*cvn*pb(i,j)*yt*yn b3 = b3+2.0D+00*cvn*pb(i,j)*yt*yt c1 = c1+2.0D+00*cvn*pb(i,j)*xn*xn c2 = c2-4.0D+00*cvn*pb(i,j)*xt*xn c3 = c3+2.0D+00*cvn*pb(i,j)*xt*xt ! ! Now compute terms from the derivatives of the weight functions. ! gt = 0.5D+00*(pd(i+1,j)-pd(i-1,j)) gn = 0.5D+00*(pd(i,j+1)-pd(i,j-1)) gx = orth*0.5D+00*(gt*yn-gn*yt)*(xt*xn+yt*yn)**2/xjac gy = orth*0.5D+00*(gn*xt-gt*xn)*(xt*xn+yt*yn)**2/xjac wt = 0.5D+00*(pb(i+1,j)-pb(i-1,j)) wn = 0.5D+00*(pb(i,j+1)-pb(i,j-1)) wx = cvn*xjac*(yn*wt-yt*wn) wy = cvn*xjac*(xt*wn-xn*wt) ! ! (RES1, RES2) is the residual that should be driven to zero. ! xtt = xc(i+1,j)-2.0D+00*xc(i,j)+xc(i-1,j) xnn = xc(i,j+1)-2.0D+00*xc(i,j)+xc(i,j-1) ytt = yc(i+1,j)-2.0D+00*yc(i,j)+yc(i-1,j) ynn = yc(i,j+1)-2.0D+00*yc(i,j)+yc(i,j-1) xtn = 0.25D+00*(xc(i+1,j+1)+xc(i-1,j-1)-xc(i-1,j+1)-xc(i+1,j-1)) ytn = 0.25D+00*(yc(i+1,j+1)+yc(i-1,j-1)-yc(i-1,j+1)-yc(i+1,j-1)) res1 = b1*xtt+b2*xtn+b3*xnn+a1*ytt+a2*ytn+a3*ynn+gx+wx res2 = a1*xtt+a2*xtn+a3*xnn+c1*ytt+c2*ytn+c3*ynn+gy+wy if ( abs(res1)+abs(res2) > rmax ) then rmax = abs(res1)+abs(res2) end if a11 = -2.0D+00*(b1+b3) a12 = -2.0D+00*(a1+a3) a22 = -2.0D+00*(c1+c3) a21 = a12 det = a11*a22-a12*a21 ! ! (DX, DY) is the pointwise solution of the linear system. ! dx = (res2*a21-res1*a22)/det dy = (res1*a12-res2*a11)/det ! ! Now look at how movement of XC(I,J), YC(I,J) affects the neighbors ! ! (I+1,J-1) (I+1,J) (I+1,J+1) ! ! (I, J-1) (I, J) (I, J+1) ! ! (I-1,J-1) (I-1,J) (I-1,J+1) ! ! to determine the linear factor SCALE that will multiply (DX,DY). ! ratio = 0.25D+00 xij = xc(i,j) yij = yc(i,j) xdisp = xc(i+1,j)-xij ydisp = yc(i+1,j)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i+1,j+1)-xij ydisp = yc(i+1,j+1)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i,j+1)-xij ydisp = yc(i,j+1)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i-1,j+1)-xij ydisp = yc(i-1,j+1)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i-1,j)-xij ydisp = yc(i-1,j)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i-1,j-1)-xij ydisp = yc(i-1,j-1)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i,j-1)-xij ydisp = yc(i,j-1)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i+1,j-1)-xij ydisp = yc(i+1,j-1)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) scale = min(0.5D+00,0.25D+00/ratio) dxc(i,j) = scale*dx dyc(i,j) = scale*dy ssmot = ssmot+(xt*xt+xn*xn+yt*yt+yn*yn)/xjac ssweg = ssweg+cvn*pb(i,j)*xjac**2 end do end do ! ! Move the points. ! do j = 3,m1-1 do i = 3,l1-1 if ( i /= icrys ) then xc(i,j) = xc(i,j)+dxc(i,j) end if if ( j /= jcrys ) then yc(i,j) = yc(i,j)+dyc(i,j) end if end do end do do i = icrys+1,l1-1 yc(i,jcrys) = 0.4D+00+0.02D+00*sin(xc(i,jcrys)*20.0D+00*3.14159D+00) end do ! ! Handle the nodes along the bottom row, from 10 positions to the right ! of the crystal, to the right hand wall. ! do j = jcrys+10,m1-1 ratio = 0.25D+00 call findp(j,xc(3,j),yc(3,j),xpp,ypp,xc,yc) dx = xpp-xc(2,j) dy = ypp-yc(2,j) xdisp = xc(2,j+1)-xc(2,j) ydisp = yc(2,j+1)-yc(2,j) vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) s1 = vmag2 xdisp = xc(2,j-1)-xc(2,j) ydisp = yc(2,j-1)-yc(2,j) vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) s2 = max(vmag2,s1) scale = min(0.3D+00,0.25D+00/ratio) s1 = (scale*dx)**2+(scale*dy)**2 if ( s1 < s2 ) then xc(2,j) = xc(2,j)+scale*dx yc(2,j) = yc(2,j)+scale*dy end if call cubic(nbot,yc(2,j),ybot,xc(2,j),xbot) yc(l1,j) = yc(l1-1,j) end do ! ! Handle the nodes along the bottom row, from the left axis ! of symmetry, to 9 positions beyond the crystal. ! do j = 2,jcrys+9 yc(2,j) = yc(3,j) call cubic(nbot,yc(2,j),ybot,xc(2,j),xbot) yc(l1,j) = yc(l1-1,j) end do do i = 3,l1 xc(i,2) = xc(i,3) xc(i,m1) = xc(i,m1-1) if ( xc(i,m1-1) <= (xc(2,m1)+0.002D+00*dble(i-2)) ) then xc(i,m1) = xc(2,m1)+0.002D+00*dble(i-2) end if end do ! ! Compute the maximum movement of all corner nodes. ! dxmax = abs(dxc(3,3)) do j = 3,m1-1 do i = 3,l1-1 dxmax = max(dxmax,abs(dxc(i,j))) dxmax = max(dxmax,abs(dyc(i,j))) end do end do ! ! Copy values to the dummy corner nodes with I = 1 or J=1. ! xc(1,1) = xc(2,2) yc(1,1) = yc(2,2) do j = 2,m1 xc(1,j) = xc(2,j) yc(1,j) = yc(2,j) end do do i = 2,l1 xc(i,1) = xc(i,2) yc(i,1) = yc(i,2) end do ! ! Check for convergence. ! if ( dxmax <= epsad)return end do return end subroutine cubic ( nbot, ynew, ybot, xnew, xbot ) ! !*****************************************************************************80 ! !! CUBIC constructs a cubic spline through data. ! ! ! Discussion: ! ! CUBIC is given NBOT points (YBOT,XBOT) which lie along the curve ! defining the bottom of the crucible. ! ! CUBIC is also given a value YNEW which lies between YBOT(1) and ! YBOT(NBOT). ! ! CUBIC constructs a cubic spline through the data, and evaluates it at ! YNEW, returning the value as XNEW. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NBOT, the number of points to be used to ! define the bottom of the crucible. ! implicit none ! integer nbot integer, parameter :: ni = 64 ! double precision dome integer i integer ist double precision p2 double precision p3 double precision shift double precision t1 double precision t2 double precision tm(ni) double precision tt(ni) double precision x1 double precision x2 double precision xbot(nbot) double precision xl(ni) double precision xnew double precision y1 double precision y2 double precision ybot(nbot) double precision yl(ni) double precision ymin double precision ynew ! if ( ynew < ybot(1) .or. ybot(nbot).lt.ynew ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CUBIC - Fatal error!' write ( *, '(a,g14.6)' ) ' YNEW = ',ynew write ( *, '(a,g14.6)' ) ' outside of range YBOT(1) = ',ybot(1) write ( *, '(a,g14.6)' ) ' through YBOT(NBOT) = ',ybot(nbot) stop end if do i = 1,ni xl(i) = 0.0D+00 yl(i) = 0.0D+00 tt(i) = 0.0D+00 tm(i) = 0.0D+00 end do yl(1) = ybot(1)+(ybot(1)-ybot(3)) yl(2) = ybot(2)+(ybot(1)-ybot(3)) do i = 3,nbot+2 yl(i) = ybot(i-2) xl(i) = xbot(i-2) end do yl(nbot+3) = yl(nbot+1)+(yl(nbot+2)-yl(nbot)) yl(nbot+4) = yl(nbot+2)+(yl(nbot+2)-yl(nbot)) shift = (xl(3)-xl(4))/(yl(3)-yl(4))-(xl(4)-xl(5))/(yl(4)-yl(5)) xl(2) = xl(3)+(yl(2)-yl(3))*(shift+(xl(3)-xl(4))/(yl(3)-yl(4))) xl(1) = xl(2)+(yl(1)-yl(2))*(shift+(xl(2)-xl(3))/(yl(2)-yl(3))) shift = (xl(nbot+2)-xl(nbot+1))/(yl(nbot+2)-yl(nbot+1)) & -(xl(nbot+1)-xl(nbot))/(yl(nbot+1)-yl(nbot)) xl(nbot+3) = xl(nbot+2)+(yl(nbot+3)-yl(nbot+2))* & (shift+(xl(nbot+2)-xl(nbot+1))/(yl(nbot+2)-yl(nbot+1))) xl(nbot+4) = xl(nbot+3)+(yl(nbot+4)-yl(nbot+3))* & (shift+(xl(nbot+3)-xl(nbot+2))/(yl(nbot+3)-yl(nbot+2))) do i = 1,nbot+3 tm(i) = (xl(i+1)-xl(i))/(yl(i+1)-yl(i)) end do do i = 3,nbot+2 dome = abs(tm(i+1)-tm(i))+abs(tm(i-1)+tm(i-2)) if ( dome < 1.0D-10) then tt(i) = 0.0D+00 else tt(i) = (abs(tm(i+1)-tm(i))*tm(i-1)+abs(tm(i-1)-tm(i-2))*tm(i))/dome end if end do ! ! Find the node YL(IST) which is nearest to YNEW. ! ymin = abs(yl(3)-ynew) ist = 3 do i = 3,nbot+2 if ( abs(yl(i)-ynew) < ymin ) then ist = i ymin = abs(yl(i)-ynew) end if end do if ( (yl(ist)-ynew) > 0.0D+00 ) then y1 = yl(ist-1) x1 = xl(ist-1) y2 = yl(ist) x2 = xl(ist) t1 = tt(ist-1) t2 = tt(ist) else y1 = yl(ist) x1 = xl(ist) y2 = yl(ist+1) x2 = xl(ist+1) t1 = tt(ist) t2 = tt(ist+1) end if ! ! Evaluate the spline at YNEW. ! p2 = (3.0D+00*(x2-x1)/(y2-y1)-2.0D+00*t1-t2)/(y2-y1) p3 = (t1+t2-2.0*(x2-x1)/(y2-y1))/(y2-y1)**2 xnew = x1+t1*(ynew-y1)+p2*(ynew-y1)**2+p3*(ynew-y1)**3 return end subroutine diflow ( acof, diff, flow ) ! !*****************************************************************************80 ! !! DIFLOW computes the convection-diffusion coefficient ACOF. ! ! ! Discussion: ! ! DIFLOW is given ! FLOW, the mass velocity RHO*U, and DIFF, the value of GAMMA/DELX. ! Several schemes are available, but currently the power law is used. ! ! See Patankar, chapter 5. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! double precision acof double precision diff double precision flow ! acof = diff if ( diff == 0.0D+00) return ! ! Power Law. ! ! Define ! FE = rho*u ! DE = gamma/delx ! Then the coefficient is ! ! AE = DE * Max(0, (1-0.1*Abs(FE)/DE)**5) + Max(0,-FE) ! acof = diff*max(1.0D-10,(1.0D+00-0.1D+00*abs(flow/diff))**5) if ( flow < 0.0D+00) acof = acof-flow return end subroutine doarea(area,areal,areas,areat,icrys,jcrys,l0,m0,ni,nj,xc,yc) ! !*****************************************************************************80 ! !! DOAREA calculates control volume areas. ! ! ! Discussion: ! ! DOAREA is given (XC,YC), the locations of the "corners" of the ! control volumes, and calculates the area of each control volume. ! ! DOAREA uses the fact that the area of a polygon which is bounded ! by the N points (X(I),Y(I)) can be computed by: ! ! AREA = 0.5 * SUM (I=1 to N) X(I) * (Y(I+1)-Y(I-1)) ! ! where Y(0) should be replaced by Y(N), and Y(N+1) by Y(1). ! ! Using this formula, AREA may come out negative, depending on ! whether the nodes are given in clockwise or counterclockwise order. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, double precision AREA(NI,NJ). ! AREA contains the area of the control volume (I,J) which is ! bounded by the corner nodes [I,J], [I+1,J], [I,J+1], and [I+1,J+1]. ! ! Output, double precision AREAL, the area of the liquid region. ! ! Output, double precision AREAS. ! AREAS is the area of the solid or crystal region. ! ! Output, double precision AREAT, the total area. ! ! Input, integer ICRYS. ! ICRYS specifies the end of the crystal in the I array direction, ! and in the vertical coordinate direction. ! ! Input, integer JCRYS. ! JCRYS specifies the end of the crystal in the J array direction, ! and in the horizontal coordinate direction. ! ! Input, integer L0. ! L0 is the extent of the grid in the vertical or "I" coordinate. ! ! Input, integer M0. ! M0 is the extent of the grid in the horizontal or "J" coordinate. ! ! Input, integer NI. ! NI is the maximum number of grid points in the "I" or vertical ! direction. ! ! Input, integer NJ. ! NJ is the maximum number of grid points in the "J" or horizontal ! direction. ! ! Input, double precision XC(NI,NJ). ! XC contains the X coordinate of nodes which are the "corners" of ! control volumes. ! ! Input, double precision YC(NI,NJ). ! YC contains the Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer ni integer nj ! double precision area(ni,nj) double precision areal double precision areas double precision areat integer i integer icrys integer j integer jcrys integer l0 integer m0 integer nbad double precision xc(ni,nj) double precision yc(ni,nj) ! do i = 1,l0 do j = 1,m0 if ( i == 1 .or. i == l0.or.j == 1.or.j == m0 ) then area(i,j) = 0.0D+00 else area(i,j) = 0.5D+00*( & xc(i, j)* (yc(i+1,j) -yc(i, j+1)) & +xc(i+1,j)* (yc(i+1,j+1)-yc(i, j)) & +xc(i+1,j+1)*(yc(i, j+1)-yc(i+1,j)) & +xc(i, j+1)*(yc(i, j) -yc(i+1,j+1))) end if end do end do ! ! Check for illegal zero length sides. ! nbad = 0 do i = 2,l0-1 do j = 2,m0-1 if ( xc(i,j) == xc(i+1,j) .and. yc(i,j) == yc(i+1,j) ) then nbad = nbad+1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,2i6)' ) ' Zero length side for cell I,J:',i,j write ( *, '(a,2g14.6)' ) & 'XC,YC(I,J) = XC,YC(I+1,J)=',xc(i,j),yc(i,j) else if ( xc(i+1,j) == xc(i+1,j+1) .and. yc(i+1,j) == yc(i+1,j+1) ) then nbad = nbad+1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,2i6)' ) ' Zero length side for cell I,J:',i,j write ( *, '(a,2g14.6)' ) & 'XC,YC(I+1,J) = XC,YC(I+1,J+1)=',xc(i+1,j),yc(i+1,j) else if ( xc(i+1,j+1) == xc(i,j+1) .and. yc(i+1,j+1) == yc(i,j+1) ) then nbad = nbad+1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,2i6)' ) ' Zero length side for cell I,J:',i,j write ( *, '(a,2g14.6)' ) & 'XC,YC(I+1,J+1) = XC,YC(I,J+1)=',xc(i+1,j+1),yc(i+1,j+1) else if ( xc(i,j+1) == xc(i,j) .and. yc(i,j+1) == yc(i,j) ) then nbad = nbad+1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,2i6)' ) ' Zero length side for cell I,J:',i,j write ( *, '(a,2g14.6)' ) & 'XC,YC(I,J+1) = XC,YC(I,J)=',xc(i,j+1),yc(i,j+1) end if end do end do if ( nbad > 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,i6,a)' ) ' A total of ',nbad,' zero length cell sides.' stop end if ! ! Check for illegal zero area cells. ! nbad = 0 do i = 2,l0-1 do j = 2,m0-1 if ( area(i,j) == 0.0D+00 ) then nbad = nbad+1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,i6,a,i6)' ) ' Zero area for cell I = ',i,' J=',j end if end do end do if ( nbad > 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,i6,a)' ) ' A total of ',nbad,' null cells.' stop end if ! ! Compute liquid, solid, and total areas. ! areal = 0.0 do i = 1,icrys-1 do j = 1,m0 areal = areal+area(i,j) end do end do areas = 0.0 do i = icrys,l0 do j = 1,jcrys areas = areas+area(i,j) end do end do areat = 0.0 do i = 1,l0 do j = 1,m0 areat = areat+area(i,j) end do end do return end subroutine findp ( j, xold, yold, xnew, ynew, xc, yc ) ! !*****************************************************************************80 ! !! FINDP finds the boundary nodes associated with a Neumann condition. ! ! ! Discussion: ! ! FINDP is given a point (XOLD,YOLD). ! ! FINDP returns a pair of values (XNEW,YNEW), which represent ??? ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! double precision a1 double precision ai double precision aj double precision alp double precision b1 double precision bet double precision c1 double precision dels double precision denm double precision dxb1 double precision dxb2 double precision dyb1 double precision dyb2 integer j double precision r0 double precision r1 double precision r11 double precision r2 double precision r22 double precision r33 double precision slpi double precision xc(ni,nj) double precision xnew double precision xnew1 double precision xnew2 double precision xold double precision yc(ni,nj) double precision yn1 double precision yn2 double precision ynew double precision yold ! ! Consider the three points P, Q and R, which have a constant ! KSI coordinate, and an increasing ETA coordinate: ! ! P = ( XC(2,J-1), YC(2,J-1) ) ! Q = ( XC(2,J), YC(2,J) ) ! R = ( XC(2,J+1), YC(2,J+1) ) ! ! ! ^ R = [2,J+1] ! | ! E ! T Q = [2,J] ! A ! | ! | P = [2,J-1] ! | ! +------KSI-------------> ! ! The slope of the line from P to Q is (YQ-YP)/(XQ-XP), and ! the slope of the line from Q to R is (YR-YQ)/(XR-XQ). ! ! If these slopes are close enough, we can assume the three points ! lie on a straight line. So look at the size of ! DENM = (YQ-YP)*(XR-XQ)-(YR-YQ)*(XQ-XP). ! dxb1 = xc(2,j)-xc(2,j-1) dxb2 = xc(2,j+1)-xc(2,j) dyb1 = yc(2,j)-yc(2,j-1) dyb2 = yc(2,j+1)-yc(2,j) denm = dyb2*dxb1-dyb1*dxb2 ! if ( abs(denm) < 0.0001D+00 ) then if ( abs(dxb1) >= 0.001D+00 ) then slpi = -dyb1/dxb1 aj = yc(2,j)+slpi*xc(2,j) ai = xold-slpi*yold ynew = (aj-ai*slpi)/(1.0+slpi*slpi) xnew = slpi*ynew+ai else ynew = yold xnew = xc(2,j) end if ! ! Use 3 points to find a circle. ! (X-ALP)**2+(Y-BET)**2 = R0 ! ! ! Find the cross point between circle and straight line ! x = ai+slpi*y get a1*y**2-2*b1*y+c1=0. ! else r11 = xc(2,j-1)**2+yc(2,j-1)**2 r22 = xc(2,j)**2+yc(2,j)**2 r33 = xc(2,j+1)**2+yc(2,j+1)**2 r1 = 0.5D+00*(r22-r11) r2 = 0.5D+00*(r33-r22) alp = (r1*dyb2-r2*dyb1)/denm bet = (r2*dxb1-r1*dxb2)/denm r0 = (alp-xc(2,j))**2+(bet-yc(2,j))**2 slpi = (xold-alp)/(yold-bet) ai = alp-slpi*bet a1 = 1.0D+00+slpi*slpi b1 = slpi*(alp-ai)+bet c1 = (alp-ai)**2+bet*bet-r0 dels = sqrt(b1*b1-a1*c1) yn1 = (b1+dels)/a1 yn2 = (b1-dels)/a1 xnew1 = slpi*yn1+ai xnew2 = slpi*yn2+ai r1 = (xnew1-xc(2,j))**2+(yn1-yc(2,j))**2 r2 = (xnew2-xc(2,j))**2+(yn2-yc(2,j))**2 if ( r1 < r2 ) then ynew = yn1 else ynew = yn2 end if xnew = slpi*ynew+ai end if return end subroutine flux(ae1,ae2,aim,aip,ajm,ajp,ak1,ak2,ap,b1jbl,b2jbl, & b3jbl,cappa,cd,cmu,con,epsil,f,fjeta,fjksi,fmax,fn,fo,gam,heta, & hksi,icrys,isol,iter,izone,jcrys,l1,lblk,lconv,lortho,m1,mode,nf, & nsolve,ntimes,r,relax,res,rueta,ruksi) ! !*****************************************************************************80 ! !! FLUX computes the flux of various quantities. ! ! ! It looks like if ISOL is 0, you just go through the big loop once. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nk = 14 integer, parameter :: nmaxij = 64 integer, parameter :: ns = 10 ! double precision acof double precision ae1(ni,nj) double precision ae2(ni,nj) double precision aim(ni,nj) double precision aip(ni,nj) double precision ajm(ni,nj) double precision ajp(ni,nj) double precision ak1(ni,nj) double precision ak2(ni,nj) double precision ap(ni,nj) double precision ap0(ni,nj) double precision b1jbl(nj) double precision b2jbl(nj) double precision b3jbl(ni) double precision cappa double precision cd double precision cmu double precision con(ni,nj) double precision con0(ni,nj) double precision diff double precision epsil double precision error double precision f(ni,nj,ns) double precision fjbi1(nj,ns) double precision fjbj1(ni,ns) double precision fjbl1(nj,ns) double precision fjbm1(ni,ns) double precision fjeta(ni,nj) double precision fjksi(ni,nj) double precision flow double precision fmax(ns) double precision fn(ni,nj,ns) double precision fo(ni,nj,ns) double precision gam(ni,nj) double precision heta(ni,nj) double precision hksi(ni,nj) integer i integer ii integer icrys integer isol integer iter integer izone integer j integer jcrys integer jj integer l1 logical lblk(ns) logical lconv logical lortho integer m1 integer mode integer nf integer nsolve(ns) integer nt integer ntimes(ns) double precision pt(nmaxij) double precision qt(nmaxij) double precision r(ni,nj) double precision relax(nk) double precision res(ns) double precision rueta(ni,nj) double precision ruksi(ni,nj) double precision t1 double precision t2 double precision t3 double precision t4 double precision tem1 double precision tem2 double precision xje double precision xjn double precision xjw double precision xjs double precision xme double precision xmn double precision xmw double precision xms ! ! Save copies of CON and AP. ! do i = 1,l1 do j = 1,m1 con0(i,j) = con(i,j) ap0(i,j) = ap(i,j) end do end do ! ! I think ITER = 0 occurs only for Temperature in the Solid zone... ! if ( iter == 0 ) then if ( isol /= 0 ) then call solve1(aim,aip,ajm,ajp,ap,cappa,cd,cmu,con,f,fo, & heta,hksi,icrys,izone,jcrys,l1,m1,nf) call solve2(aim,aip,ajm,ajp,ap,con,f,l1,lblk,m1,nf,nsolve) end if return end if do nt = 1,ntimes(nf) ! ! Find the flux on the first step only. ! if ( nt == 1 ) then do j = 1,m1 diff = gam(2,j)/(0.5D+00*hksi(2,j)) flow = -ruksi(2,j) call diflow(acof,diff,flow) qt(2) = ruksi(2,j)*f(2,j,nf)+acof*(f(1,j,nf)-f(2,j,nf)) do i = 3,l1-1 jj = j if ( j == 1) jj = 2 if ( j == m1) jj = m1-1 diff = gam(i,jj)*gam(i-1,jj)/(0.5D+00*hksi(i,jj)*gam(i-1,jj) & +0.5D+00*hksi(i-1,jj)*gam(i,jj)) flow = -ruksi(i,j) call diflow(acof,diff,flow) qt(i) = ruksi(i,j)*f(i,j,nf)+acof*(f(i-1,j,nf)-f(i,j,nf)) end do diff = gam(l1-1,j)/(0.5D+00*hksi(l1-1,j)) flow = ruksi(l1,j) call diflow(acof,diff,flow) qt(l1) = ruksi(l1,j)*f(l1-1,j,nf)+acof*(f(l1-1,j,nf)-f(l1,j,nf)) do i = 2,l1 fjksi(i,j) = qt(i) end do end do ! ! Along J. ! do i = 1,l1 diff = gam(i,2)/(0.5D+00*heta(i,2)) flow = -rueta(i,2) call diflow(acof,diff,flow) qt(2) = rueta(i,2)*f(i,2,nf)+acof*(f(i,1,nf)-f(i,2,nf)) do j = 3,m1-1 if ( i == 1 ) then ii = 2 else if ( i == l1 ) then ii = l1-1 else ii = i end if diff = gam(ii,j)*gam(ii,j-1)/(0.5D+00*heta(ii,j)*gam(ii,j-1) & +0.5D+00*heta(ii,j-1)*gam(ii,j)) flow = -rueta(i,j) call diflow(acof,diff,flow) qt(j) = rueta(i,j)*f(i,j,nf)+acof*(f(i,j-1,nf)-f(i,j,nf)) end do diff = gam(i,m1-1)/(0.5D+00*heta(i,m1-1)) flow = rueta(i,m1) call diflow(acof,diff,flow) qt(m1) = rueta(i,m1)*f(i,m1-1,nf)+acof*(f(i,m1-1,nf)-f(i,m1,nf)) do j = 2,m1 fjeta(i,j) = qt(j) end do end do ! ! Finish JKSI and JETA calculation. ! Construct boundary flux for J type boundary condition. ! Boundary condition con(1,j),*,*,* go into jksi and jeta ! ak2/ak1,* is project on xi - direction from eta - direction ! variable, in program let ak2 = 0.0, because sym. orth. ! do j = 2,m1-1 if ( gam(1,j) == 0.0D+00 ) then fjksi(2,j) = (con(1,j)*heta(1,j)+0.5D+00*ak2(2,j) & *(fjeta(1,j)+fjeta(1,j+1)))/ak1(2,j) end if if ( gam(l1,j) == 0.0D+00 ) then fjksi(l1,j) = (con(l1,j)*heta(l1,j)+0.5D+00*ak2(l1,j) & *(fjeta(l1,j)+fjeta(l1,j+1)))/ak1(l1,j) end if end do do i = 2,l1-1 if ( gam(i,1) == 0.0D+00 ) then fjeta(i,2) = (con(i,1)*hksi(i,1)+0.5D+00*ae2(i,2) & *(fjksi(i,1)+fjksi(i+1,1)))/ae1(i,2) end if if ( gam(i,m1) == 0.0D+00 ) then fjeta(i,m1) = (con(i,m1)*hksi(i,m1)+0.5D+00*ae2(i,m1) & *(fjksi(i,m1)+fjksi(i+1,m1)))/ae1(i,m1) end if end do ! ! (2-31) and (2-42) find out vector J cdot vector n or J_n in bound. ! do i = 2,l1-1 t1 = 0.5D+00*(fjksi(i,1)+fjksi(i+1,1)) t2 = 0.5D+00*(fjksi(i,m1)+fjksi(i+1,m1)) fjbj1(i,nf) = (fjeta(i,2)*ae1(i,2)-t1*ae2(i,2))/hksi(i,1) fjbm1(i,nf) = (fjeta(i,m1)*ae1(i,m1)-t2*ae2(i,m1))/hksi(i,m1) if ( mode == 1 ) then fjbj1(i,nf) = fjbj1(i,nf)/r(i,2) fjbm1(i,nf) = fjbm1(i,nf)/r(i,m1) end if end do do j = 2,m1-1 t1 = 0.5D+00*(fjeta(1,j)+fjeta(1,j+1)) t2 = 0.5D+00*(fjeta(l1,j)+fjeta(l1,j+1)) fjbi1(j,nf) = (fjksi(2,j)*ak1(2,j)-t1*ak2(2,j))/heta(1,j) fjbl1(j,nf) = (fjksi(l1,j)*ak1(l1,j)-t2*ak2(l1,j))/heta(l1,j) if ( mode == 1 ) then fjbi1(j,nf) = fjbi1(j,nf)/r(2,j) fjbl1(j,nf) = fjbl1(j,nf)/r(l1,j) end if end do ! ! This code is only for the temperature variable. ! We compute ! B1JBL, the heat flux between rows (?) ICRYS-1 and ICRYS, ! for columns 2 through JCRYS and column M1. ! B2JBL, the heat flux between rows (?) ICRYS and ICRYS+1, ! for columns 2 through JCRYS and column M1. ! B3JBL, the heat flux between rows (?) and (?), ! for columns JCRYS+1 to M1-1. ! if ( nf == 5 ) then do j = 2,jcrys-1 t1 = gam(icrys-1,j)*(f(icrys-1,j,5)-f(icrys,j,5))/hksi(icrys-1,j) b1jbl(j) = t1*ak1(icrys,j)/heta(icrys,j) if ( mode == 1 ) then b1jbl(j) = b1jbl(j)/r(icrys,j) end if end do b1jbl(jcrys) = b1jbl(jcrys+1) b1jbl(m1) = b1jbl(m1-1) do j = 2,jcrys-1 t1 = gam(icrys+1,j)*(f(icrys,j,5)-f(icrys+1,j,5))/hksi(icrys+1,j) b2jbl(j) = t1*ak1(icrys,j)/heta(icrys,j) if ( mode == 1 ) then b2jbl(j) = b2jbl(j)/r(icrys,j) end if end do b2jbl(jcrys) = b2jbl(jcrys+1) b2jbl(m1) = b2jbl(m1-1) do j = jcrys+1,m1-1 t1 = gam(icrys-1,j)*(f(icrys-1,j,5)-f(icrys,j,5))/hksi(icrys-1,j) t2 = gam(icrys-1,j)*(f(icrys,j+1,5)-f(icrys,j,5))/heta(icrys,j) b3jbl(j) = (t1*ak1(icrys,j)-t2*ak2(icrys,j))/heta(icrys,j) if ( mode == 1 ) then b3jbl(j) = b3jbl(j)/r(icrys,j) end if end do b3jbl(jcrys) = b3jbl(jcrys+1) b3jbl(m1) = b3jbl(m1-1) end if else ! ! This code is done if this is NOT the first iteration. ! ! Correct J from known PHI and J. ! Correct JKSI. ! do j = 2,m1-1 if ( gam(1,j) /= 0.0D+00 ) then fjksi(2,j) = fjksi(2,j)+aim(2,j)/ak1(2,j)* & (f(1,j,nf)-f(2,j,nf))+ruksi(2,j)*f(2,j,nf) end if do i = 3,l1 if ( gam(l1,j) /= 0.0D+00 .or. i.ne.l1 ) then fjksi(i,j) = fjksi(i,j)+aip(i-1,j)/ak1(i,j)*(f(i-1,j,nf) & -f(i,j,nf))+ruksi(i,j)*f(i-1,j,nf) end if end do end do ! ! Correct FJETA. ! do i = 2,l1-1 if ( gam(i,1) /= 0.0D+00 ) then fjeta(i,2) = fjeta(i,2)+ajm(i,2)/ae1(i,2)* & (f(i,1,nf)-f(i,2,nf))+rueta(i,2)*f(i,2,nf) end if do j = 3,m1 if ( gam(i,m1) /= 0.0D+00 .or. j.ne.m1 ) then fjeta(i,j) = fjeta(i,j)+ajp(i,j-1)/ae1(i,j)* & (f(i,j-1,nf)-f(i,j,nf))+rueta(i,j)*f(i,j-1,nf) end if end do end do end if ! ! End of "Is this the first iteration or not" block. ! ! Restore old values of CON and AP. ! do i = 1,l1 do j = 1,m1 con(i,j) = con0(i,j) ap(i,j) = ap0(i,j) end do end do if ( .not.lortho ) then do j = 2,m1-1 do i = 1,l1 pt(i) = 0.5D+00*(fjeta(i,j)+fjeta(i,j+1)) qt(i) = 0.5D+00*(rueta(i,j)+rueta(i,j+1)) end do do i = 2,l1-1 t1 = hksi(i,j)/(hksi(i+1,j)+hksi(i,j)) t2 = hksi(i+1,j)/(hksi(i+1,j)+hksi(i,j)) t3 = hksi(i,j)/(hksi(i-1,j)+hksi(i,j)) t4 = hksi(i-1,j)/(hksi(i-1,j)+hksi(i,j)) xje = t2*pt(i)+t1*pt(i+1) xme = t2*qt(i)+t1*qt(i+1) xjw = t4*pt(i)+t3*pt(i-1) xmw = t4*qt(i)+t3*qt(i-1) con(i,j) = con(i,j)+(xje-f(i,j,nf)*xme)*ak2(i+1,j) & -(xjw-f(i,j,nf)*xmw)*ak2(i,j) end do end do do i = 2,l1-1 do j = 1,m1 pt(j) = 0.5D+00*(fjksi(i,j)+fjksi(i+1,j)) qt(j) = 0.5D+00*(ruksi(i,j)+ruksi(i+1,j)) end do do j = 2,m1-1 tem1 = heta(i,j+1)+heta(i,j) tem2 = heta(i,j-1)+heta(i,j) t1 = heta(i,j)/tem1 t2 = heta(i,j+1)/tem1 t3 = heta(i,j)/tem2 t4 = heta(i,j-1)/tem2 xjn = t2*pt(j)+t1*pt(j+1) xmn = t2*qt(j)+t1*qt(j+1) xjs = t4*pt(j)+t3*pt(j-1) xms = t4*qt(j)+t3*qt(j-1) con(i,j) = con(i,j)+(xjn-f(i,j,nf)*xmn)*ae2(i,j+1) & -(xjs-f(i,j,nf)*xms)*ae2(i,j) end do end do end if ! ! Calculate JHAT. ! Step 5 use equation (2-97) and (2-102), (2-103), (2-104) ! calculate boundary data from J and other purpose for B_SP. ! from here to end. pc = JKSI HAT or JETA HAT. ! ! jksi hat ! do j = 2,m1-1 f(2,j,4) = 0.0D+00 if ( gam(1,j) == 0.0D+00 ) then t1 = fjksi(2,j)-ruksi(2,j)*f(2,j,nf) diff = gam(2,j)/(0.5D+00*hksi(2,j)) flow = -ruksi(2,j) call diflow(acof,diff,flow) if ( acof /= 0.0D+00 ) then f(1,j,nf) = f(2,j,nf)+t1/acof else f(1,j,nf) = f(2,j,nf) end if f(2,j,4) = fjksi(2,j) end if f(l1,j,4) = 0.0D+00 if ( gam(l1,j) == 0.0D+00 ) then t1 = fjksi(l1,j)-ruksi(l1,j)*f(l1-1,j,nf) diff = gam(l1-1,j)/(0.5D+00*hksi(l1-1,j)) flow = ruksi(l1,j) call diflow(acof,diff,flow) if ( acof /= 0.0D+00 ) then f(l1,j,nf) = f(l1-1,j,nf)-t1/acof else f(l1,j,nf) = f(l1-1,j,nf) end if f(l1,j,4) = fjksi(l1,j) end if do i = 3,l1-1 f(i,j,4) = 0.0D+00 end do end do do j = 2,m1-1 do i = 2,l1 fjksi(i,j) = f(i,j,4) end do end do ! ! jeta hat ! do i = 2,l1-1 f(i,2,4) = 0.0D+00 if ( gam(i,1) == 0.0D+00 ) then t1 = fjeta(i,2)-rueta(i,2)*f(i,2,nf) diff = gam(i,2)/(0.5D+00*heta(i,2)) flow = -rueta(i,2) call diflow(acof,diff,flow) if ( acof /= 0.0D+00 ) then f(i,1,nf) = f(i,2,nf)+t1/acof else f(i,1,nf) = f(i,2,nf) end if f(i,2,4) = fjeta(i,2) end if f(i,m1,4) = 0.0D+00 if ( gam(i,m1) == 0.0D+00 ) then t1 = fjeta(i,m1)-rueta(i,m1)*f(i,m1-1,nf) diff = gam(i,m1-1)/(0.5D+00*heta(i,m1-1)) flow = rueta(i,m1) call diflow(acof,diff,flow) if ( acof /= 0.0D+00 ) then f(i,m1,nf) = f(i,m1-1,nf)-t1/acof else f(i,m1,nf) = f(i,m1-1,nf) end if f(i,m1,4) = fjeta(i,m1) end if do j = 3,m1-1 f(i,j,4) = 0.0D+00 end do end do do i = 2,l1-1 do j = 2,m1 fjeta(i,j) = f(i,j,4) end do end do ! ! Calculate and solve phy equation ! (2-106) and (2-107) b = b_s + b_no + b_sp, t1=b_sp ! do i = 2,l1-1 do j = 2,m1-1 con(i,j) = con(i,j)+fjksi(i,j)*ak1(i,j)-fjksi(i+1,j)*ak1(i+1,j) & +fjeta(i,j)*ae1(i,j)-fjeta(i,j+1)*ae1(i,j+1) end do end do if ( isol == 0) return do i = 2,l1-1 do j = 2,m1-1 ap(i,j) = ap(i,j)/relax(nf) con(i,j) = con(i,j)+(1.0D+00-relax(nf))*ap(i,j)*f(i,j,nf) end do end do call solve1(aim,aip,ajm,ajp,ap,cappa,cd,cmu,con,f,fo, & heta,hksi,icrys,izone,jcrys,l1,m1,nf) call solve2(aim,aip,ajm,ajp,ap,con,f,l1,lblk,m1,nf,nsolve) ! ! BEGIN NEW. ! ! Compute the largest solution component. ! fmax(nf) = 0.0 do i = 1,l1 do j = 1,m1 fmax(nf) = max(fmax(nf),abs(f(i,j,nf))) end do end do ! ! Compare current and previous iterates. ! if ( nt > 1 ) then res(nf) = 0.0D+00 do i = 1,l1 do j = 1,m1 error = abs(f(i,j,nf)-fn(i,j,nf)) if ( fmax(nf) > 0.0D+00 ) then error = error/fmax(nf) end if res(nf) = max(res(nf),error) end do end do end if ! ! Save the new solution in FN. ! do i = 1,l1 do j = 1,m1 fn(i,j,nf) = f(i,j,nf) end do end do if ( res(nf) < epsil .and. nt > 1 ) then ! ! write ( *, '(a,i6)' ) 'FLUX - Convergence for variable ',nf ! write ( *, '(a,i6)' ) ' on step ',nt ! write ( *, '(a,g14.6)' ) ' RES(NF) = ',res(nf) ! write ( *, '(a,g14.6)' ) ' FMAX(NF)=',fmax(nf) ! return end if end do ! ! End of big iteration ! ! Compute the largest solution component. ! ! fmax(nf) = 0.0D+00 ! do i = 1,l1 ! do j = 1,m1 ! fmax(nf) = max(fmax(nf),abs(f(i,j,nf))) ! end do ! end do ! ! Compute the largest relative change in the solution. ! ! res(nf) = 0.0D+00 ! do i = 1,l1 ! do j = 1,m1 ! ! error = abs(f(i,j,nf)-fn(i,j,nf)) ! ! if ( fmax(nf) > 0.0D+00 ) then ! error = error/fmax(nf) ! end if ! ! res(nf) = max(res(nf),error) ! ! end do ! end do ! ! Decide whether to declare convergence or not. ! if ( res(nf) > epsil ) then lconv = .false. end if ! ! Save the new solution in FN. ! ! do i = 1,l1 ! do j = 1,m1 ! fn(i,j,nf) = f(i,j,nf) ! end do ! end do ! nt = ntimes(nf) ! write ( *, '(a,i6)' ) 'FLUX - No convergence for variable ',nf ! write ( *, '(a,i6)' ) ' on step ',nt ! write ( *, '(a,g14.6)' ) ' RES(NF) = ',res(nf) ! write ( *, '(a,g14.6)' ) ' FMAX(NF)=',fmax(nf) return end subroutine gamsor(ap,birad,ce1,ce2,cmu,con,ewall,f,fcsl,fksl, & fma,fo,frsl,gam,gamt,grash,hamag,heta,hksi,icrys,inturb, & izone,jcrys,l1,lsolve,m1,mode,nf,pr,r,rdtm,re,recb,rect, & rho,rhocon,rpr,sige,sigk,sigt,stel,stes,tal,tas,tf,tw,x, & xc,y,yc) ! !*****************************************************************************80 ! !! GAMSOR sets the coefficients for the transport problems. ! ! ! Discussion: ! ! These coefficients are associated with U, V, T, W, TK, TE and E. ! ! Note that GAM(I,1) = 0 means the gradient is zero at wall J=1. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nmaxij = 64 integer, parameter :: ns = 10 ! double precision ap(ni,nj) double precision ar double precision birad double precision ce1 double precision ce2 double precision cm4 double precision cmu double precision con(ni,nj) double precision dcdr double precision depdr double precision depdx double precision dudx(ni,nj) double precision dudy(ni,nj) double precision dvdx(ni,nj) double precision dvdy(ni,nj) double precision dwdx(ni,nj) double precision dwdy(ni,nj) double precision ewall double precision f(ni,nj,ns) double precision fcsl double precision fksl double precision fma double precision fo(ni,nj,ns) double precision frsl double precision gam(ni,nj) double precision gamt(ni,nj) double precision gdudx(ni,nj) double precision gdudy(ni,nj) double precision gdvdx(ni,nj) double precision gdvdy(ni,nj) double precision grash double precision hamag double precision heta(ni,nj) double precision hksi(ni,nj) integer i integer icrys integer inturb integer izone integer j integer jcrys integer l1 logical lsolve(ns) integer m1 integer mode integer modelm integer nf double precision p11 double precision p12 double precision p21 double precision p22 double precision p31 double precision p32 double precision p41 double precision p42 double precision pbig double precision pr double precision prod(ni,nj) double precision r(ni,nj) double precision rdtm double precision re double precision recb double precision rect double precision rho(ni,nj) double precision rhocon double precision rpr double precision sige double precision sigk double precision sigt double precision stel double precision stes double precision tal double precision tas double precision tf double precision tw double precision x(ni,nj) double precision xc(ni,nj) double precision xplus(nmaxij) double precision y(ni,nj) double precision yc(ni,nj) double precision yplus(nmaxij) ! do i = 1,l1 do j = 1,m1 ap(i,j) = -rdtm*rho(i,j) con(i,j) = 0.0D+00 gam(i,j) = rhocon/re+gamt(i,j) end do end do ! ! Horizontal velocity. ! if ( nf == 1 ) then do i = 2,l1-1 do j = 2,m1-1 con(i,j) = fo(i,j,1)*rho(i,j)*rdtm+rho(i,j)*f(i,j,5)*grash/re**2 end do end do if ( inturb == 1 ) then do i = 1,l1 do j = 1,m1 gdudx(i,j) = 0.0D+00 gdudy(i,j) = 0.0D+00 gdvdx(i,j) = 0.0D+00 gdvdy(i,j) = 0.0D+00 end do end do do i = 2,l1-1 do j = 2,m1-1 call gradnt(heta,hksi,i,j,f(1,1,1),dudx(i,j),dudy(i,j),xc,yc) call gradnt(heta,hksi,i,j,f(1,1,2),dvdx(i,j),dvdy(i,j),xc,yc) gdudx(i,j) = gamt(i,j)*dudx(i,j) gdudy(i,j) = gamt(i,j)*dudy(i,j) gdvdx(i,j) = gamt(i,j)*dvdx(i,j) gdvdy(i,j) = gamt(i,j)*dvdy(i,j) end do end do do i = 2,l1-1 do j = 2,m1-1 call gradnt(heta,hksi,i,j,gdudx,p11,p12,xc,yc) call gradnt(heta,hksi,i,j,gdudy,p21,p22,xc,yc) call gradnt(heta,hksi,i,j,gdvdx,p31,p32,xc,yc) call gradnt(heta,hksi,i,j,gdvdy,p41,p42,xc,yc) con(i,j) = con(i,j)+p11+p32 end do end do cm4 = cmu**0.25 if ( mode == 0 ) then do i = 2,l1-1 yplus(i) = re*rho(i,2)*sqrt(f(i,2,7))*cm4*0.5*heta(i,2)/rhocon if ( yplus(i) <= 11.63D+00 ) then gam(i,1) = rhocon/re else gam(i,1) = yplus(i)/(2.5D+00*log(ewall*yplus(i))) end if end do end if do i = 2,l1-1 yplus(i) = re*rho(i,m1-1)*sqrt(f(i,m1-1,7)) & *cm4*0.5D+00*heta(i,m1-1)/rhocon if ( yplus(i) <= 11.63D+00 ) then gam(i,m1) = rhocon/re else gam(i,m1) = yplus(i)/(2.5D+00*log(ewall*yplus(i))) end if end do do j = 2,m1-1 xplus(j) = re*rho(2,j)*sqrt(f(2,j,7))*cm4*0.5D+00*hksi(2,j)/rhocon if ( xplus(j) <= 11.63D+00 ) then gam(1,j) = rhocon/re else gam(1,j) = xplus(j)/(2.5D+00*log(ewall*xplus(j))) end if end do do j = 2,m1-1 xplus(j) = re*rho(l1-1,j)*sqrt(f(l1-1,j,7)) & *cm4*0.5D+00*hksi(l1-1,j)/rhocon if ( xplus(j) <= 11.63D+00 ) then gam(l1,j) = rhocon/re else gam(l1,j) = xplus(j)/(2.5D+00*log(ewall*xplus(j))) end if end do end if do i = 2,l1-1 gam(i,1) = 0.0D+00 end do do j = 2,m1-1 f(1,j,1) = 0.0D+00 f(l1,j,1) = 0.0D+00 end do f(l1,1,1) = 0.0D+00 do j = jcrys+1,m1-1 f(l1,j,1) = fo(l1,j,2)*(xc(l1,j+1)-xc(l1,j))/(yc(l1,j+1)-yc(l1,j)) end do f(l1,jcrys+1,1) = 0.0D+00 ! ! Vertical velocity. ! else if ( nf == 2 ) then do i = 2,l1-1 do j = 2,m1-1 con(i,j) = fo(i,j,2)*rho(i,j)*rdtm ap(i,j) = ap(i,j)-rho(i,j)*hamag**2/re if ( mode == 1 ) then con(i,j) = con(i,j)+rho(i,j)*f(i,j,6)**2/r(i,j)**3 ap(i,j) = ap(i,j)-rho(i,j)/r(i,j)**2 end if end do end do if ( inturb == 1 ) then do i = 1,l1 do j = 1,m1 gdudx(i,j) = 0.0D+00 gdudy(i,j) = 0.0D+00 gdvdx(i,j) = 0.0D+00 gdvdy(i,j) = 0.0D+00 end do end do do i = 2,l1-1 do j = 2,m1-1 call gradnt(heta,hksi,i,j,f(1,1,1),dudx(i,j),dudy(i,j),xc,yc) call gradnt(heta,hksi,i,j,f(1,1,2),dvdx(i,j),dvdy(i,j),xc,yc) gdudx(i,j) = gamt(i,j)*dudx(i,j) gdudy(i,j) = gamt(i,j)*dudy(i,j) gdvdx(i,j) = gamt(i,j)*dvdx(i,j) gdvdy(i,j) = gamt(i,j)*dvdy(i,j) end do end do do i = 2,l1-1 do j = 2,m1-1 call gradnt(heta,hksi,i,j,gdudx,p11,p12,xc,yc) call gradnt(heta,hksi,i,j,gdudy,p21,p22,xc,yc) call gradnt(heta,hksi,i,j,gdvdx,p31,p32,xc,yc) call gradnt(heta,hksi,i,j,gdvdy,p41,p42,xc,yc) con(i,j) = con(i,j)+p21+p42 end do end do cm4 = cmu**0.25 if ( mode == 0 ) then do i = 2,l1-1 yplus(i) = re*rho(i,2)*sqrt(f(i,2,7))*cm4*0.5D+00*heta(i,2)/rhocon gam(i,1) = rhocon/re if ( yplus(i) > 11.63D+00) then gam(i,1) = yplus(i)/(2.5D+00*log(ewall*yplus(i))) end if end do end if do i = 2,l1-1 yplus(i) = re*rho(i,m1-1)*sqrt(f(i,m1-1,7)) & *cm4*0.5*heta(i,m1-1)/rhocon gam(i,m1) = rhocon/re if ( yplus(i) > 11.63D+00 ) then gam(i,m1) = yplus(i)/(2.5D+00*log(ewall*yplus(i))) end if end do do j = 2,m1-1 xplus(j) = re*rho(2,j)*sqrt(f(2,j,7))*cm4*0.5D+00*hksi(2,j)/rhocon gam(1,j) = rhocon/re if ( xplus(j) > 11.63D+00 ) then gam(1,j) = xplus(j)/(2.5D+00*log(ewall*xplus(j))) end if end do do j = 2,m1-1 xplus(j) = re*rho(l1-1,j)*sqrt(f(l1-1,j,7)) & *cm4*0.5D+00*hksi(l1-1,j)/rhocon gam(l1,j) = rhocon/re if ( xplus(j) > 11.63D+00 ) then gam(l1,j) = xplus(j)/(2.5D+00*log(ewall*xplus(j))) end if end do end if do i = 2,l1-1 f(i,1,2) = 0.0D+00 end do do j = jcrys+1,m1-1 dcdr = (f(l1,j+1,5)-f(l1,j,5))/(y(l1,j+1)-y(l1,j)) f(l1,j,2) = f(l1-1,j,2)+(xc(l1,j)-xc(l1-1,j))*dcdr*(fma/(re*pr)) end do f(l1,jcrys+1,2) = 0.0D+00 do j = 2,jcrys f(1,j,2) = 0.0D+00 f(l1,j,2) = 0.0D+00 end do f(l1,m1-1,2) = 0.0D+00 f(l1,m1,2) = 0.0D+00 ! ! Temperature computations in the solid zone. ! else if ( nf == 5 .and. izone == 1 ) then do i = 1,l1 do j = 1,m1 gam(i,j) = rpr*fksl/fcsl/frsl end do end do do i = icrys+1,l1 gam(i,1) = 0.0D+00 end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = fo(i,j,5)*rho(i,j)*rdtm end do end do do j = 2,jcrys f(l1,j,5) = -stes/stel end do do i = icrys+1,l1 f(i,m1,5) = f(i,m1-1,5)-birad*((f(i,m1,5)*(tw-tf)+tf)**4-tas**4) & *(y(i,m1)-y(i,m1-1))/100000000.0D+00 end do ! ! Temperature calculations in the liquid zone. ! else if ( nf == 5 .and. izone == 2 ) then do i = 1,l1 do j = 1,m1 gam(i,j) = rpr+gamt(i,j)/sigt end do end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = fo(i,j,5)*rho(i,j)*rdtm end do end do if ( inturb == 1 ) then cm4 = cmu**0.25 pbig = 9.24D+00*((pr/sigt)**0.75-1.0D+00) & *(1.0D+00+0.28D+00*exp(-0.007D+00*pr/sigt)) if ( mode == 0) then do i = 2,l1-1 yplus(i) = re*rho(i,2)*sqrt(f(i,2,7))*cm4*0.5D+00*heta(i,2)/rhocon gam(i,1) = rpr if ( yplus(i) > 11.63D+00 ) then gam(i,1) = max(rhocon/re, & yplus(i)/(sigt*(2.5D+00*log(ewall*yplus(i))+pbig))) end if end do end if do i = 2,l1-1 yplus(i) = re*rho(i,m1-1)*sqrt(f(i,m1-1,7)) & *cm4*0.5D+00*heta(i,m1-1)/rhocon gam(i,m1) = rpr if ( yplus(i) > 11.63D+00 ) then gam(i,m1) = max(rhocon/re, & yplus(i)/(sigt*(2.5D+00*log(ewall*yplus(i))+pbig))) end if end do do j = 2,m1-1 xplus(j) = re*rho(2,j)*sqrt(f(2,j,7))*cm4*0.5D+00*hksi(2,j)/rhocon gam(1,j) = rpr if ( xplus(j) > 11.63D+00 ) then gam(1,j) = max(rhocon/re, & xplus(j)/(sigt*(2.5*log(ewall*xplus(j))+pbig))) end if end do do j = 2,m1-1 xplus(j) = re*rho(l1-1,j)*sqrt(f(l1-1,j,7)) & *cm4*0.5*hksi(l1-1,j)/rhocon gam(l1,j) = rpr if ( xplus(j) > 11.63D+00 ) then gam(l1,j) = max(rhocon/re, & xplus(j)/(sigt*(2.5D+00*log(ewall*xplus(j))+pbig))) end if end do end if do j = 2,jcrys f(l1,j,5) = 0.0D+00 end do do j = jcrys+1,m1 f(l1,j,5) = f(l1-1,j,5)-0.1D+00*birad*((f(l1,j,5)*(tw-tf)+tf)**4-tal**4) & *(x(l1,j)-x(l1-1,j))/100000000.0D+00 end do do j = 1,m1 f(1,j,5) = 0.5D+00 + 0.5D+00*yc(2,j) end do do i = 1,l1 f(i,m1,5) = 1.0D+00 end do do i = 1,l1 gam(i,1) = 0.0D+00 end do ! ! Angular momentum. ! else if ( nf == 6 ) then ! ! modelm = 1. based on rectangular de/dx of Sabhapathy and Salcudean ! modelm = 2. based on curvilinear de/dx of Sabhapathy and Salcudean. ! modelm = 3. same as 2 but different numerical treatment for source ! terms based on strong magnetic. 1,2 is good for weak Ha ! 3 is good for moderate. ! modelm = 4. based on strong magnetic. rotation is eliminated. ! modelm = 5. is a test. ! if ( hamag == 0.0D+00 ) then modelm = 0.0D+00 else if ( hamag > 0.0D+00 .and. hamag <= 20.0D+00 ) then modelm = 2 else if ( hamag > 20.0D+00 .and. hamag <= 40.0D+00 ) then modelm = 5 else if ( hamag > 40.0D+00 ) then modelm = 4 end if do i = 2,l1-1 do j = 2,m1-1 con(i,j) = fo(i,j,6)*rho(i,j)*rdtm if ( modelm == 1 ) then con(i,j) = con(i,j)-rho(i,j)*hamag**2/re/r(i,j)* & (f(i+1,j,9)-f(i,j,9))/(x(i+1,j)-x(i,j)) else if ( modelm == 2 ) then call gradnt(heta,hksi,i,j,f(1,1,9),depdx,depdr,xc,yc) if ( j >= 3 ) then con(i,j) = con(i,j)-rho(i,j)*hamag**2/re/r(i,j)*depdx end if else if ( modelm == 3 ) then call gradnt(heta,hksi,i,j,f(1,1,9),depdx,depdr,xc,yc) if ( j >= 3 ) then con(i,j) = con(i,j)-rho(i,j)*hamag**2/re/r(i,j)*depdx & +rho(i,j)*hamag**2/re*f(i,j,6) ap(i,j) = ap(i,j)-rho(i,j)*hamag**2/re end if else if ( modelm == 4 ) then ap(i,j) = ap(i,j)-rho(i,j)*hamag**2/re else if ( modelm == 5 ) then call gradnt(heta,hksi,i,j,f(1,1,9),depdx,depdr,xc,yc) if ( j >= 20 ) then con(i,j) = con(i,j)-rho(i,j)*hamag**2/re/r(i,j)*depdx end if if ( j <= 20 ) then ap(i,j) = ap(i,j)-rho(i,j)*hamag**2/re end if end if ! ! -2/r*dOmega/dr /re ! ar = 2.0D+00/r(i,j)/(y(i,j)-y(i,j-1))/re con(i,j) = con(i,j)+rho(i,j)*ar*f(i,j-1,6) ap(i,j) = ap(i,j)-rho(i,j)*ar end do end do do i = 1,l1 f(i,m1,6) = recb f(i,1,6) = 0.0D+00 end do do j = jcrys+1,m1 gam(l1,j) = 0.0D+00 end do do j = 1,m1 f(1,j,6) = recb*r(2,j)**2 if ( j > jcrys ) then f(l1,j,6) = f(l1-1,j,6) else f(l1,j,6) = rect*r(l1,j)**2 end if end do ! ! Turbulent kinetic energy. ! else if ( nf == 7 ) then do i = 2,l1-1 do j = 2,m1-1 con(i,j) = fo(i,j,7)*rdtm gam(i,j) = rhocon/re+gamt(i,j)/sigk end do end do do i = 2,l1-1 do j = 2,m1-1 call gradnt(heta,hksi,i,j,f(1,1,1),dudx(i,j),dudy(i,j),xc,yc) call gradnt(heta,hksi,i,j,f(1,1,2),dvdx(i,j),dvdy(i,j),xc,yc) prod(i,j) = gamt(i,j)*(2.0D+00*(dudx(i,j)**2+dvdy(i,j)**2) & +(dudy(i,j)+dvdx(i,j))**2) if ( mode == 1 ) then prod(i,j) = prod(i,j)+gamt(i,j)*2.0D+00*(f(i,j,2)/r(i,j))**2 end if if ( lsolve(6) ) then call gradnt(heta,hksi,i,j,f(1,1,6),dwdx(i,j),dwdy(i,j),xc,yc) prod(i,j) = prod(i,j)+gamt(i,j)*(dwdx(i,j)**2 & +(dwdy(i,j)-2.0D+00*f(i,j,6)/r(i,j))**2)/r(i,j)**2 end if con(i,j) = con(i,j)+prod(i,j) ap(i,j) = ap(i,j)-cmu*rho(i,j)**2*abs(f(i,j,7))/(gamt(i,j)+1.0D-05) end do end do ! ! Turbulent dissipation. ! else if ( nf == 8 ) then do i = 2,l1-1 do j = 2,m1-1 con(i,j) = fo(i,j,8)*rdtm gam(i,j) = rhocon/re+gamt(i,j)/sige end do end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = con(i,j)+ce1*cmu*rho(i,j)*abs(f(i,j,7))*prod(i,j) & /((gamt(i,j)+1.0D-05)) ap(i,j) = ap(i,j)-ce2*cmu*rho(i,j)**2*abs(f(i,j,7))/(gamt(i,j)+1.0D-05) end do end do end if return end subroutine gradnt(heta,hksi,i,j,phi,dphidx,dphidy,xc,yc) ! !*****************************************************************************80 ! !! GRADNT estimats the gradient of PHI at a primary node index (I,J). ! ! ! Discussion: ! ! The routine is given: ! ! PHI, the value of a state quantity at the primary nodes; ! HETA, HKSI, the two dimensions of the cell, as measured ! in the X, Y coordinate system; ! I, J, the indices of the primary node at which a gradient ! is desired; ! XC, YC, the X, Y coordinates of the corner nodes that define the ! shape of the cells, and the location of the midside nodes, and ! primary nodes. ! ! GRADNT must estimate ! ! DPHIDX, DPHIDY, an estimate for the gradient (d PHI/d X, d PHI/d Y) ! of the physical quantity PHI at the primary node (I,J). ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, double precision HETA(NI,NJ). ! HETA is the length in the (X,Y) coordinate system, of the cell ! in the ETA direction. This is the actual length of the line ! connecting the mid-side nodes that separate primary node (I,J) ! from primary nodes (I,J-1) and (I,J+1). ! ! Input, double precision HKSI(NI,NJ). ! HKSI is the length in the (X,Y) coordinate system of the cell ! in the KSI direction. This is the actual length of the line ! connecting the mid-side nodes that separate primary node (I,J) ! from primary nodes (I-1,J) and (I+1,J). ! ! Input, integer I, J. ! I and J are the row and column indices of the ! primary node at which the gradient of PHI is desired. ! ! Input, double precision PHI(NI,NJ). ! PHI is the physical quantity whose value is given at each of ! the primary nodes. ! ! Output, double precision DPHIDX, DPHIDY. ! DPHIDX and DPHIDY are the computed approximations to the values ! of the partial derivatives d PHI/d X and d PHI/d Y. ! ! Input, double precision XC(NI,NJ). ! XC contains the X coordinate of nodes which are the "corners" of ! control volumes. ! ! Input, double precision YC(NI,NJ). ! YC contains the Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! double precision detadx double precision detady double precision dpdeta double precision dpdksi double precision dphidx double precision dphidy double precision dxdeta double precision dxdksi double precision dksidx double precision dksidy double precision dydeta double precision dydksi double precision heta(ni,nj) double precision hksi(ni,nj) integer i integer j double precision phi(ni,nj) double precision phie double precision phin double precision phis double precision phiw double precision t1 double precision t2 double precision volume double precision xc(ni,nj) double precision xe double precision xn double precision xs double precision xw double precision yc(ni,nj) double precision ye double precision yn double precision ys double precision yw ! ! We are given the value of some quantity PHI at the primary nodes. ! ! We want to estimate the gradients dPHI/dKSI and dPHI/dETA at ! those primary nodes. ! ! To do so, we must first estimate the value of PHI at the midpoint nodes, ! which lie between the given primary node and each of its four neighbors, ! above, below, to the right and to the left. ! ! Our task is made a little easier because of the fact that, in the ! KSI and ETA coordinate system, primary nodes are separated by ! exactly one unit, (and so are corner nodes, and midside nodes). ! ! This will make it easy to estimate PHI as a function of KSI and ! ETA, and hence to get dPHI/dKSI and dPHI/dETA by using a difference ! of two values. ! ! However, since we really want to compute dPHI/dX and dPHI/dY, ! we must make some adjustments to this problem: ! ! 1) the interpolated value of PHI at the midside node must not be ! the average of the two nearby primary node values, but the WEIGHTED ! average. The weights come from the actual X, Y distances between ! the primary nodes and the midside node in the "real world". ! ! 2) in the X, Y coordinate system, some nodes are separated by ! zero distance. This is just to help us handle boundary conditions, ! but we have to be careful to avoid getting a zero divisor in ! our weights. ! ! 3) we start out by calculating dPHI/dKSI and dPHI/dETA. To get ! dPHI/dX and dPHI/dY, we must compute the jacobian of the ! transformation that maps KSI and ETA to X and Y. Then we ! have to compute the INVERSE of this 2 by 2 matrix, to get ! the terms dKSI/dX, dKSI/dY, dETA/dX and dETA/dY. Then we ! can use the chain rule to calculate terms like: ! ! dPHI/dX = dPHI/dKSI * dKSI/dX + dPHI/dETA * dETA/dX ! ! ! Now let us look at our primary node in the simple (KSI,ETA) ! coordinate system, where: ! ! the (I,J) primary node has (KSI,ETA) coordinates (I,J). ! ! the midside nodes that neighbor primary node (I,J) have ! (KSI,ETA) coordinates as follows: ! ! North: (I, J+1/2) ! South: (I, J-1/2) ! East: (I+1/2, J) ! West: (I-1/2, J) ! ! the corner nodes that form the corners of the cell containing primary ! node (I,J) have (KSI,ETA) coordinates as follows: ! ! Northwest: (I-1/2, J+1/2) ! Northeast: (I+1/2, J+1/2) ! Southwest: (I-1/2, J-1/2) ! Southeast: (I+1/2, J-1/2) ! ! However, if we want to get the X, Y coordinates of a corner node, ! we can't use these (KSI,ETA) subscripts as indices, because FORTRAN ! only accepts whole number indices. So the appropriate indices to ! use to access the XC, YC arrays are: ! ! Northwest: (I, J+1) ! Northeast: (I+1, J+1) ! Southwest: (I, J) ! Southeast: (I+1, J) ! ! In the following diagram, we are interested primary node (I,J). ! ! We show the cell containing (I,J), the neighbor four primary nodes, ! the corner nodes, and the midside nodes. ! ! ! ^ | (I, J+1) | ! | | | ! | | | ! | ----------[I-1/2, J+1/2]--{I, J+1/2}--[I+1/2, J+1/2]---- ! | | | ! E | | ! T (I-1, J) {I-1/2, J] (I, J) {I+1/2, J} (I+1/2, J) ! A | | ! | | | ! | ----------[I-1/2, J-1/2]--{I, J-1/2}--[I+1/2, J-1/2]------ ! | | | ! | | | ! | | (I, J-1) | ! | ! | ! +-------------------KSI ( = I for primary nodes) -------------> ! ! ! Make a weighted average of PHI at the east and at the center, ! to estimate PHIE, the value of PHI at the east midside node: ! if ( hksi(i,j)+hksi(i+1,j) /= 0.0D+00 ) then t1 = hksi(i,j)/(hksi(i,j)+hksi(i+1,j)) t2 = hksi(i+1,j)/(hksi(i,j)+hksi(i+1,j)) else t1 = 0.5D+00 t2 = 0.5D+00 end if phie = t1*phi(i+1,j)+t2*phi(i,j) ! ! Make a weighted average of PHI at the north and at the center, ! to estimate PHIN, the value of PHI at the north midside node. ! if ( heta(i,j)+heta(i,j+1) /= 0.0D+00 ) then t1 = heta(i,j)/(heta(i,j)+heta(i,j+1)) t2 = heta(i,j+1)/(heta(i,j)+heta(i,j+1)) else t1 = 0.5D+00 t2 = 0.5D+00 end if phin = t1*phi(i,j+1)+t2*phi(i,j) ! ! Make a weighted average of PHI at the west and at the center, ! to estimate PHIW, the value of PHI at the west midside node. ! if ( hksi(i,j)+hksi(i-1,j) /= 0.0D+00 ) then t1 = hksi(i,j)/(hksi(i,j)+hksi(i-1,j)) t2 = hksi(i-1,j)/(hksi(i,j)+hksi(i-1,j)) else t1 = 0.5D+00 t2 = 0.5D+00 end if phiw = t1*phi(i-1,j)+t2*phi(i,j) ! ! Make a weighted average of PHI at the south and at the center, ! to estimate PHIS, the value of PHI at the south midside node. ! if ( heta(i,j)+heta(i,j-1) /= 0.0D+00 ) then t1 = heta(i,j)/(heta(i,j)+heta(i,j-1)) t2 = heta(i,j-1)/(heta(i,j)+heta(i,j-1)) else t1 = 0.5D+00 t2 = 0.5D+00 end if phis = t1*phi(i,j-1)+t2*phi(i,j) ! ! Now subtract opposing values, to estimate dPHI/dKSI and dPHI/dETA ! at the primary node. We are assuming that DELTA KSI = DELTA ETA = 1. ! dpdksi = phie-phiw dpdeta = phin-phis ! ! Now we have to convert our results, which are in terms of KSI and ! ETA, into the X and Y coordinate system. ! ! Using the X, Y coordinates of the "corner nodes", compute ! the X, Y coordinates of the midside nodes of the control volume ! interfaces that surround the primary node with coordinates (I,J). ! ! ^ ! | [I-1/2, J+1/2]-----{I, J+1/2}----[I+1/2, J+1/2] ! E | | ! T {I-1/2, J} (I, J) {I+1/2, J} ! A | | ! | [I-1/2, J-1/2]-----{I, J-1/2}----[I+1/2, J] ! | ! +----------------------XSI--------------------------> ! xe = 0.5D+00*(xc(i+1,j+1)+xc(i+1,j)) ye = 0.5D+00*(yc(i+1,j+1)+yc(i+1,j)) xw = 0.5D+00*(xc(i,j+1)+xc(i,j)) yw = 0.5D+00*(yc(i,j+1)+yc(i,j)) xn = 0.5D+00*(xc(i,j+1)+xc(i+1,j+1)) yn = 0.5D+00*(yc(i,j+1)+yc(i+1,j+1)) xs = 0.5D+00*(xc(i,j)+xc(i+1,j)) ys = 0.5D+00*(yc(i,j)+yc(i+1,j)) ! ! From the X and Y coordinates of the midside nodes, ! estimate dX/dKSI, dX/dETA, dY/dKSI and dY/dETA at the primary node, ! again assuming a spacing delta KSI or delta ETA = 1 ! for opposing midside nodes. ! dxdksi = xe-xw dxdeta = xn-xs dydksi = ye-yw dydeta = yn-ys ! ! Compute the determinant of the Jacobian (XSI,ETA)-->(X,Y). ! ! J(XSI,ETA) = ( dX/dKSI dX/dETA) ! ( dY/dKSI dY/dETA) ! volume = dxdksi*dydeta-dxdeta*dydksi ! ! If the determinant is zero, then the control volume has zero ! volume, which should never happen. ! if ( volume == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRADNT - Fatal error!' write ( *, '(a)' ) ' The control volume has zero volume!' stop end if ! ! Now compute the elements of the inverse Jacobian. ! ! J_Inverse(XSI,ETA) = ( dKSI/dX dKSI/dY) ! ( dETA/dX dETA/dY) ! ! using the simple fact that, for a 2 by 2 matrix, ! ! A_Inverse = (a22 -a21) / determinant(A) ! (-a12 a11) ! dksidx = dydeta/volume dksidy = -dxdeta/volume detadx = -dydksi/volume detady = dxdksi/volume ! ! Now we can finally compute dPHI/dX and dPHI/dY using the chain rule. ! dphidx = dpdksi*dksidx + dpdeta*detadx dphidy = dpdeta*detady + dpdksi*dksidy return end subroutine hello ( liwork, lwork ) ! !*****************************************************************************80 ! !! HELLO prints a brief message identifying the program. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer liwork integer lwork ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CRYSTAL_QED:' write ( *, '(a)' ) ' Version of 29 October 1996' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' CRYSTAL_QED seeks a set of parameters that minimize a ' write ( *, '(a)' ) ' cost functional associated with a crystallization ' write ( *, '(a)' ) ' problem.' write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' LIWORK, integer workspace size = ',liwork write ( *, '(a,i6)' ) ' LWORK, real workspace size = ',lwork return end subroutine inidat(ae1,ae2,ak1,ak2,area,b,b1jbl,b2jbl,b3jbl, & birad,bo,cappa,cd,ce1,ce2,cfo,cinc,cmu,cost,cvn,delt,dtm, & epsad,epsil,ewall,f,fcsl,fks,fksl,fma,fn,fnu,fo,fr,frsl, & gam,gamt,grash,hamag,heta,hf,hksi,icost,icrys,inturb,iplot, & ipref,iprint,jcrys,jpref,l0,l1,last,lastt,lblk,lortho, & lsolve,m0,m1,maxbot,mode,nbot,ndt,ni,nj,nk,np,npc,ns, & nsolve,ntimes,orth,pr,ra,rdtm, & re,recb,rect,relax,rho,rhocon,rpr,rueta,ruksi,sige,sigk, & sigma,sigt,smooth,stel,stes,tal,tanca,tanca2,tas,tend,tf,tinit, & title,tnow,tw,vave,vol,vort,x,xbot,xc,xlen,y,ybot,yc,ylen) ! !*****************************************************************************80 ! !! INIDAT sets the initial values of certain data. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, double precision XC(NI,NJ). ! The X coordinate of nodes which are the "corners" of ! control volumes. ! ! Output, double precision YC(NI,NJ). ! The Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer maxbot integer ni integer nj integer nk integer ns ! double precision ae1(ni,nj) double precision ae2(ni,nj) double precision ak1(ni,nj) double precision ak2(ni,nj) double precision area(ni,nj) double precision b double precision b1jbl(nj) double precision b2jbl(nj) double precision b3jbl(ni) double precision birad double precision bo double precision cappa double precision cd double precision ce1 double precision ce2 double precision cfo double precision cinc double precision cl double precision cmu double precision cost double precision cs double precision cvn double precision delt double precision dtm double precision epsad double precision epsil double precision ewall double precision f(ni,nj,ns) double precision fcsl double precision fkl double precision fks double precision fksl double precision fma double precision fn(ni,nj,ns) double precision fnu double precision fo(ni,nj,ns) double precision fr double precision frsl double precision gam(ni,nj) double precision gamt(ni,nj) double precision grash double precision grav double precision hamag double precision heta(ni,nj) double precision hf double precision hksi(ni,nj) integer i integer icost integer icrys integer inturb integer iplot integer ipref integer iprint integer j integer jcrys integer jpref integer k integer l0 integer l1 integer last integer lastt logical lblk(ns) logical lortho logical lsolve(ns) integer m0 integer m1 integer mode integer nbot integer ndt integer np integer npc integer nsolve(ns) integer ntimes(ns) double precision orth double precision pr double precision ra double precision rdtm double precision re double precision recb double precision rect double precision relax(nk) double precision rho(ni,nj) double precision rhocon double precision rhol double precision rhos double precision rpr double precision rueta(ni,nj) double precision ruksi(ni,nj) double precision sige double precision sigk double precision sigma double precision sigt double precision smooth double precision stel double precision stes double precision tal double precision tanca double precision tanca2 double precision tas double precision tend double precision tf double precision tin double precision tinit character ( len = 25 ) title(ns) double precision tnow double precision tw double precision vave double precision vol(ni,nj) double precision vort(ni,nj) double precision x(ni,nj) double precision xbot(maxbot) double precision xc(ni,nj) double precision xlen double precision y(ni,nj) double precision ybot(maxbot) double precision yc(ni,nj) double precision ylen ! ! Set quantities whose values, initial values, or dummy values ! do not depend on other quantities. ! do i = 1,ni do j = 1,nj ae1(i,j) = 0.0D+00 ae2(i,j) = 0.0D+00 ak1(i,j) = 0.0D+00 ak2(i,j) = 0.0D+00 end do end do do i = 1,ni do j = 1,nj area(i,j) = 0.0D+00 end do end do b = 0.1778D+00 do i = 1,nj b1jbl(i) = 0.0D+00 b2jbl(i) = 0.0D+00 end do do i = 1,ni b3jbl(i) = 0.0D+00 end do cappa = 0.4187D+00 cd = 1.0D+00 ce1 = 1.44D+00 ce2 = 1.92D+00 cinc = 0.0D+00 cl = 1000.0D+00 cmu = 0.09D+00 cost = 0.0D+00 cs = 1000.0D+00 cvn = 200000.0D+00 dtm = 100.0D+00 epsad = 0.0001D+00 epsil = 0.0001D+00 ewall = 9.793D+00 do i = 1,ni do j = 1,nj do k = 1,ns f(i,j,k) = 0.0D+00 end do end do end do fkl = 64.0D+00 fks = 22.0D+00 fma = -1000.0D+00 do i = 1,ni do j = 1,nj do k = 1,ns fn(i,j,k) = 0.0D+00 end do end do end do fnu = 3.0D-07 do i = 1,ni do j = 1,nj do k = 1,ns fo(i,j,k) = 0.0D+00 end do end do end do fr = 1.0D-10 do i = 1,ni do j = 1,nj gam(i,j) = 0.0D+00 end do end do grash = 10000000.0D+00 grav = 9.81D+00 hamag = 0.0D+00 heta(1:ni,1:nj) = 0.0D+00 hf = 1800000.0D+00 do i = 1,ni do j = 1,nj hksi(i,j) = 0.0D+00 end do end do ! ! ICOST = 1, U**2+V**2 ! 2, (T-TF)**2 ! 3, (VAVE-SQRT(U**2+V**2)) ! icost = 1 icrys = 42 inturb = 0 ! ! IPLOT = 0, do not create a plot file. ! 1, create plot file. ! iplot = 1 ipref = 11 ! ! IPRINT = 0, don't print very much. ! 1, print intermediate information. ! iprint = 0 jcrys = 22 jpref = 30 l0 = 62 l1 = 62 ! ! Interstep iteration number. ! Use LAST = 40 for more accurate results. ! Use LAST = 10 for quicker results. ! last = 40 ! ! Set the number of timesteps. ! ! For tests, set LASTT = 2. ! For a real run, try LASTT = 200. ! lastt = 2 do i = 1,4 lblk(i) = .false. end do do i = 5,10 lblk(i) = .true. end do lortho = .false. do i = 1,ns lsolve(i) = .false. end do m0 = 42 m1 = 42 mode = 1 ndt = 0 np = 3 npc = 4 ! ! Number of linear iterations. ! nsolve(1) = 3 nsolve(2) = 3 nsolve(3) = 1 nsolve(4) = 1 nsolve(5) = 3 nsolve(6) = 3 nsolve(7) = 3 nsolve(8) = 3 nsolve(9) = 1 nsolve(10) = 1 ! ! Number of nonlinear iterations. ! ntimes(1) = 5 ntimes(2) = 5 ntimes(3) = 1 ntimes(4) = 1 ntimes(5) = 3 ntimes(6) = 3 ntimes(7) = 3 ntimes(8) = 3 ntimes(9) = 1 ntimes(10) = 1 orth = 50000.0D+00 pr = 0.015D+00 rdtm = 0.0D+00 re = 1.0D+00 recb = 0.0D+00 rect = 0.0D+00 relax(1) = 0.3D+00 relax(2) = 0.3D+00 relax(3) = 0.8D+00 relax(4) = 1.0D+00 relax(5) = 0.7D+00 relax(6) = 0.2D+00 relax(7) = 0.5D+00 relax(8) = 0.5D+00 relax(9) = 1.0D+00 relax(10) = 1.0D+00 relax(11) = 1.0D+00 relax(12) = 0.6D+00 relax(13) = 1.0D+00 relax(14) = 1.0D+00 rhocon = 1.0D+00 rhol = 2490.0D+00 rhos = 2490.0D+00 do i = 1,ni do j = 1,nj rueta(i,j) = 0.0D+00 ruksi(i,j) = 0.0D+00 end do end do sige = 1.3D+00 sigk = 1.0D+00 sigma = 0.72D+00 sigt = 0.9D+00 smooth = 0.0001D+00 tal = 1523.0D+00 tanca = 1.0D+00 tanca2 = 1.0D+00 tas = 1523.0D+00 tf = 1683.0D+00 tin = 1523.0D+00 tinit = 0.0D+00 title(1) = 'U velocity' title(2) = 'V velocity' title(3) = 'Pressure' title(4) = 'Corrected pressure' title(5) = 'Temperature' title(6) = 'Rotational velocity' title(7) = 'Turbulent dissipation' title(8) = 'Turbulent energy' title(9) = 'Magnetic stream function' title(10) = 'Stream function' tw = 1713.0D+00 vave = 1850.0D+00 do i = 1,ni do j = 1,nj vol(i,j) = 0.0D+00 end do end do do i = 1,ni do j = 1,nj vort(i,j) = 0.0D+00 end do end do do i = 1,ni do j = 1,nj x(i,j) = 0.0D+00 end do end do do i = 1,ni do j = 1,nj xc(i,j) = 0.0D+00 end do end do xlen = 1.2 do i = 1,ni do j = 1,nj y(i,j) = 0.0D+00 yc(i,j) = 0.0D+00 end do end do ylen = 1.0 ! ! Set things that depend on things. ! birad = 0.25D+00 * 5.67D+00 * 1.5D+00 * 0.0254D+00 /fks/(tw-tf)/re bo = (rhol*grav*b**2)/sigma cfo = fnu/(b*b) if ( inturb == 1 ) then do i = 1,ni do j = 1,nj f(i,j,7) = 0.005D+00 end do end do end if fcsl = cs/cl fksl = fks/fkl frsl = rhos/rhol ra = grash*pr rho(1:ni,1:nj) = rhocon rpr = rhocon/(pr*re) stel = cl*(tw-tf)/hf stes = cs*(tf-tin)/hf tend = tinit+dtm*lastt tnow = tinit ! ! Set things that depend on things that depend on things. ! delt = cfo*dtm if ( inturb == 1 ) then do i = 1,ni do j = 1,nj f(i,j,8) = f(i,j,7)**1.5 / 0.006D+00 gamt(i,j) = 0.006D+00 *cmu*rho(i,j)*f(i,j,7)**0.5 end do end do end if ! ! Set the crucible shape. ! xbot(1) = 0.0D+00 ybot(1) = 0.0D+00 xbot(2) = 0.015D+00 ybot(2) = 0.3D+00 xbot(3) = 0.03D+00 ybot(3) = 0.5D+00 xbot(4) = 0.046D+00 ybot(4) = 0.6D+00 xbot(5) = 0.07D+00 ybot(5) = 0.7D+00 xbot(6) = 0.10D+00 ybot(6) = 0.8D+00 xbot(7) = 0.14D+00 ybot(7) = 0.9D+00 xbot(8) = 0.18D+00 ybot(8) = 0.95D+00 xbot(9) = 0.25D+00 ybot(9) = 1.0D+00 nbot = 9 return end subroutine inigrd ( icrys, jcrys, l0, m0, nbot, xbot, xc, xlen, ybot, & yc, ylen ) ! !*****************************************************************************80 ! !! INIGRD makes an initial assignment of the grid points XC, YC. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, double precision XC(NI,NJ). ! The X coordinate of nodes which are the "corners" of ! control volumes. ! ! Output, double precision YC(NI,NJ). ! The Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer nbot integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! double precision free integer i integer icrys integer j integer jcrys integer l0 integer licrys integer ll1 integer llst1 integer m0 double precision t1 double precision t2 double precision wave double precision xbot(nbot) double precision xc(ni,nj) double precision xlen double precision ybot(nbot) double precision yc(ni,nj) double precision ylen ! ll1 = icrys/2+1 llst1 = jcrys/2+1 licrys = (jcrys+m0)/2+1 free = 0.5D+00*xlen ! ! Set the Y coordinates ! do i = 2, l0 if ( i < ll1 ) then t2 = free*(0.5D+00*(dble(i-2)/dble(ll1-2))**1.5) else if ( i < icrys ) then t2 = free*(1.0D+00-0.5D+00*(dble(icrys-i)/dble(icrys-ll1))**1.5) else t2 = free+(xlen-free)*(dble(i-icrys)/dble(l0-icrys))**1.2 end if do j = 2, m0 wave = 0.4D+00 + 0.02D+00 * sin ( t2 * 20.0D+00 * 3.14159D+00 ) if ( j < llst1 ) then t1 = wave/2.0D+00*(dble(j-2)/dble(llst1-2))**1.5 else if ( j < jcrys ) then t1 = wave-wave/2.0D+00*(dble(jcrys-j)/dble(jcrys-llst1))**1.5 else if ( j < licrys ) then t1 = 0.4D+00 + 0.3D+00*(dble(j-jcrys)/dble(licrys-jcrys))**1.5 else t1 = 1.0D+00 - 0.3D+00*(dble(m0-j)/dble(m0-licrys))**1.5 end if yc(i,j) = t1*ylen end do end do ! ! Set the XC coordinates along the bottom boundary. ! xc(2,2) = xbot(1) do j = 3, m0-1 call cubic ( nbot, yc(2,j), ybot, xc(2,j), xbot ) end do xc(2,m0) = xbot(nbot) ! ! Now assign the XC coordinates of the other nodes. ! Here, we explicitly assume that the free surface and crystal ! boundaries occur at a coordinate value of 0.6. ! do j = 2, m0 do i = 3, l0 if ( i < ll1 ) then t1 = 0.5D+00 * (dble(i-2)/dble(ll1-2))**1.5 xc(i,j) = (1.0D+00 - t1)*xc(2,j)+t1*free else if ( i <= icrys ) then t1 = 1.0D+00 - 0.5D+00 * (dble(icrys-i)/dble(icrys-ll1))**1.5 xc(i,j) = ( 1.0D+00 - t1)*xc(2,j)+t1*free else t1 = (dble(i-icrys)/dble(l0-icrys))**1.2 if ( j < jcrys ) then xc(i,j) = xc(icrys,j)+t1*((xlen-free)+0.5D+00*yc(l0,jcrys)**2 & -0.5D+00 * yc(l0,j)**2) else xc(i,j) = xc(icrys,j)+t1*(xlen-free) end if end if end do end do ! ! Copy values to dummy nodes with I = 1 or J=1. ! xc(1,1) = xc(2,2) yc(1,1) = yc(2,2) do j = 2, m0 xc(1,j) = xc(2,j) yc(1,j) = yc(2,j) end do do i = 2, l0 xc(i,1) = xc(i,2) yc(i,1) = yc(i,2) end do return end subroutine movgrd ( b1jbl, b2jbl, bo, delt, fksl, fr, frsl, icrys, & iprint, jcrys, l0, m0, p, pr, re, stel, tanca, tanca2, xc, xlen, yc ) ! !*****************************************************************************80 ! !! MOVGRD calculates the new position of the interface and free surface. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, double precision XC(NI,NJ). ! The X coordinate of nodes which are the "corners" of ! control volumes. ! ! Input/output, double precision YC(NI,NJ). ! The Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! double precision ac double precision as double precision b1jbl(nj) double precision b2jbl(nj) double precision bo double precision coeff double precision coff double precision d2hdy2(nj) double precision delt double precision dhdy(nj) double precision dhdym(nj) double precision dhdyp(nj) double precision dlen double precision dx1 double precision dxm1 double precision dy1 double precision dym1 double precision fksl double precision flow(ni) double precision flomax double precision flomin double precision fr double precision frsl double precision high(ni) double precision hilev double precision himax double precision himin integer i integer icrys integer interl integer iprint integer j integer jcrys integer k integer l0 integer m0 double precision omega double precision omega2 double precision p(ni,nj) double precision pcsurf double precision pr double precision pres(ni) double precision pullv double precision re double precision rr(nj) double precision stel double precision tanca double precision tanca2 double precision xb(ni) double precision xc(ni,nj) double precision xlen double precision yc(ni,nj) double precision zb1(ni) ! do j = 2,m0 zb1(j) = xc(icrys,j) end do do k = 1,5 do j = jcrys+1,m0-1 dhdy(j) = (xc(icrys,j+1)-xc(icrys,j))/(yc(icrys,j+1)-yc(icrys,j)) dhdym(j) = (xc(icrys,j)-xc(icrys,j-1))/(yc(icrys,j)-yc(icrys,j-1)) dhdyp(j) = (xc(icrys,j+1)-xc(icrys,j))/(yc(icrys,j+1)-yc(icrys,j)) end do do j = jcrys+1,m0-1 d2hdy2(j) = (dhdyp(j)-dhdym(j))/(yc(icrys,j+1)-yc(icrys,j)) rr(j) = d2hdy2(j)/(1.0D+00 + (dhdy(j)*dhdy(j)))**1.5 & +dhdy(j)/yc(icrys,j)/(1.0D+00 + (dhdy(j)*dhdy(j)))**0.5 end do pres(jcrys+1) = 0.0D+00 do j = jcrys+1,m0-2 pres(j+1) = pres(j)+(p(icrys-1,j+1)-p(icrys-1,j))/yc(icrys-1,j+1) end do himax = 0.0D+00 himin = 0.0D+00 interl = (jcrys+m0)/2 do j = jcrys+1,m0-1 pcsurf = pres(j)-pres(interl) high(j) = fr*pcsurf+rr(j)/bo end do hilev = high(interl) omega = 0.0D+00 do j = jcrys+1,m0-1 high(j) = omega*((high(j)-hilev)-(xc(icrys,j)-xc(icrys,interl))) himax = max(himax,high(j)) himin = min(himin,high(j)) xc(icrys,j) = xc(icrys,j)+high(j) end do if ( omega /= 0.0 ) then xc(icrys,jcrys) = xc(icrys,jcrys+1) & +tanca*(yc(icrys,jcrys+1)-yc(icrys,jcrys)) xc(icrys,m0) = xc(icrys,m0-1)+tanca2*(yc(icrys,m0)-yc(icrys,m0-1)) end if flomax = 0.0D+00 flomin = 0.0D+00 if ( iprint > 0 ) then if ( k == 5 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MOVGRD:' write ( *, '(a,g14.6)' ) ' Himin = ',himin write ( *, '(a,g14.6)' ) ' Himax = ',himax write ( *, '(a,g14.6)' ) ' Omega = ',omega write ( *, '(a)' ) ' ' write(*,'(7f10.4)')high(jcrys+1),high(jcrys+2),high(m0-10), & high(m0-8),high(m0-4),high(m0-2),high(m0-1) write(*,'(7f9.4)')xc(icrys,jcrys),xc(icrys,jcrys+1), & xc(icrys,jcrys+4),xc(icrys,m0-8),xc(icrys,m0-4), & xc(icrys,m0-1),xc(icrys,m0) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MOVGRD:' write ( *, '(a)' ) ' Free surface movement' end if end if end do omega2 = 0.0D+00 do j = 2,jcrys-1 dxm1 = xc(icrys,j+1)-xc(icrys,j) dym1 = yc(icrys,j+1)-yc(icrys,j) dlen = sqrt(dxm1*dxm1+dym1*dym1) as = -dxm1/dlen ac = dym1/dlen dx1 = delt*frsl*(b1jbl(j)-fksl*b2jbl(j))*stel/(re*pr)*ac dy1 = delt*frsl*(b1jbl(j)-fksl*b2jbl(j))*stel/(re*pr)*as if ( j == jcrys-1 ) then flow(j) = dx1 else flow(j) = dx1+dy1**2/dx1 end if end do ! ! What is PULLV? ! pullv = flow(jcrys-1) do j = 3,jcrys-1 flomax = max(flomax,(flow(j)-pullv)) flomin = min(flomin,(flow(j)-pullv)) end do omega2 = min(omega2, 0.003D+00/(flomax+1.0D-06)) omega2 = min(omega2, 0.003D+00/(-flomin+1.0D-06)) do j = 3,jcrys-1 flow(j) = omega2*(flow(j)-pullv)+(xc(icrys,jcrys)-zb1(jcrys)) xb(j) = xc(icrys,j)+flow(j) xc(icrys,j) = xb(j) end do xc(icrys,2) = xc(icrys,3) if ( iprint > 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MOVGRD:' write ( *, '(a,g14.6)' ) ' Flomin = ',flomin write ( *, '(a,g14.6)' ) ' Flomax = ',flomax write ( *, '(a,g14.6)' ) ' Omega = ',omega2 write ( *, '(a)' ) ' ' write(*,'(2x,7f9.4)')xb(jcrys-1),xb(jcrys-2),xb(8),xb(6),xb(4),xb(3),xb(2) write(*,'(2x,7f9.4)')b1jbl(jcrys-1),-fksl*b2jbl(jcrys-1), & b1jbl(jcrys-2),-fksl*b2jbl(jcrys-2),b1jbl(jcrys-5),-fksl*b2jbl(jcrys-5) write(*,'(2x,7f9.4)') b1jbl(jcrys-1)-fksl*b2jbl(jcrys-1), & b1jbl(jcrys-2)-fksl*b2jbl(jcrys-2),b1jbl(jcrys-5)-fksl*b2jbl(jcrys-5), & b1jbl(jcrys-10)-fksl*b2jbl(jcrys-10),b1jbl(jcrys-15)-fksl*b2jbl(jcrys-15), & b1jbl(jcrys-17)-fksl*b2jbl(jcrys-17),b1jbl(jcrys-19)-fksl*b2jbl(jcrys-19) write(*,'(2x,7f9.4)')flow(21),flow(19),flow(17),flow(15), & flow(13),flow(11),flow(10) write(*,'(2x,7f9.4)')flow(9),flow(8),flow(7),flow(6),flow(5),flow(4),flow(3) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MOVGRD:' write ( *, '(a)' ) ' Solid-liquid interface movement' write ( *, '(a)' ) ' ' end if do j = 2,m0 coeff = (xc(2,j)-xc(icrys,j))/(xc(2,j)-zb1(j)) coff = (xlen-xc(icrys,j))/(xlen-zb1(j)) do i = 2,icrys-1 xc(i,j) = xc(2,j)-(xc(2,j)-xc(i,j))*coeff end do do i = icrys+1,l0 xc(i,j) = xlen-(xlen-xc(i,j))*coff if ( j <= jcrys .and. i == l0 ) then xc(l0,j) = xlen+0.5D+00*yc(l0,jcrys)**2 - 0.5D+00*yc(l0,j)**2 end if end do end do return end subroutine output ( cfo, iter, izone, res, smax, ssum, t, tnow, u, w ) ! !*****************************************************************************80 ! !! OUTPUT prints some sample data from the ongoing solution. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: ns = 10 ! double precision cfo integer iter integer izone double precision res(ns) double precision smax double precision ssum double precision t(ni,nj) double precision tnow double precision u(ni,nj) double precision w(ni,nj) ! ! Solid zone output. ! if ( izone == 1 ) then if ( iter == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'OUTPUT:' write ( *, '(a,g14.6)' ) ' Current time = ',tnow write ( *, '(a,g14.6,a)' ) ' = ',tnow*cfo,' seconds.' write ( *, '(a)' ) ' Solid zone results.' write ( *, '(a)' ) ' Iter Res(5) T(50,21) T(50,22)' write ( *, '(a)' ) ' ' end if write(*,'(i4,2x,10g12.4)')iter,res(5),t(50,21),t(50,22) ! ! Liquid zone output. ! else if ( izone == 2 ) then if ( iter == 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'OUTPUT:' write ( *, '(a)' ) ' Liquid zone results.' write ( *, '(a)' ) ' SMAX is the maximum local mass imbalance.' write ( *, '(a)' ) ' SSUM is the total mass imbalance.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Iter SMAX SSUM U(5,5) W(5,5) T(5,5)' write ( *, '(a)' ) ' ' end if write(*,'(i4,2x,10g12.4)')iter,smax,ssum,u(5,5),w(5,5),t(5,5) end if return end subroutine pmod ( ipref, jcrys, jpref, l1, m1, p ) ! !*****************************************************************************80 ! !! PMOD interpolates the pressure at boundary corners. ! ! ! Discussion: ! ! PMOD carries out some simple operations to extend the value of the ! pressure to primary nodes at the corners, and to normalize the ! pressure by subtracting off its value at the reference point. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! integer i integer ipref integer j integer jcrys integer jpref integer l1 integer m1 double precision p(ni,nj) double precision pref ! ! Extrapolate to get pressures on the boundary corners. ! p(1,1) = p(2,1)+p(1,2)-p(2,2) p(l1,1) = p(l1-1,1)+p(l1,2)-p(l1-1,2) p(1,jcrys) = p(2,jcrys)+p(1,jcrys-1)-p(2,jcrys-1) p(l1,jcrys) = p(l1-1,jcrys)+p(l1,jcrys-1)-p(l1-1,jcrys-1) ! ! Subtract off the reference pressure. ! pref = p(ipref,jpref) do j = 1,m1 do i = 1,l1 p(i,j) = p(i,j)-pref end do end do return end subroutine prdat(b,birad,bo,cappa,cfo,cvn,delt,dtm,fcsl,fksl, & fma,fnu,fr,frsl,grash,icost,icrys,inturb,iprint, & jcrys,l0,last,lastt,m0,mode,nbot,ns,nsolve,ntimes,orth, & pr,ra,rdtm,recb,rect,rhocon,smooth,stel,stes,tanca,tanca2,tend, & tf,tinit,title,tw,vave,xlen,ylen) ! !*****************************************************************************80 ! !! PRDAT prints out the initial values of certain data. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer ns ! double precision b double precision birad double precision bo double precision cappa double precision cfo double precision cvn double precision delt double precision dtm double precision fcsl double precision fksl double precision fma double precision fnu double precision fr double precision frsl double precision grash integer i integer icost integer icrys integer inturb integer iprint integer jcrys integer l0 integer last integer lastt integer m0 integer mode integer nbot integer nsolve(ns) integer ntimes(ns) double precision orth double precision pr double precision ra double precision rdtm double precision recb double precision rect double precision rhocon double precision smooth double precision stel double precision stes double precision tanca double precision tanca2 double precision tend double precision tf double precision tinit character ( len = 25 ) title(ns) double precision tw double precision vave double precision xlen double precision ylen ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PRDAT:' write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' B, Crucible radius = ',b write ( *, '(a,g14.6)' ) ' BIRAD, thermal coefficient = ',birad write ( *, '(a,g14.6)' ) ' BO, the Bond number = ',bo write ( *, '(a,g14.6)' ) ' CAPPA, wall function constant = ',cappa write ( *, '(a,g14.6)' ) ' CFO, time scale = ',cfo write ( *, '(a,g14.6)' ) ' CVN, ADAPT cell volume weight = ',cvn write ( *, '(a,g14.6)' ) ' DELT, time step in seconds = ',delt write ( *, '(a,g14.6)' ) ' DTM, dimensionless time step = ',dtm write ( *, '(a,g14.6)' ) ' FCSL = CS/CL = ',fcsl write ( *, '(a,g14.6)' ) ' FKSL = FKS/FKL = ',fksl write ( *, '(a,g14.6)' ) ' FMA, Marangoni number FMA = ',fma write ( *, '(a,g14.6)' ) ' FR, Froude number = ',fr write ( *, '(a,g14.6)' ) ' FRSL = RHOS/RHOL = ',frsl write ( *, '(a,g14.6)' ) ' GRASH, Grashof number = ',grash write ( *, '(a,i6)' ) ' ICOST, cost functional = ',icost write ( *, '(a,i6)' ) ' ICRYS, maximum I of crystal = ',icrys write ( *, '(a,i6)' ) ' INTURB, turbulence option = ',inturb write ( *, '(a,i6)' ) ' IPRINT, printing option = ',iprint write ( *, '(a,i6)' ) ' JCRYS, maximum J of crystal = ',jcrys write ( *, '(a,i6)' ) ' L0, number of I nodes = ',l0 write ( *, '(a,i6)' ) ' LAST, number of zone iterations on ' write ( *, '(a,i6)' ) ' each time step = ',last write ( *, '(a,i6)' ) ' LASTT, number of time steps = ',lastt write ( *, '(a,i6)' ) ' M0, number of J nodes = ',m0 write ( *, '(a,i6)' ) ' MODE, 0 cartesian, 1 axisymmetric = ',mode write ( *, '(a,i6)' ) ' NBOT, number of boundary points = ',nbot write ( *, '(a,g14.6)' ) ' ORTH, ADAPT orthogonality weight = ',orth write ( *, '(a,g14.6)' ) ' PR, Prandtl number PR = ',pr write ( *, '(a,g14.6)' ) ' RA, Rayleigh number = ',ra write ( *, '(a,g14.6)' ) ' RBD = FMA/RA = ',fma/ra write ( *, '(a,g14.6)' ) ' RDTM, = 1/DTM or 0 = ',rdtm write ( *, '(a,g14.6)' ) ' RECB, crucible Reynolds number = ',recb write ( *, '(a,g14.6)' ) ' RECT, crystal Reynolds number = ',rect write ( *, '(a,g14.6)' ) ' RHOCON, density constant = ',rhocon write ( *, '(a,g14.6)' ) ' SMOOTH, ADAPT smoothness weight = ',smooth write ( *, '(a,g14.6)' ) ' STEL, liquid Stefan number STEL = ',stel write ( *, '(a,g14.6)' ) ' STES, solid Stefan number STES = ',stes write ( *, '(a,g14.6)' ) ' TANCA = ',tanca write ( *, '(a,g14.6)' ) ' TANCA2 = ',tanca2 write ( *, '(a,g14.6)' ) ' TEND, dimensionless end time = ',tend write ( *, '(a,g14.6)' ) ' TF, crystal melting temperature = ',tf write ( *, '(a,g14.6)' ) ' TINIT, dimensionless start time = ',tinit write ( *, '(a,g14.6)' ) ' TW, the wall temperature = ',tw write ( *, '(a,g14.6)' ) ' VAVE, desired average velocity = ',vave write ( *, '(a,g14.6)' ) ' WEBER = FR*BO = ',fr*bo write ( *, '(a,g14.6)' ) ' XLEN = problem region length = ',xlen write ( *, '(a,g14.6)' ) ' YLEN = problem region height = ',ylen write ( *, '(a,g14.6)' ) ' Characteristic velocity FNU/B = ',fnu/b write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Variable Nonlinear Linear' write ( *, '(a)' ) ' Iterations Iterations' write ( *, '(a)' ) ' ' do i = 1,ns write ( *, '(a,i6,i6)' ) title(i),ntimes(i),nsolve(i) end do return end subroutine resid ( aim, aip, ajm, ajp, ap, con, f, l1, m1, nf ) ! !*****************************************************************************80 ! !! RESID computes the residual of the linear equations. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: ns = 10 ! double precision aim(ni,nj) double precision aip(ni,nj) double precision ajm(ni,nj) double precision ajp(ni,nj) double precision ap(ni,nj) double precision con(ni,nj) double precision conmax double precision diamax double precision f(ni,nj,ns) integer i integer imax integer iprint integer j integer jmax integer l1 integer m1 integer nf double precision offmax double precision ratio double precision ratmin double precision res double precision resmax double precision sum double precision xmax ! iprint = 0 ! conmax = 0.0D+00 diamax = 0.0D+00 offmax = 0.0D+00 ratmin = 100000.0D+00 resmax = 0.0D+00 xmax = 0.0D+00 imax = 0 jmax = 0 do i = 2,l1-1 do j = 2,m1-1 res = con(i,j)-ap(i,j)*f(i,j,nf)+aim(i,j)*f(i-1,j,nf) & +aip(i,j)*f(i+1,j,nf) & +ajm(i,j)*f(i,j-1,nf)+ajp(i,j)*f(i,j+1,nf) if ( abs(res) >= resmax ) then resmax = abs(res) imax = i jmax = j end if if ( abs(con(i,j)) >= conmax ) then conmax = abs(con(i,j)) end if if ( abs(f(i,j,nf)) >= xmax ) then xmax = abs(f(i,j,nf)) end if if ( abs(ap(i,j)) >= diamax ) then diamax = abs(ap(i,j)) end if sum = abs(aim(i,j))+abs(aip(i,j))+abs(ajm(i,j))+abs(ajp(i,j)) if ( sum >= offmax ) then offmax = sum end if if ( sum /= 0.0 ) then ratio = abs(ap(i,j))/sum if ( ratio < ratmin ) then ratmin = ratio end if end if end do end do if ( iprint == 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RESID' write ( *, '(a,i6)' ) ' Maxixum residual for variable ',nf write ( *, '(a,g14.6)' ) ' is ',resmax write ( *, '(a,2i6)' ) ' at I,J = ',imax,jmax i = imax j = jmax write ( *, '(a)' ) ' ' write ( *, '(a,2g14.6)' ) ' A(I,J)*X(I,J) = ',ap(i,j),f(i,j,nf) write ( *, '(a,2g14.6)' ) ' A(I-1,J)*X(I-1,J) = ',aim(i,j),f(i-1,j,nf) write ( *, '(a,2g14.6)' ) ' A(I+1,J)*X(I+1,J) = ',aip(i,j),f(i+1,j,nf) write ( *, '(a,2g14.6)' ) ' A(I,J-1)*X(I,J-1) = ',ajm(i,j),f(i,j-1,nf) write ( *, '(a,2g14.6)' ) ' A(I,J+1)*X(I,J+1) = ',ajp(i,j),f(i,j+1,nf) write ( *, '(a,g14.6)' ) ' CON(I,J) = ',con(i,j) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' Norm of variable is ',xmax write ( *, '(a,g14.6)' ) ' Norm of RHS is ',conmax write ( *, '(a,g14.6)' ) ' Norm of diagonal is ',diamax write ( *, '(a,g14.6)' ) ' Norm of off-diag is ',offmax write ( *, '(a,g14.6)' ) ' Min diag/off-diag ',ratmin end if return end subroutine rswrit(b1jbl,b2jbl,b3jbl,cost,e,gamt,icrys,jcrys, & l0,m0,nbot,p,pc,psi,rueta,ruksi,t,te,tk,tnow,u,v,vort,w,x, & xbot,xc,y,ybot,yc) ! !*****************************************************************************80 ! !! RSWRIT writes out information which can be used for restarts or plots. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real COST, the value of the cost functional ! associated with this solution. ! ! Input, real E(NI,NJ), the magnetic stream function. ! ! Input, real GAMT(NI,NJ), the diffusion coefficient. ! ! Input, integer ICRYS, JCRYS, the I and J coordinates of ! the lower left corner of the crystal. ! ! Input, integer L0, the number of rows of data. ! ! Input, integer M0, the number of columns of data. ! ! Input, integer NBOT, the number of crucible nodes. ! ! Input, real P(NI,NJ), the pressure. ! ! Input, real PC(NI,NJ), the corrected pressure. ! ! Input, real PSI(NI,NJ), the stream function. ! ! Input, real RUETA(NI,NJ), RUKSI(NI,NJ), the ! the momentum in the ETA and KSI directions. ! ! Input, real T(NI,NJ), the temperature. ! ! Input, real TE(NI,NJ), the turbulent epsilon. ! ! Input, real TK(NI,NJ), the turbulent K. ! ! Output, real TNOW, the current time. ! ! Input, real U(NI,NJ), the horizontal velocity. ! ! Input, real V(NI,NJ), the vertical velocity. ! ! Input, real VMAG(NI,NJ), the velocity magnitude. ! ! Input, real VORT(NI,NJ), the vorticity. ! ! Input, real W(NI,NJ), the axial velocity. ! ! Input, real X(NI,NJ), the X coordinates of the primary ! nodes. ! ! Input, real XBOT(NBOT), the X coordinates of the crucible ! bottom nodes. ! ! Input, real XC(NI,NJ), the X coordinates of the corner ! nodes. ! ! Input, real Y(NI,NJ), the Y coordinates of the primary ! nodes. ! ! Input, real YBOT(NBOT), the Y coordinates of the crucible ! bottom nodes. ! ! Input, real YC(NI,NJ), the Y coordinates of the corner ! nodes. ! implicit none ! integer nbot integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! double precision b1jbl(nj) double precision b2jbl(nj) double precision b3jbl(ni) double precision cost logical, parameter :: debug = .false. double precision e(ni,nj) double precision gamt(ni,nj) integer i integer icrys integer j integer jcrys integer l0 integer m0 double precision p(ni,nj) double precision pc(ni,nj) double precision psi(ni,nj) double precision rueta(ni,nj) double precision ruksi(ni,nj) double precision t(ni,nj) double precision te(ni,nj) double precision tk(ni,nj) double precision tnow double precision u(ni,nj) double precision v(ni,nj) double precision vort(ni,nj) double precision w(ni,nj) double precision x(ni,nj) double precision xbot(nbot) double precision xc(ni,nj) double precision y(ni,nj) double precision ybot(nbot) double precision yc(ni,nj) ! integer imax integer imin integer jmax integer jmin double precision vmax double precision vmin ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RSWRIT:' write ( *, '(a)' ) ' Writing restart information.' write ( *, '(a)' ) ' ' ! ! Temporary check that VORT contains legitimate data. ! if ( debug ) then vmax = vort(1,1) imax = 1 jmax = 1 vmin = vort(1,1) imin = 1 jmin = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I, J, VORT(I,J)' write ( *, '(a)' ) ' ' do i = 1,l0 do j = 1,m0 if ( vort(i,j) > vmax ) then vmax = vort(i,j) imax = i jmax = j end if if ( vort(i,j) < vmin ) then vmin = vort(i,j) imin = i jmin = j end if write(*,'(2i6,g14.5)')i,j,vort(i,j) end do end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RSWRIT - Temporary note:' write ( *, '(a,g14.6)' ) ' Minimum vorticity is ',vmin write ( *, '(a,i6,a,i6)' ) ' IMIN = ',imin,' JMIN=',jmin write ( *, '(a,g14.6)' ) ' Maximum vorticity is ',vmax write ( *, '(a,i6,a,i6)' ) ' IMAX = ',imax,' JMAX=',jmax write ( *, '(a,g14.6)' ) ' VMIN = ', vmin write ( *, '(a,g14.6)' ) ' VMAX = ', vmax end if ! ! Delete any old copy of the restart file. ! open(unit = 10,file='rswrit.txt',status='old',err=10) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RSWRIT - Note:' write ( *, '(a)' ) ' Deleting a previous copy of the restart file.' close(unit = 10,status='delete') 10 continue open(unit = 10,file='rswrit.txt',status='unknown') write(10,*)cost write(10,*)l0 write(10,*)jcrys write(10,*)icrys write(10,*)m0 write(10,*)nbot write(10,*)tnow ! ! NEW NEW NEW ! write(10,'(6g14.5)')(b1jbl(i),i = 1,nj) write(10,'(6g14.5)')(b2jbl(i),i = 1,nj) write(10,'(6g14.5)')(b3jbl(i),i = 1,ni) write(10,'(6g14.5)')((e(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((gamt(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((p(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((pc(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((psi(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((rueta(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((ruksi(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((t(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((te(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((tk(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((u(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((v(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((vort(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((w(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((x(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')(xbot(i),i = 1,nbot) write(10,'(6g14.5)')((xc(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((y(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')(ybot(i),i = 1,nbot) write(10,'(6g14.5)')((yc(i,j),i = 1,l0),j=1,m0) close(unit = 10) return end subroutine setcst ( area, cinc, gam, icost, icrys, jcrys, l0, m0, ni, & nj, t, tf, u, v, vave, vort, xc, yc ) ! !*****************************************************************************80 ! !! SETCST computes a portion on the cost functional. ! ! ! Discussion: ! ! The cost functional is currently the integral over time and ! space of the norm of the fluid velocity. ! ! Currently, a crude approximation is made to the integrals. ! For instance, the area of the flow region, LIQUID(T), is computed ! by summing up the estimated areas of each control volume. These ! areas, in turn, are estimated by computing the jacobian at the ! center of the control volume. ! ! +-----V-----+ ! | | ! | | ! U N U ! | | ! | | ! +-----V-----+ ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, double precision AREA(NI,NJ). ! AREA(I,J) is the area of the control volume (I,J), ! centered on primary node (I,J), with corner nodes ! (I,J), (I+1,J), (I,J+1) and (I+1,J+1). ! ! Output, double precision CINC. ! CINC is the value of F(T) at the current time, that is, the ! integral of the cost function over the current region. ! ! Input, double precision GAM(NI,NJ). ! GAM(I,J) contains the transport coefficient for the current equation ! at primary node (I,J). GAM(I,J) is general the fluid viscosity ! (equations for variables 1, 2, 3, or 4), or the thermal diffusion ! coefficient (equations for variable 5). ! GAM is only needed for cost function ICOST = 5. If SETCST is ! called with ICOST = 5, then routine GAMSOR