c============================================================================== c c This file contains the routines associated with the nonlinear c problem "Swirling Flow III", [Ascher, Mattheij, and Russell, c "Numerical Solution of Boundary Value Problems for Ordinary c Differential Equations", Classics in Applied Mathematics Series, c Society for Industrial and Applied Mathematics, Philadelphia, 1995]. c This is the special case of counter rotating disks - thus Sigma_0=-1 c and Sigma_1=1. c c============================================================================== c ***Model of main program to use MIRKDC, for Swirling Flow III test problem*** c============================================================================== c c ***declaration of constants*** c integer neqns, MxNsub, MxNprocs, ldY parameter(neqns=6,MxNsub=12000,MxNprocs=8,ldY=6) c c `neqns' Number of differential equations. c `MxNsub' Maximum number of subintervals. c `MxNprocs' Maximum number of processors/partitions. c `ldY' Leading dimension of the matrix Y. c c ***declaration of remaining parameters to mirkdc*** c integer method double precision tol integer Nsub integer mesh_input double precision mesh(0:MxNsub) integer sol_input double precision Y(ldY,(MxNsub+1)) integer output_control integer info double precision work(20*MxNsub*neqns + + 5*neqns*(MxNsub+1) + 2*neqns**2*(MxNsub+1) + + 8*(neqns + 22*neqns**2) + c c additional workspace needed for parallel RSCALE is c (MxNsub+2*MxNprocs+1)*neqns**2 c + (MxNsub+2*MxNprocs+1)*neqns**2 + c c additional workspace needed for parallel COUPLE is c MxNprocs*neqns**2 c + MxNprocs*neqns**2) c c (Actual dimension of work depends on Mxs=10 and nprocs=8. c double precision work(2*MxNsub*Mxs*neqns + 5*neqns*(MxNsub+1)+ c + 2*neqns**2*(MxNsub+1) + nprocs*(neqns + 2*neqns**2*(Mxs+1))) ) c integer iwork(neqns*(MxNsub+1)) integer nprocs c c `method' MIRK method to be used. c `tol' User-defined tolerance applied to the defect. c `Nsub' Number of subintervals. c `mesh_input' Communication parameter for initial mesh info. c `mesh' Set of points that partition the problem interval. c `sol_input' Communication parameter for initial solution info. c `Y' Discrete approximation to the solution. c `output_control' Controls the amount of information output. c `info' Communication flag for termination info. c `nprocs' Number of processors/partitions. c c ***declaration of user-supplied subroutines*** c external init_de,init_Y,fsub,gsub,dfsub,dgsub c c See descriptions below. c c ***declaration of local variables*** c c integer i,j,fail c double precision x,h c double precision z(neqns),zp(neqns) c c `i,j' Loop indexes c `fail' Flag to control execution of `sol_eval' c `x' Local variable to store a meshpoint c `h' Local variable to store the size of a subinterval c `z' Value of the interpolant evaluated at a mesh point c `zp' Value of the derivative of the interpolant evaluated at a c mesh point c c ***declaration of problem dependent parameter*** c double precision epsilon, lft, rgt common /ode/ epsilon, lft, rgt c---------------------------------------------------------------------------- c***declaration of variables used in parameter continuation c integer mx_num_conts parameter (mx_num_conts = 10) double precision cont_eps(mx_num_conts) integer k, nconts logical all_iters c------------------------------------------------------------------------- c***nonLAstats (non-[Linear Algebra] Statistics) common block c integer resid_count, defect_count, newmat_count real resid_time, defect_time, newmat_time common /nonLAstats/ resid_time, defect_time, newmat_time, * resid_count, defect_count, newmat_count c------------------------------------------------------------------------- c***LAstats (Linear Algebra Statistics) common block c integer mx_num_mesh, mx_num_fails parameter (mx_num_mesh = 25, mx_num_fails = 10) real counts(mx_num_mesh,5) integer fails(mx_num_mesh,mx_num_fails) common /LAstats/ counts, fails c------------------------------------------------------------------------- c***TOTALstats (nonLA + LA + [segments not monitored]) common block c real total_time common /TOTALstats/ total_time c------------------------------------------------------------------------- c***hybrid common block (used with hybrid LA kernel only) c integer LAselCriteria common /hybrid/ LAselCriteria c------------------------------------------------------------------------- c***declaration of parameters associated with the `graphData' array c (This array holds data for all runs in a compressed format; c the data is later processed and graphed using Matlab software.) c integer nGraph, mxnGraph, nGraphFields parameter (mxnGraph = 64, nGraphFields = 20) real graphData(mxnGraph,nGraphFields) c------------------------------------------------------------------------- write(6,*) ' swf_iii begin execution ... ' c c ******Setup parameters****** c c Prior to each execution, input the order of the MIRK scheme, c defect tolerance, number of subintervals in the initial mesh, c epsilon, left and right endpoints of the interval of integration, c # of processors (# of partitions defaults to the # of processors), c output control, linear algebra sequential-to-parallel control c (this switch is used with the hybrid LA kernel only), and number c of parameter continuations to be used. Parameter continuation c specifications are conveyed in the input as follows: c c (a) if (nconts .eq. 1 .or. nconts .eq. 0) parameter continuation c is not used; i.e. we attempt to solve the problem for the c given value of epsilon with a single call to `mirkdc' c c (b) if (iabs(nconts) .gt. 1) parameter continuation is used; c the next line of input contains `nconts' values for epsilon, c in descending order of magnitude, with the last value the c desired epsilon for the final solution, and sign(nconts) c specifies how statistics will be recorded: c (1) if (nconts .gt. 0) statistics are accumulated over c all `nconts' continuation problems c (2) if (nconts .lt. 0) statistics are recorded for the c last continuation problem only c nGraph = 0 10 read(5,*,end=30) method, tol, Nsub, epsilon, lft, rgt, + nprocs, output_control, LAselCriteria, nconts write(6,9999) method, tol, Nsub, epsilon, lft, rgt, + nprocs, output_control, LAselCriteria 9999 format (/1x, 'METHOD', 2x, 'DEFECT_TOL', 2x, 'NSUB_0', + 2x, 'EPSILON', 6x, 'LFT', 8x, 'RGT', 4x, 'NPROCS', + 2x, 'PROF', 2x, 'L/A_CTL' /1x, i4, 4x, d10.3, + 3x, i4, 3(1x, d10.3), 3x, i2, 6x, i1, 6x, i1/) if (nconts .eq. 1 .or. nconts .eq. 0) then nconts = 1 write(6,9997) 9997 format (' Parameter continuation will not be used.'/) cont_eps(nconts) = epsilon all_iters = .true. else if (nconts .gt. 0) then all_iters = .true. else all_iters = .false. endif nconts = iabs(nconts) read(5,*) (cont_eps(k),k=1,nconts) write(6,9996) nconts, (cont_eps(k),k=1,nconts) if (all_iters) then write(6,9994) nconts else write(6,9992) nconts endif write(6,*) 9996 format (' Parameter continuation will be invoked ', i2, + ' times;'/ ' selected epsilon values: ', 10d10.3) 9994 format (' Reported statistics are accumulated over all ', + i2, ' problems.') 9992 format (' Reported statistics apply to problem ', + i2, ' only.') end if c c Set mesh_input to indicate that mirkdc should c begin with a uniform mesh. mesh_input = 0 c c Set indicator for type of initial solution estimate. sol_input = 1 c------------------------------------------------------------------------- c ******Initialize processors****** c call init_procs(nprocs) c------------------------------------------------------------------------- c ******Initialize non-LA, LA and TOTAL stats****** c if (all_iters) then call init_nonLAstats call init_LAstats call init_TOTALstats endif c------------------------------------------------------------------------- c ******Parameter continuation loop. Statistics are either c (1) accumulated over all `nconts' iterations, or c (2) apply to the last iteration only. c if (all_iters) call update_TOTALstats(0) do 20 k = 1, nconts epsilon = cont_eps(k) write(6,9990) epsilon 9990 format (' ... attempting solution for epsilon = ', d10.3) c------------------------------------------------------------------------- c ******Initialize non-LA, LA and TOTAL stats****** c if (.not. all_iters) then call init_nonLAstats call init_LAstats call init_TOTALstats endif c------------------------------------------------------------------------- c ******Call mirkdc****** c if (.not. all_iters) call update_TOTALstats(0) call mirkdc(method,tol,neqns,Nsub,MxNsub,mesh_input, + mesh,sol_input,Y,ldY,output_control,info,iwork,work, + init_de,init_Y,fsub,gsub,dfsub,dgsub,nprocs) if (.not. all_iters) call update_TOTALstats(1) c c ******Check return from mirkdc****** c if (info .NE. 0) then write(6,*) ' Solution not computed: too many mesh ', + 'subintervals required.' if (nconts .gt. 1) then write(6,*) ' Parameter continuation aborted!' go to 25 endif else write(6,*) ' Solution computed to within specified ', + 'defect tolerance.' endif mesh_input = 2 sol_input = 2 20 continue 25 if (all_iters) call update_TOTALstats(1) c------------------------------------------------------------------------- c ******Display non-LA, LA and TOTAL stats for this run****** c call display_nonLAstats call display_LAstats call display_TOTALstats(0,method,tol,Nsub,epsilon,lft,rgt, + nprocs,LAselCriteria,nGraph,mxnGraph,nGraphFields,graphData) c------------------------------------------------------------------------- go to 10 30 continue c------------------------------------------------------------------------- c ******Write out the `graphData' array to be processed c and graphed later using Matlab software c call display_TOTALstats(1,method,tol,Nsub,epsilon,lft,rgt, + nprocs,LAselCriteria,nGraph,mxnGraph,nGraphFields,graphData) c c Upon successful return from mirkdc, print out solution. c Solution interpolant is available through the sol_eval routine. c c Compute solution approximations at 11 points. c c h = (mesh(Nsub) - mesh(0))/10.0d0 c do 40 i = 1, 11 c x = (mesh(0) + (i-1)*h) c call sol_eval(x,z,zp,fail,iwork,work) c if (fail .EQ. 0) then c write(6,*) 'At meshpoint x=', x,' the solution is' c do 35 j = 1, neqns c write(6,100) z(j) c100 format (1x, d30.15) c35 continue c else c write(6,*) 'Attempted to evaluate the interpolant at a point', c + 'that was outside of the problem interval. No information', c + 'is available for that point.' c end if c40 continue write(6,*) ' ... swf_iii end execution. ' stop end c c============================================================================== c ***Models of user-defined subroutines for Swirling Flow III test problem*** c============================================================================== c c c Original Form First Order System Form c ------------- ----------------------- c y1' = y2, c g'' = gf' - fg' c --------- , y2' = y1*y4 - y3*y2 c epsilon ------------- , c epsilon c f'''' = -ff''' - gg' y3' = y4 , c ------------ , c epsilon y4' = y5 , y5' = y6 , c c y6' = -y3*y6 - y1*y2 c -------------- , c epsilon c with c c g(0) = -1, g(1) = 1, y1(0) = -1, y1(1) = 1, c c f(0) = f'(0) = f(1) = f'(1) = 0, y3(0) = y4(0) = 0, c y3(1) = y4(1) = 0, c and c epsilon = {a small parameter}. c c No closed form solution is available. c============================================================================== c subroutine init_de(leftbc,a,b) c c The purpose of this routine is to initialize some of the problem c dependent parameters. c declaration of parameters: c exports: integer leftbc double precision a,b c c `leftbc', Number of boundary conditions at the left c endpoint of the problem interval. c `a', Left endpoint of the problem interval. c `b', Right endpoint of the problem interval. c c*** lft' and `rgt' now are input by the driver and stored, c*** along with `epsilon', in common block /ode/. c double precision epsilon, lft, rgt common /ode/ epsilon, lft, rgt c------------------------------------------------------------------------------ c leftbc = 3 c a = lft b = rgt c return end c c============================================================================= c subroutine init_Y(Nsub,neqns,mesh,Y,ldY) c c The purpose of this routine is to initialize the discrete solution c approximation. c c declaration of parameters: c imports: integer Nsub, neqns, ldY double precision mesh(0:Nsub) c c `Nsub' Number of subintervals. c `neqns' Number of differential equations. c `ldY' Leading dimension of `Y'. c `mesh' Set of points that partition the problem interval. c c exports: double precision Y(ldY,(Nsub+1)) c c `Y', Discrete solution approximation at meshpoints. The i-th c column of `Y' is the solution approximation vector for the c mesh point, `mesh(i-1)'. c c declaration of local variables: integer i c c `i' loop index from 0 to Nsub c c------------------------------------------------------------------------------ c c The initial approximation is based on straight lines joining the c boundary conditions for each differential equation in the first c order system. Boundary conditions are given for the first, c third and fourth equations, defining the lines y = t (or mesh(i)) c y = 0 and y = 0. Since there are no boundary conditions for the second, c fifth, or sixth equations, the line drawn for y2 is y = 1 (the c derivative of y = t since y1' = y2 ) and y5 and y6 are initialized c by the line y = 0. c do 10 i = 0, Nsub Y(1, i+1) = mesh(i) Y(2, i+1) = 1.d0 Y(3, i+1) = 0.d0 Y(4, i+1) = 0.d0 Y(5, i+1) = 0.d0 Y(6, i+1) = 0.d0 10 continue c return end c c============================================================================= c subroutine fsub(neqns,t,y,f) c c This routine defines f(t,y) for the first order system of c differential equations. c c declaration of parameters: c imports: integer neqns double precision t, y(neqns) c c `neqns' Number of differential equations. c `t' A point in the problem interval where f(t,y(t)) is to be c evaluated. c `y' Current solution approximation at `t', used in the c evaluation of f(t,y(t)). c c exports: double precision f(neqns) c c `f' Value of f(t,y(t)) for given `t' and `y'. c c declaration for parameter 'epsilon' double precision epsilon, lft, rgt common /ode/ epsilon, lft, rgt c c------------------------------------------------------------------------------ c f(1) = y(2) f(2) = (y(1)*y(4) - y(3)*y(2))/epsilon f(3) = y(4) f(4) = y(5) f(5) = y(6) f(6) = (-y(3)*y(6) - y(1)*y(2))/epsilon c return end c c============================================================================== c subroutine dfsub(neqns,t,y,JJ) c c This routine defines the Jacobian, d f(t,y(t))/dy, of the system of c differential equations. Since `JJ' has already been set to zero, it is c only necessary to set the non-zero elements. c c------------------------------------------------------------------------------ c c declaration of parameters: c imports: integer neqns double precision t, y(neqns) c c `neqns' Number of differential equations. c `t' A point in the problem interval where d f(t,y(t))/dy is to c be evaluated. c `y' Current solution approximation at `t', used in the c evaluation of d f(t,y(t))/dy. c c exports: double precision JJ(neqns,neqns) c c `JJ' Value of d f(t,y(t))/dy for given `t' and `y'. c c declaration for parameter 'epsilon' double precision epsilon, lft, rgt common /ode/ epsilon, lft, rgt c c------------------------------------------------------------------------------ c JJ(1,2) = 1.0d0 c JJ(2,1) = y(4)/epsilon JJ(2,2) = -y(3)/epsilon JJ(2,3) = -y(2)/epsilon JJ(2,4) = y(1)/epsilon c JJ(3,4) = 1.0d0 c JJ(4,5) = 1.0d0 c JJ(5,6) = 1.0d0 c JJ(6,1) = -y(2)/epsilon JJ(6,2) = -y(1)/epsilon JJ(6,3) = -y(6)/epsilon JJ(6,6) = -y(3)/epsilon c return end c c============================================================================== c subroutine gsub(neqns,ya,yb,bc) c c This routine defines the boundary condition equations. c (MIRKDC requires the boundary conditions equations to be written with c zero on the right hand side of each equation; thus, e.g. c yb(1) = 1 becomes yb(1) - 1 = 0 and bc(4) = yb(1) - 1 .) c c------------------------------------------------------------------------------ c c declaration of parameters: c imports: integer neqns double precision ya(neqns),yb(neqns) c c `neqns' Number of differential equations. c `ya' Current solution approximation at `a', used in the c evaluation of g_a(y(a)). c `yb' Current solution approximation at `b', used in the c evaluation of g_b(y(b)). c c exports: double precision bc(neqns) c c `bc' The first `leftbc' locations contain g_a(y(a)). c The remaining `neqns-leftbc' locations contain g_b(y(b)). c c------------------------------------------------------------------------------ c bc(1) = ya(1) + 1.0d0 bc(2) = ya(3) bc(3) = ya(4) c bc(4) = yb(1) - 1.0d0 bc(5) = yb(3) bc(6) = yb(4) c return end c c============================================================================== c subroutine dgsub(neqns,ya,yb,dbc) c c This routine defines the Jacobian, d bc/dy, of the boundary conditions. c `dbc' has already been set to zero so it is only necessary to set c the non-zero elements. c c------------------------------------------------------------------------------ c c declaration of parameters: c imports: integer neqns double precision ya(neqns),yb(neqns) c c `neqns' Number of differential equations. c `ya' Current solution approximation at `a', used in the c evaluation of d g_a(y(a))/d ya. c `yb' Current solution approximation at `b', used in the c evaluation of d g_b(y(b))/ d yb. c c exports: double precision dbc(neqns,neqns) c c `dbc' Value of d bc/dy for given `t',`ya', and `yb'. c c--------------------------------------------------------------------------- c dbc(1,1) = 1.0d0 dbc(2,3) = 1.0d0 dbc(3,4) = 1.0d0 dbc(4,1) = 1.0d0 dbc(5,3) = 1.0d0 dbc(6,4) = 1.0d0 c return end