program main !*****************************************************************************80 ! !! MAIN is the main program for FLOW7. ! ! Discussion: ! ! FLOW7 solves a 2D steady incompressible flow using finite elements. ! ! FLOW7 is a version of FLOW6 that is being modified to include the option ! of iterative solution of the linear system. ! ! FLOW6 is an experimental version of FLOW5. A new free slip boundary ! condition is being added. ! ! FLOW5 is a pared down and simplified version of the research code FLOW4. ! ! FLOW4 includes parameterization, sensitivity analysis, the use of ! isoparametric elements, a curved bump on the bottom of the channel region, ! more choices for quadrature points and other features which have ! been removed from FLOW5. ! ! This program was developed by John Burkardt, based on programs ! of Janet Peterson and Max Gunzburger. ! ! ! The program works on an underlying fluid flow problem, whose ! behavior is determined by a particular version of the Navier Stokes ! equations. ! ! The fluid flow in the region is described by three state functions: ! ! U(X,Y), the horizontal velocity, ! V(X,Y), the vertical velocity, ! P(X,Y), the pressure. ! ! In theory, these functions may be determined once we know the partial ! differential equations that govern them within the region, and ! the value of the functions or certain derivatives of them along the ! boundary of the region. ! ! For our work, we assume that at every point within the flow region, ! the functions obey the Navier Stokes equations for stationary, ! incompressible, viscous flow: ! ! - nu * ( ddU/dxdx + ddU/dydy ) + U * dU/dx + V * dU/dy + dP/dx = 0 ! ! - nu * ( ddV/dxdx + ddV/dydy ) + U * dV/dx + V * dV/dy + dP/dy = 0 ! ! dU/dx + dV/dy = 0 ! ! The first two equations are sometimes called the "momentum" equations, ! and the third the "continuity" equation. Nu is a physical parameter called ! the "dynamic viscosity". ! ! We prefer the equivalent formulation (when nu is nonzero): ! ! - ddU/dxdx - ddU/dydy + nu_inv * ( U * dU/dx + V * dU/dy + dP/dx ) = 0 ! ! - ddV/dxdx - ddV/dydy + nu_inv * ( U * dV/dx + V * dV/dy + dP/dy ) = 0 ! ! dU/dx + dV/dy = 0 ! ! where nu_inv = ( 1 / nu ). ! ! ! To complete the specification of the problem, we specify boundary ! conditions for the flow functions. ! !*****************************************************************************80 ! ! THE ROLE OF THE DYNAMIC VISCOSITY: ! ! Nu is a physical parameter called the "dynamic viscosity". We explicitly ! assume that nu is not zero. In the momentum equations, Nu multiplies ! the Laplacian operator, and thus, determines the relative weight of ! the smoothing or diffusion terms, as compared to the nonlinear terms. ! ! For a particular fluid, Nu is a physical parameter, which is usually ! taken to be a constant. When comparing two fluids, the fluid with ! the higher value of Nu will behave more like syrup; velocity tends ! to diffuse; a small moving particle will slow down, feeling a heavy drag, ! while neighboring particles begin to move along. A fluid with a lower ! value of Nu is more "slippery", making it possible for large variations ! in velocity to occur over a short distance, with less of a tendency ! to be damped out. ! ! Engineers use a combination of the dynamic viscosity, a typical ! length, and a typical velocity, to define the Reynolds number for ! a given flow: ! ! Re = v * L / nu ! ! As the Reynolds number increases, (because, perhaps, the average ! velocity is increasing), the character of the flow will tend to ! change from a regular laminar flow to an unsteady turbulent flow. ! A large value of Nu delays this transition. ! ! Correspondingly, as the Reynolds number gets larger, the mathematical ! problem becomes more difficult to solve, because the nonlinear part ! of the momentum equations becomes predominant. For large enough ! Reynolds number there may be no solution, or multiple solutions. ! !*****************************************************************************80 ! ! DERIVATION OF FINITE ELEMENT EQUATIONS ! ! Except for special cases, such as the Poiseuille flow solution ! discussed elsewhere, there are no methods of producing the exact ! solution functions U, V and P for a general Navier Stokes problem. ! In order to get any insight into flow problems, we must replace the ! original problem by one that is much weaker. It's important that the ! weaker problem be possible to solve, and that the solutions produced ! are in general close to solutions of the original problem, and that ! these solutions can be made even closer,if desired. ! ! A standard method of doing this uses finite elements. ! ! To do so, we assume that instead of being smooth but otherwise ! completely arbitrary functions, that U, V and P are representable ! as linear combinations of a finite set of basis functions. ! ! We multiply the first two equations by an arbitrary velocity basis ! function Wi, and the third equation by an arbitrary pressure basis ! function Qi, and integrate over the region. The integrand of the ! resulting finite element equations is then transformed, using ! integration by parts, into: ! ! UM-Eqn(I): ! ! ( dU/dx * dW(I)/dx + dU/dy * dW(I)/dy ) ! + nu_inv * ( U * dU/dx + V * dU/dy + dP/dx ) * W(I) ! ! VM-Eqn(I): ! ! ( dV/dx * dW(I)/dx + dV/dy * dW(I)/dy ) ! + nu_inv * ( U * dV/dx + V * dV/dy + dP/dy ) * W(I) ! ! PC-Eqn(I): ! ! ( dU/dx + dV/dy ) * Q(I) ! ! These integrands may be rewritten using the program's variable names: ! ! dUdx * dwidx + dUdy * dwidy ! + nu_inv * ( U * dUdx + V * dUdy + dPdx ) * wi ! ! dVdx * dwidx + dVdy * dwidy ! + nu_inv * ( U * dVdx + V * dVdy + dPdy ) * wi ! ! ( dUdx + dVdy ) * qi ! ! This system of nonlinear equations is then solved by Newton's method. ! That means that we have to differentiate each nonlinear equation ! with respect to the unknowns, getting the Jacobian matrix, and ! solving DF(X) * DEL(X) = -F(X). If we abuse notation, we can ! consider the linear system DF(X) * DEL(X). ! ! Here, variables U, V and P in capital letters are to be solved for, ! but the same variable names in lowercase represent the current ! values of those same variables. ! ! d UM-Eqn(I) / d U-coefficient * U coefficient: ! ! dUdx * dwidx + dUdy * dwidy ! + nu_inv * ( U * dudx + u * dUdx + v * dUdy ) * wi ! ! d UM-Eqn(I) / d V coefficient * V coefficient: ! ! nu_inv * V * dudy * wi ! ! d UM-Eqn(I) / d P coefficient * P coefficient: ! ! nu_inv * dPdx * wi ! ! d VM-Eqn(I) / d U coefficient * U coefficient: ! ! nu_inv * U * dvdx * wi ! ! d VM-Eqn(I) / d V coefficient * V coefficient: ! ! dVdx * dwidx + dVdy * dwidy ! + nu_inv * ( u * dVdx + v * dVdy + V * dvdy ) * wi ! ! d VM-Eqn(I) / d P coefficient * P coefficient: ! ! nu_inv * dPdy * wi ! ! d PC-Eqn(I) / d U coefficient * U coefficient: ! ! dUdx * qi ! ! d PC-Eqn(I) / d V coefficient * V coefficient: ! ! dVdx * qi ! !*****************************************************************************80 ! ! Modified: ! ! 03 August 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Max Gunzburger, ! Finite Element Methods for Viscous Incompressible Flows, ! A Guide to Theory, Practice, and Algorithms, ! Academic Press, 1989, ! LC: TA357.G86 ! ! Parameters: ! ! real A(NROW,MAXEQN), the value of D F(I)/D X(J) ! for each of the NEQN residual functions F(I) with respect to each ! of the unknown coefficients X(J). The information is stored as ! a LINPACK/LAPACK general band format matrix. ! ! real DETMAP(MAXQUAD2,MAXELM), the reference map determinant. ! ! character ( len = 3 ) EQN(MAXEQN), the type of each equation. ! 'UM' the horizontal momentum equation; ! 'UI' horizontal velocity specified at node. ! 'U0' U = 0 at node. ! 'VM' the vertical momentum equation; ! 'VI' V = specified value at node. ! 'V0' V = 0 at node. ! 'N0' the velocity component normal to the wall is 0. ! 'T0' the velocity component tangential to the wall is 0. ! 'TMP' momentum equation for the velocity component tangential ! to the wall, plus penalty term. ! 'PC' the continuity equation; ! 'PI' P = specified value at node. ! 'P0' P = 0 at node. ! ! 'FS' free slip, Integral (U,V) dot (boundary tangent) is penalized; ! ! real ETAQUAD(MAXQUAD2), the ETA quadrature coordinates. ! ! real G(MAXEQN), the finite element coefficients. ! ! integer IERROR, 0/nonzero, an error DID NOT/DID occur. ! ! integer INDX(3,MAXNP), the mapping from nodes to degrees of freedom. ! INDEX(1,J) = horizontal velocity degree of freedom index. ! INDEX(2,J) = vertical velocity degree of freedom index. ! INDEX(3,J) = pressure degree of freedom index, or 0 if no pressure. ! ! integer IPIVOT(MAXEQN), pivoting space for the linear system solver. ! ! integer JAC, analytic/finite difference jacobian option. ! 0, use analytic jacobian. ! 1, use finite difference jacobian. ! ! integer MAXELM, the maximum number of elements, = 2 * MAXNX * MAXNY. ! ! integer MAXEQN, the maximum number of equations. ! ! integer NEWTON_MAX, the maximum number of Newton iterations. ! ! integer MAXNP, the maximum number of nodes. ! ! integer MAXNX, the maximum number of elements in the X direction. ! ! integer MAXNY, the maximum number of elements in the Y direction. ! ! integer MAXQUAD1, the maximum number of 1D quadrature points. ! ! integer MAXQUAD2, the maximum number of 2D quadrature points. ! ! integer MAXROW, the first dimension of the matrix A. ! ! integer NELEM, the number of elements. ! ! integer NEQN, the number of finite element equations. ! ! integer NLBAND, the lower bandwidth of the matrix A. ! ! integer NODE(6,MAXELM), the nodes that make up each element. ! ! integer NP, the number of nodes. ! ! integer NQUAD1, the number of 1D quadrature points. ! ! integer NQUAD2, the number of 2D quadrature points. ! ! integer NROW, the number of rows needed to store the matrix A. ! ! integer NX, the number of elements along a horizontal line. ! ! integer NY, the number of elements along a vertical line. ! ! real PENALTY1, the TMP equation friction penalty coefficient. ! ! real PENALTY2, the TMP momentum penalty coefficient. ! ! real PHI(MAXQUAD2,6,6,MAXELM), basis function values. ! PHI contains the value of a basis function, its derivative, ! or other information, evaluated at a quadrature point. ! For a particular element I, quadrature point J, and basis ! function K, we use the following shorthand for the ! entries of PHI: ! W, dWdX, dWdY ! Q, dQdX, dQdY ! W is the quadratic basis function associated with velocity, ! Q the linear basis function associated with pressure, ! Xsi and Eta the reference coordinates for the point. ! In particular, PHI(J,K,1,I) is the value of the quadratic ! basis function associated with local node K in element I, ! evaluated at quadrature point J. ! ! real PNORM_LMAX, the max-norm of the pressure. ! ! real PNORM_L2, the L2-norm of the pressure. ! ! character ( len = 20 ) REGION, the flow problem, 'CAVITY', 'CHANNEL', ! 'FREESLIP', or 'STEP'. ! ! real RES(MAXEQN), the finite element Newton residual. ! ! real NU_INV, the inverse viscosity. ! ! character ( len = 20 ) SOLVER, the linear system solver, 'GAUSS' or 'CGS'. ! ! real SQUAD1(MAXQUAD1), the 1D reference quadrature coordinates. ! ! real NEWTON_TOLERANCE, the Newton convergence tolerance. ! ! real UVNORM_LMAX, the max-norm of the velocity magnitude. ! ! real UVNORM_L2, the L2 norm of the velocity magnitude. ! ! real WQUAD1(MAXQUAD1), the 1D quadrature weights. ! ! real WQUAD2(MAXQUAD2), the 2D quadrature weights. ! ! real XC(MAXNP), the X coordinates of the nodes. ! ! real XRANGE, the width of the region. ! ! real XSIQUAD(MAXQUAD2), the XSI quadrature coordinates. ! ! real YC(MAXNP), the Y coordinates of the nodes. ! ! real YRANGE, the height of the region. ! implicit none integer, parameter :: maxnx = 20 integer, parameter :: maxny = 20 integer, parameter :: maxquad1 = 3 integer, parameter :: maxquad2 = 7 ! ! The assignment of MAXROW assumes that the nodes are ordered starting ! at the bottom left corner, proceeding upwards in a column, and then ! moving back to the bottom of the next column to the right. For ! our choice of elements, this allows us to estimate the greatest ! difference between indices of two variables which occur in the ! same equation. ! integer, parameter :: maxrow = 29 * ( maxny + 1 ) integer, parameter :: maxelm = 2 * maxnx * maxny integer, parameter :: maxeqn = 2 * ( 2 * maxnx + 1 ) & * ( 2 * maxny + 1 ) + ( maxnx + 1 ) * ( maxny + 1 ) integer, parameter :: maxnp = ( 2 * maxnx + 1 ) & * ( 2 * maxny + 1 ) integer, parameter :: maxside = 4 * maxnx * maxny integer, parameter :: node_nabor_max = 19 integer, parameter :: eqn_nabor_max = 45 real ( kind = 8 ) a(maxrow,maxeqn) integer bc_tag(maxnp) real ( kind = 8 ) detmap(maxquad2,maxelm) character ( len = 3 ) eqn(maxeqn) integer eqn_nabor(eqn_nabor_max,maxeqn) integer eqn_nabor_num(maxeqn) real ( kind = 8 ) etaquad(maxquad2) character ( len = 20 ) file_name real ( kind = 8 ) g(maxeqn) real ( kind = 8 ) g2(maxeqn) integer i integer ierror integer indx(3,maxnp) integer ipivot(maxeqn) integer iplot integer iprint integer jac integer ncol integer node_nabor(node_nabor_max,maxnp) integer node_nabor_num(maxnp) integer nelem integer neqn integer neqn2 integer newton_max character ( len = 10 ) newton_start integer newton_stutter real ( kind = 8 ) newton_tolerance integer nlband integer node(6,maxelm) integer np integer nquad1 integer nquad2 integer nrow integer nside integer nx integer nx2 integer ny integer ny2 real ( kind = 8 ) p(maxnp) real ( kind = 8 ) penalty1 real ( kind = 8 ) penalty2 real ( kind = 8 ) phi(maxquad2,6,6,maxelm) real ( kind = 8 ) pnorm_h1 real ( kind = 8 ) pnorm_l2 real ( kind = 8 ) pnorm_lmax character ( len = 20 ) region character ( len = 20 ) region2 real ( kind = 8 ) region_xmax real ( kind = 8 ) region_ymax real ( kind = 8 ) res(maxeqn) real ( kind = 8 ) res2(maxeqn) real ( kind = 8 ) nu_inv real ( kind = 8 ) nu_inv2 logical s_eqi integer side_basis(3,maxside) integer side_elem(maxside) character ( len = 3 ) side_eqn(maxside) real ( kind = 8 ) side_etam(maxside) real ( kind = 8 ) side_etap(maxside) integer side_indx(3,maxside) real ( kind = 8 ) side_xsim(maxside) real ( kind = 8 ) side_xsip(maxside) character ( len = 20 ) solver real ( kind = 8 ) squad1(maxquad1) character ( len = 20 ) system real ( kind = 8 ) uvnorm_h1 real ( kind = 8 ) uvnorm_l2 real ( kind = 8 ) uvnorm_lmax real ( kind = 8 ) wquad1(maxquad1) real ( kind = 8 ) wquad2(maxquad2) real ( kind = 8 ) xc(maxnp) real ( kind = 8 ) xsiquad(maxquad2) real ( kind = 8 ) yc(maxnp) call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' Last modified on 09 May 2000.' ! ! Say hello. ! call hello ( maxelm, maxeqn, maxnp, maxnx, maxny, maxquad2, maxrow, maxside ) ! ! Zero out the variables. ! call init ( a, bc_tag, detmap, eqn, etaquad, g, ierror, indx, ipivot, jac, & maxelm, maxeqn, maxnp, maxquad1, maxquad2, maxrow, maxside, nelem, neqn, & newton_max, newton_stutter, newton_tolerance, nlband, node, np, nquad1, & nquad2, nrow, nside, nx, ny, p, penalty1, penalty2, phi, pnorm_l2, & pnorm_lmax, region, res, nu_inv, side_basis, side_elem, side_eqn, & side_etam, side_etap, side_indx, side_xsim, side_xsip, solver, squad1, & uvnorm_l2, uvnorm_lmax, wquad1, wquad2, xc, xsiquad, yc ) ! ! Now we set some values which will define all the others. ! These values could be read in from a file. ! region = 'CAVITY' ! ! Specify: ! ! JAC, 0 to use jacobian routine, nonzero to estimate the jacobian. ! NEWTON_STUTTER, number of times jacobian is used before recomputed; ! NX, the number of elements in the horizontal direction, and ! NY, the number of elements in the vertical direction; ! PENALTY1, the tangential flow penalty coefficient applied to TMP nodes; ! PENALTY2, the momentum penalty coefficient applied to TMP nodes; ! REGION_XMAX, the width and ! REGION_YMAX, the height of the region; ! NU_INV, the value of the inverse of the viscosity; ! SYSTEM, 'NAVIER-STOKES' or 'STOKES'. ! ! For the cavity, it would be nice if ! ! NX = NY. ! if ( region == 'CAVITY' ) then jac = 0 newton_stutter = 1 nx = 3 ny = 3 penalty1 = 0.0D+00 penalty2 = 1.0D+00 region_xmax = 1.0D+00 region_ymax = 1.0D+00 nu_inv = 10.0D+00 solver = 'GAUSS' ! solver = 'CGS' system = 'NAVIER-STOKES' ! ! For the channel, it would be nice if ! ! NX / NY = REGION_XMAX / REGION_YMAX ! else if ( region == 'CHANNEL' ) then jac = 0 newton_stutter = 1 nx = 10 ny = 3 penalty1 = 0.0D+00 penalty2 = 1.0D+00 region_xmax = 10.0D+00 region_ymax = 3.0D+00 nu_inv = 10.0D+00 solver = 'GAUSS' system = 'NAVIER-STOKES' ! ! For the freeslip channel, it would be nice if ! ! NX / NY = REGION_XMAX / REGION_YMAX ! else if ( region == 'FREESLIP' ) then jac = 0 newton_stutter = 1 nx = 10 ny = 3 penalty1 = 1.0D+00 penalty2 = 0.0D+00 region_xmax = 10.0D+00 region_ymax = 3.0D+00 nu_inv = 10.0D+00 solver = 'GAUSS' system = 'NAVIER-STOKES' ! ! For the step, it would be nice if ! ! NX / NY = REGION_XMAX / REGION_YMAX ! ! and if ! ! NX is divisible by 8. ! else if ( region == 'STEP' ) then jac = 0 newton_stutter = 1 nx = 32 ny = 8 penalty1 = 1.0D+00 penalty2 = 0.0D+00 region_xmax = 8.0D+00 region_ymax = 2.0D+00 nu_inv = 10.0D+00 solver = 'GAUSS' system = 'NAVIER-STOKES' else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7 - Fatal error!' write ( *, '(a)' ) ' An unexpected region was specified.' write ( *, '(a)' ) ' REGION = ' // trim ( region ) stop end if ! ! Set up the problem geometry. ! call geometry ( bc_tag, eqn, indx, maxelm, maxeqn, maxnp, nelem, neqn, node, & np, nx, ny, region, region_xmax, region_ymax, xc, yc ) ! ! Determine the node neighbor array. ! if ( solver == 'CGS' ) then call node_nabor_set ( nelem, node, np, node_nabor, node_nabor_max, & node_nabor_num ) call node_nabor_print ( np, node_nabor, node_nabor_max, node_nabor_num ) end if ! ! Determine the equation neighbor array. ! if ( solver == 'CGS' ) then call eqn_nabor_set ( eqn_nabor, eqn_nabor_max, eqn_nabor_num, indx, & neqn, np, node_nabor, node_nabor_max, node_nabor_num ) call eqn_nabor_print ( eqn_nabor, eqn_nabor_max, eqn_nabor_num, neqn ) end if ! ! Compute boundary integral information. ! call boundary_integral ( maxside, nside, nx, ny, region, side_basis, & side_elem, side_eqn, side_etam, side_etap, side_indx, side_xsim, side_xsip ) ! ! Print data. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The flow region is ' // trim ( region ) write ( *, '(a)' ) ' The state equations are ' // trim ( system ) write ( *, '(a)' ) ' The linear solver is ' // trim ( solver ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The number of:' write ( *, '(a,i8)' ) ' Horizontal elements NX = ', nx write ( *, '(a,i8)' ) ' Vertical elements NY= ', ny write ( *, '(a,i8)' ) ' Elements, NELEM = ', nelem write ( *, '(a,i8)' ) ' Nodes, NP = ', np write ( *, '(a,i8)' ) ' Unknowns, NEQN = ', neqn write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Newton method parameters:' write ( *, '(a)' ) ' ' if ( jac == 0 ) then write ( *, '(a)' ) ' JAC = 0, use analytic jacobian.' else write ( *, '(a)' ) ' JAC nonzero, use finite difference jacobian.' end if write ( *, '(a,i8)' ) ' Iteration limit, NEWTON_MAX =', newton_max write ( *, '(a,i8)' ) ' Jacobian reuse, NEWTON_STUTTER =', newton_stutter write ( *, '(a,g14.6)' ) & ' Iteration tolerance, NEWTON_TOLERANCE =', newton_tolerance write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) & ' TMP friction term penalty coefficient PENALTY1 = ', penalty1 write ( *, '(a,g14.6)' ) & ' TMP momentum term penalty coefficient PENALTY2 = ', penalty2 iprint = 0 if ( iprint == 1 ) then call xy_print ( maxnp, np, xc, yc ) end if iprint = 0 if ( iprint == 1 ) then call element_print ( maxelm, nelem, node ) end if iprint = 0 if ( iprint == 1 ) then call equation_print ( eqn, indx, neqn, np ) end if iprint = 0 if ( iprint == 1 ) then call boundary_integral_print ( maxside, nside, side_basis, side_elem, & side_eqn, side_etam, side_etap, side_indx, side_xsim, side_xsip ) end if ! ! Determine the lower matrix bandwidth, and NROW, ! the number of rows we will use in the matrix A. ! if ( s_eqi ( solver, 'GAUSS' ) ) then call fp_band_width ( indx, maxelm, maxnp, nelem, nlband, node, nrow ) ncol = neqn if ( maxrow < nrow ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7 - Fatal error!' write ( *, '(a,i8)' ) ' NROW is too large! NROW = ', nrow write ( *, '(a,i8)' ) ' The maximum allowed is MAXROW = ', maxrow stop end if write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Lower bandwidth NLBAND = ', nlband write ( *, '(a,i8)' ) ' Available matrix rows, MAXROW = ', maxrow else if ( s_eqi ( solver, 'CGS' ) ) then nrow = neqn ncol = node_nabor_max end if write ( *, '(a,i8)' ) ' Required matrix rows NROW = ', nrow write ( *, '(a,i8)' ) ' Required matrix columns NCOL = ', ncol if ( solver == 'CGS' ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7 - Warning!' write ( *, '(a)' ) ' The CGS solver has been requested.' write ( *, '(a)' ) ' But this solver is not available.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Abnormal end of execution.' stop end if ! ! Set the boundary quadrature rule. ! call ref_quad1 ( maxquad1, nquad1, squad1, wquad1 ) ! ! Set the reference element quadrature rule. ! call ref_quad2 ( etaquad, maxquad2, nquad2, wquad2, xsiquad ) ! ! Set the value of the basis functions at all quadrature points. ! call setbas ( detmap, etaquad, maxelm, maxnp, maxquad2, nelem, node, & nquad2, phi, xc, xsiquad, yc ) ! ! Solve the linear STOKES system, or the nonlinear NAVIER-STOKES system. ! if ( system == 'STOKES' ) then call stokes ( a, bc_tag, detmap, eqn, g, ierror, indx, ipivot, jac, & maxelm, maxeqn, maxnp, maxquad1, maxquad2, maxside, ncol, nelem, neqn, & nlband, node, np, nquad1, nquad2, nrow, nside, penalty1, penalty2, phi, & region, region_ymax, res, res2, nu_inv, side_basis, side_elem, & side_eqn, side_etam, side_etap, side_indx, side_xsim, side_xsip, & solver, squad1, wquad1, wquad2, xc, yc ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7:' write ( *, '(a)' ) ' STOKES computed the Stokes solution.' else if ( system == 'NAVIER-STOKES' ) then ! ! Get an initial estimate of the solution. ! newton_start = 'ZERO' if ( newton_start == 'STOKES' ) then write ( *, '(a)' ) 'Get Newton start via Stokes solution.' call stokes ( a, bc_tag, detmap, eqn, g, ierror, indx, ipivot, jac, & maxelm, maxeqn, maxnp, maxquad1, maxquad2, maxside, ncol, nelem, & neqn, nlband, node, np, nquad1, nquad2, nrow, nside, penalty1, & penalty2, phi, region, region_ymax, res, res2, nu_inv, side_basis, & side_elem, side_eqn, side_etam, side_etap, side_indx, side_xsim, & side_xsip, solver, squad1, wquad1, wquad2, xc, yc ) else if ( newton_start == 'ZERO' ) then write ( *, '(a)' ) 'Set Newton starting point to zero.' g(1:neqn) = 0.0D+00 else if ( newton_start == 'FILE' ) then write ( *, '(a)' ) 'Get Newton starting point from file.' file_name = 'g.txt' call coef_read ( file_name, g, ierror, maxeqn, neqn2, nx2, ny2, & region2, nu_inv2 ) if ( ierror == 0 ) then write ( *, '(a,g14.6)' ) ' Saved solution has NU_INV = ', nu_inv2 end if if ( neqn2 /= neqn ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7 - Warning!' write ( *, '(a,i8)' ) ' Saved solution file has NEQN = ', neqn2 write ( *, '(a,i8)' ) ' But we need NEQN = ', neqn ierror = 1 end if if ( ierror /= 0 ) then g(1:neqn) = 0.0D+00 write ( *, '(a)' ) 'File reading failed. Setting G to zero.' end if end if ! ! Solve the nonlinear system. ! call newton ( a, bc_tag, detmap, eqn, g, ierror, indx, ipivot, jac, & maxelm, maxeqn, maxnp, maxquad1, maxquad2, maxside, ncol, nelem, & neqn, newton_max, newton_stutter, newton_tolerance, nlband, node, & np, nquad1, nquad2, nrow, nside, penalty1, penalty2, phi, region, & region_ymax, res, res2, nu_inv, side_basis, side_elem, side_eqn, & side_etam, side_etap, side_indx, side_xsim, side_xsip, solver, squad1, & wquad1, wquad2, xc, yc ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7 - Fatal error!' write ( *, '(a)' ) ' The Newton iteration failed!' stop end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7:' write ( *, '(a)' ) ' NEWTON computed the Navier Stokes solution.' else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7 - Fatal error!' write ( *, '(a)' ) ' Unrecognized state equation name.' stop end if call uvp_norm_lmax ( g, indx, maxeqn, maxnp, np, pnorm_lmax, uvnorm_lmax ) write ( *, '(a,g14.6)' ) ' Max-norm of velocity magnitude = ', uvnorm_lmax write ( *, '(a,g14.6)' ) ' Max-norm of pressure = ', pnorm_lmax call uvp_norm_l2 ( detmap, g, indx, maxelm, maxeqn, maxnp, maxquad2, & nelem, node, nquad2, phi, pnorm_l2, uvnorm_l2, wquad2 ) write ( *, '(a,g14.6)' ) ' L2-norm of velocity magnitude = ', uvnorm_l2 write ( *, '(a,g14.6)' ) ' L2-norm of pressure = ', pnorm_l2 call uvp_norm_h1 ( detmap, g, indx, maxelm, maxeqn, maxnp, maxquad2, & nelem, node, nquad2, phi, pnorm_h1, uvnorm_h1, wquad2 ) write ( *, '(a,g14.6)' ) ' H1-norm of velocity magnitude = ', uvnorm_h1 write ( *, '(a,g14.6)' ) ' H1-norm of pressure = ', pnorm_h1 ! ! Interpolate values of the pressure at nodes which do not ! have an associated pressure unknown, and create the P array. ! call press_interp ( g, indx, maxelm, maxeqn, maxnp, nelem, node, p ) ! ! Print the solution. ! iprint = 0 if ( iprint == 1 ) then call uvp_print ( g, indx, maxeqn, maxnp, np, p, xc, yc ) end if ! ! Write element and solution data. ! iplot = 1 if ( iplot == 1 ) then file_name = 'elements.txt' call element_write ( file_name, ierror, maxelm, nelem, node ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7:' write ( *, '(a)' ) ' Wrote the element data file "' // & trim ( file_name ) // '".' file_name = 'uvp.txt' call node_write ( file_name, g, ierror, indx, maxeqn, maxnp, np, p, xc, yc ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7:' write ( *, '(a)' ) ' Wrote the UVP(XY) data file "' // & trim ( file_name ) // '".' file_name = 'g.txt' call coef_write ( file_name, g, ierror, maxeqn, neqn, nx, ny, region, & nu_inv ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7:' write ( *, '(a)' ) ' Wrote the coefficient data file "' // & trim ( file_name ) // '".' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW7:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine boundary_integral ( maxside, nside, nx, ny, region, side_basis, & side_elem, side_eqn, side_etam, side_etap, side_indx, side_xsim, side_xsip ) !*****************************************************************************80 ! !! BOUNDARY_INTEGRAL sets data useful for boundary integrals. ! ! Discussion: ! ! This routine catalogs each element side that is along the boundary. ! For each side, the routine records: ! ! * the element to which the side belongs; ! * the three basis functions, in counter-clockwise order, along the side; ! * the values of ETA and XSI at the beginning and end of the side. ! * the local degree of freedom associated with the boundary condition. ! * the boundary condition being applied. ! ! Modified: ! ! 08 June 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXSIDE, the maximum number of boundary element sides. ! ! Output, integer NSIDE, the number of boundary element sides. ! ! Input, integer NX, NY, the number of elements along a horizontal, or ! vertical line. ! ! Input, character ( len = 20 ) REGION, the flow problem. ! ! Output, integer SIDE_BASIS(3,MAXSIDE), the local indices of the ! three basis functions along the boundary side. These are given in ! an order consistent with the counter clockwise ordering of the ! entire boundary. ! ! Output, integer SIDE_ELEM(MAXSIDE), the element to which the boundary ! side belongs. ! ! Output, character ( len = 3 ) SIDE_EQN(MAXSIDE), indicates the boundary ! condition being applied: ! ! '???', no condition is being applied; ! 'UI', a specified value of horizontal velocity; ! 'VI', a specified value of vertical velocity; ! 'FS' free slip, Integral (U,V) dot (boundary tangent) is penalized. ! ! Output, real ( kind = 8 ) SIDE_ETAM(MAXSIDE), SIDE_ETAP(MAXSIDE), the ! ETA value of the first and last points on the boundary side, consistent ! with the counter-clockwise ordering of the boundary. ! ! Output, integer SIDE_INDX(3,MAXSIDE), the degree of freedom ( 1, 2, or ! 3 ) associated with the local basis function, which is being controlled ! by this boundary condition. ! ! Output, real SIDE_XSIM(MAXSIDE), SIDE_XSIP(MAXSIDE), the ! XSI value of the first and last points on the boundary side, consistent ! with the counter-clockwise ordering of the boundary. ! implicit none integer maxside integer nside integer nx integer ny character ( len = 20 ) region integer side_basis(3,maxside) integer side_elem(maxside) character ( len = 3 ) side_eqn(maxside) real ( kind = 8 ) side_etam(maxside) real ( kind = 8 ) side_etap(maxside) integer side_indx(3,maxside) real ( kind = 8 ) side_xsim(maxside) real ( kind = 8 ) side_xsip(maxside) nside = 0 if ( region == 'CAVITY' ) then else if ( region == 'CHANNEL' ) then else if ( region == 'FREESLIP' ) then call boundary_integral_freeslip ( maxside, nside, nx, ny, side_basis, & side_elem, side_eqn, side_etam, side_etap, side_indx, side_xsim, & side_xsip ) else if ( region == 'STEP' ) then call boundary_integral_step ( maxside, nside, nx, ny, side_basis, & side_elem, side_eqn, side_etam, side_etap, side_indx, side_xsim, & side_xsip ) end if return end subroutine boundary_integral_freeslip ( maxside, nside, nx, ny, side_basis, & side_elem, side_eqn, side_etam, side_etap, side_indx, side_xsim, side_xsip ) !*****************************************************************************80 ! !! BOUNDARY_INTEGRAL_FREESLIP sets boundary integral data for the freeslip. ! ! Discussion: ! ! This routine catalogs each element side that is along the boundary. ! For each side, the routine records: ! ! * the element to which the side belongs; ! * the three basis functions, in counter-clockwise order, along the side; ! * the values of ETA and XSI at the beginning and end of the side. ! * the local degree of freedom associated with the boundary condition. ! * the boundary condition being applied. ! ! Modified: ! ! 08 June 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXSIDE, the maximum number of boundary element sides. ! ! Output, integer NSIDE, the number of boundary element sides. ! ! Input, integer NX, NY, the number of elements along a horizontal, or ! vertical line. ! ! Output, integer SIDE_BASIS(3,MAXSIDE), the local indices of the ! three basis functions along the boundary side. These are given in ! an order consistent with the counter clockwise ordering of the ! entire boundary. ! ! Output, integer SIDE_ELEM(MAXSIDE), the element to which the boundary ! side belongs. ! ! Output, character ( len = 3 ) SIDE_EQN(MAXSIDE), indicates the boundary ! condition being applied: ! '???', no condition is being applied; ! 'UI', a specified value of horizontal velocity; ! 'VI', a specified value of vertical velocity; ! 'FS' free slip, Integral (U,V) dot (boundary tangent) is penalized. ! ! Output, real SIDE_ETAM(MAXSIDE), SIDE_ETAP(MAXSIDE), the ! ETA value of the first and last points on the boundary side, consistent ! with the counter-clockwise ordering of the boundary. ! ! Output, integer SIDE_INDX(3,MAXSIDE), the degree of freedom ( 1, 2, or ! 3 ) associated with the local basis function, which is being controlled ! by this boundary condition. ! ! Output, real SIDE_XSIM(MAXSIDE), SIDE_XSIP(MAXSIDE), the ! XSI value of the first and last points on the boundary side, consistent ! with the counter-clockwise ordering of the boundary. ! implicit none integer maxside integer icol integer ielem integer iq integer irow integer nside integer nx integer ny integer side_basis(3,maxside) integer side_elem(maxside) character ( len = 3 ) side_eqn(maxside) real ( kind = 8 ) side_etam(maxside) real ( kind = 8 ) side_etap(maxside) integer side_indx(3,maxside) real ( kind = 8 ) side_xsim(maxside) real ( kind = 8 ) side_xsip(maxside) ! ! Bottom boundary. ! Freeslip. ! do icol = 1, nx ielem = 2 * ( icol - 1 ) * ny + 1 nside = nside + 1 side_eqn(nside) = 'FS' side_elem(nside) = ielem iq = 3 side_basis(1,nside) = iq side_indx(1,nside) = 2 iq = 4 side_basis(2,nside) = iq side_indx(2,nside) = 2 iq = 1 side_basis(3,nside) = iq side_indx(3,nside) = 2 side_etam(nside) = 0.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 1.0D+00 end do ! ! Right boundary. ! Vertical velocity is 0. ! do irow = 1, ny ielem = 2 * ny * ( nx - 1 ) + 2 * irow nside = nside + 1 side_eqn(nside) = 'V0' side_elem(nside) = ielem iq = 2 side_basis(1,nside) = iq side_indx(1,nside) = 2 iq = 6 side_basis(2,nside) = iq side_indx(2,nside) = 2 iq = 3 side_basis(3,nside) = iq side_indx(3,nside) = 2 side_etam(nside) = 1.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 0.0D+00 end do ! ! Top boundary. ! Free slip ! do icol = 1, nx ielem = 2 * icol * ny nside = nside + 1 side_eqn(nside) = 'FS' side_elem(nside) = ielem iq = 3 side_basis(1,nside) = iq side_indx(1,nside) = 2 iq = 4 side_basis(2,nside) = iq side_indx(2,nside) = 2 iq = 1 side_basis(3,nside) = iq side_indx(3,nside) = 2 side_etam(nside) = 0.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 1.0D+00 end do ! ! Left boundary. ! Horizontal velocity is specified. ! Vertical velocity is 0. ! do irow = 1, ny ielem = 2 * irow - 1 nside = nside + 1 side_eqn(nside) = 'UI' side_elem(nside) = ielem iq = 2 side_basis(1,nside) = iq side_indx(1,nside) = 1 iq = 6 side_basis(2,nside) = iq side_indx(2,nside) = 1 iq = 3 side_basis(3,nside) = iq side_indx(3,nside) = 1 side_etam(nside) = 1.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 0.0D+00 nside = nside + 1 side_eqn(nside) = 'V0' side_elem(nside) = ielem iq = 2 side_basis(1,nside) = iq side_indx(1,nside) = 2 iq = 6 side_basis(2,nside) = iq side_indx(2,nside) = 2 iq = 3 side_basis(3,nside) = iq side_indx(3,nside) = 2 side_etam(nside) = 1.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 0.0D+00 end do return end subroutine boundary_integral_step ( maxside, nside, nx, ny, side_basis, & side_elem, side_eqn, side_etam, side_etap, side_indx, side_xsim, side_xsip ) !*****************************************************************************80 ! !! BOUNDARY_INTEGRAL_STEP sets boundary integral data for the step. ! ! Discussion: ! ! This routine catalogs each element side that is along the boundary. ! For each side, the routine records: ! ! * the element to which the side belongs; ! * the three basis functions, in counter-clockwise order, along the side; ! * the values of ETA and XSI at the beginning and end of the side. ! * the local degree of freedom associated with the boundary condition. ! * the boundary condition being applied. ! ! Modified: ! ! 08 June 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXSIDE, the maximum number of boundary element sides. ! ! Output, integer NSIDE, the number of boundary element sides. ! ! Input, integer NX, NY, the number of elements along a horizontal, or ! vertical line. ! ! Output, integer SIDE_BASIS(3,MAXSIDE), the local indices of the ! three basis functions along the boundary side. These are given in ! an order consistent with the counter clockwise ordering of the ! entire boundary. ! ! Output, integer SIDE_ELEM(MAXSIDE), the element to which the boundary ! side belongs. ! ! Output, character ( len = 3 ) SIDE_EQN(MAXSIDE), indicates the boundary ! condition being applied: ! '???', no condition is being applied; ! 'UI', a specified value of horizontal velocity; ! 'VI', a specified value of vertical velocity; ! 'FS' free slip, Integral (U,V) dot (boundary tangent) is penalized. ! ! Output, real SIDE_ETAM(MAXSIDE), SIDE_ETAP(MAXSIDE), the ! ETA value of the first and last points on the boundary side, consistent ! with the counter-clockwise ordering of the boundary. ! ! Output, integer SIDE_INDX(3,MAXSIDE), the degree of freedom ( 1, 2, or ! 3 ) associated with the local basis function, which is being controlled ! by this boundary condition. ! ! Output, real SIDE_XSIM(MAXSIDE), SIDE_XSIP(MAXSIDE), the ! XSI value of the first and last points on the boundary side, consistent ! with the counter-clockwise ordering of the boundary. ! implicit none integer maxside integer ielem integer ihi integer ilo integer iq integer nside integer nx integer nx1 integer nx2 integer nx3 integer ny integer ny1 integer ny2 integer side_basis(3,maxside) integer side_elem(maxside) character ( len = 3 ) side_eqn(maxside) real ( kind = 8 ) side_etam(maxside) real ( kind = 8 ) side_etap(maxside) integer side_indx(3,maxside) real ( kind = 8 ) side_xsim(maxside) real ( kind = 8 ) side_xsip(maxside) nside = 0 nx1 = nint ( 0.375D+00 * real ( nx, kind = 8 ) ) nx1 = max ( 1, nx1 ) nx2 = nint ( 0.125D+00 * real ( nx, kind = 8 ) ) nx2 = max ( 1, nx2 ) nx3 = nx - nx1 - nx2 ny1 = max ( 1, nint ( 0.500D+00 * real ( ny, kind = 8 ) ) ) ny2 = ny - ny1 ! ! Segment 1: horizontal bottom near inflow. ! ilo = 1 ihi = 2 * ( nx1 - 1 ) * ny + 1 do ielem = ilo, ihi, 2 * ny nside = nside + 1 side_eqn(nside) = 'FS' side_elem(nside) = ielem iq = 3 side_basis(1,nside) = iq side_indx(1,nside) = 2 iq = 4 side_basis(2,nside) = iq side_indx(2,nside) = 2 iq = 1 side_basis(3,nside) = iq side_indx(3,nside) = 2 side_etam(nside) = 0.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 1.0D+00 end do ! ! Segment 2: vertical riser of step. ! ilo = 2 * ( nx1 - 1 ) * ny + 2 ihi = 2 * ( nx1 - 1 ) * ny + 2 * ny1 do ielem = ilo, ihi, 2 nside = nside + 1 side_eqn(nside) = 'FS' side_elem(nside) = ielem iq = 2 side_basis(1,nside) = iq side_indx(1,nside) = 2 iq = 6 side_basis(2,nside) = iq side_indx(2,nside) = 2 iq = 3 side_basis(3,nside) = iq side_indx(3,nside) = 2 side_etam(nside) = 1.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 0.0D+00 end do ! ! Segment 3: horizontal step. ! ilo = 2 * nx1 * ny + 1 ihi = 2 * nx1 * ny + 1 + 2 * ny2 * ( nx2 - 1 ) do ielem = ilo, ihi, 2 * ny2 nside = nside + 1 side_eqn(nside) = 'FS' side_elem(nside) = ielem iq = 3 side_basis(1,nside) = iq side_indx(1,nside) = 2 iq = 4 side_basis(2,nside) = iq side_indx(2,nside) = 2 iq = 1 side_basis(3,nside) = iq side_indx(3,nside) = 2 side_etam(nside) = 0.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 1.0D+00 end do ! ! Segment 4: vertical descender of step. ! ilo = 2 * nx1 * ny + 2 * nx2 * ny2 + 2 * ( ny1 - 1 ) + 1 ihi = 2 * nx1 * ny + 2 * nx2 * ny2 + 1 do ielem = ilo, ihi, -2 nside = nside + 1 side_eqn(nside) = 'FS' side_elem(nside) = ielem iq = 2 side_basis(1,nside) = iq side_indx(1,nside) = 2 iq = 6 side_basis(2,nside) = iq side_indx(2,nside) = 2 iq = 3 side_basis(3,nside) = iq side_indx(3,nside) = 2 side_etam(nside) = 1.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 0.0D+00 end do ! ! Segment 5: horizontal bottom near outflow ! ilo = 2 * nx1 * ny + 2 * nx2 * ny2 + 1 ihi = 2 * nx1 * ny + 2 * nx2 * ny2 + 2 * ( nx3 - 1 ) * ny + 1 do ielem = ilo, ihi, 2 * ny nside = nside + 1 side_eqn(nside) = 'FS' side_elem(nside) = ielem iq = 3 side_basis(1,nside) = iq side_indx(1,nside) = 2 iq = 4 side_basis(2,nside) = iq side_indx(2,nside) = 2 iq = 1 side_basis(3,nside) = iq side_indx(3,nside) = 2 side_etam(nside) = 0.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 1.0D+00 end do ! ! Segment 6: top horizontal near outflow ! ilo = 2 * ( nx1 * ( ny1 + ny2 ) + nx2 * ny2 + nx3 * ( ny1 + ny2 ) ) ihi = 2 * ( nx1 * ( ny1 + ny2 ) + nx2 * ny2 + ny1 + ny2 ) do ielem = ilo, ihi, - 2 * ny nside = nside + 1 side_eqn(nside) = 'FS' side_elem(nside) = ielem iq = 3 side_basis(1,nside) = iq side_indx(1,nside) = 2 iq = 4 side_basis(2,nside) = iq side_indx(2,nside) = 2 iq = 1 side_basis(3,nside) = iq side_indx(3,nside) = 2 side_etam(nside) = 0.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 1.0D+00 end do ! ! Segment 7: top horizontal over step ! ilo = 2 * nx1 * ny + 2 * nx2 * ny2 ihi = 2 * nx1 * ny + 2 * ny2 do ielem = ilo, ihi, - 2 * ny2 nside = nside + 1 side_eqn(nside) = 'FS' side_elem(nside) = ielem iq = 3 side_basis(1,nside) = iq side_indx(1,nside) = 2 iq = 4 side_basis(2,nside) = iq side_indx(2,nside) = 2 iq = 1 side_basis(3,nside) = iq side_indx(3,nside) = 2 side_etam(nside) = 0.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 1.0D+00 end do ! ! Segment 8: top horizontal near inflow ! ilo = 2 * ny * nx1 ihi = 2 * ny do ielem = ilo, ihi, - 2 * ny nside = nside + 1 side_eqn(nside) = 'FS' side_elem(nside) = ielem iq = 3 side_basis(1,nside) = iq side_indx(1,nside) = 2 iq = 4 side_basis(2,nside) = iq side_indx(2,nside) = 2 iq = 1 side_basis(3,nside) = iq side_indx(3,nside) = 2 side_etam(nside) = 0.0D+00 side_etap(nside) = 0.0D+00 side_xsim(nside) = 0.0D+00 side_xsip(nside) = 1.0D+00 end do return end subroutine boundary_integral_print ( maxside, nside, side_basis, side_elem, & side_eqn, side_etam, side_etap, side_indx, side_xsim, side_xsip ) !*****************************************************************************80 ! !! BOUNDARY_INTEGRAL_PRINT prints data useful for boundary integrals. ! ! Modified: ! ! 12 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXSIDE, the maximum number of boundary element sides. ! ! Input, integer NSIDE, the number of boundary element sides. ! ! Input, integer SIDE_BASIS(3,MAXSIDE), the local indices of the ! three basis functions along the boundary side. These are given in ! an order consistent with the counter clockwise ordering of the ! entire boundary. ! ! Input, integer SIDE_ELEM(MAXSIDE), the element to which the boundary ! side belongs. ! ! Input, character ( len = 3 ) SIDE_EQN(MAXSIDE), indicates the boundary ! condition being applied: ! ! '???', no condition is being applied; ! 'UI', a specified value of horizontal velocity; ! 'VI', a specified value of vertical velocity; ! 'FS', a no slip condition, applied to velocity. ! ! Input, real SIDE_ETAM(MAXSIDE), SIDE_ETAP(MAXSIDE), the ! ETA value of the first and last points on the boundary side, consistent ! with the counter-clockwise ordering of the boundary. ! ! Input, integer SIDE_INDX(3,MAXSIDE), the degree of freedom ( 1, 2, or ! 3 ) associated with the local basis function, which is being controlled ! by this boundary condition. ! ! Input, real SIDE_XSIM(MAXSIDE), SIDE_XSIP(MAXSIDE), the ! XSI value of the first and last points on the boundary side, consistent ! with the counter-clockwise ordering of the boundary. ! implicit none integer maxside integer iside integer nside integer side_basis(3,maxside) integer side_elem(maxside) character ( len = 3 ) side_eqn(maxside) real ( kind = 8 ) side_etam(maxside) real ( kind = 8 ) side_etap(maxside) integer side_indx(3,maxside) real ( kind = 8 ) side_xsim(maxside) real ( kind = 8 ) side_xsip(maxside) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BOUNDARY_INTEGRAL_PRINT:' write ( *, '(a,i8)' ) ' Number of boundary sides is ', nside write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' ISIDE ELEMENT EQN_TYPE' write ( *, '(a)' ) ' IQ1, INDX1, XSI, ETA' write ( *, '(a)' ) ' IQ2, INDX2' write ( *, '(a)' ) ' IQ3, INDX3, XSI, ETA' write ( *, '(a)' ) ' ' do iside = 1, nside write ( *, '(a)' ) ' ' write ( *, '(i8,3x,i8,3x,a3)' ) iside, side_elem(iside), side_eqn(iside) write ( *, '(6x,i8,3x,i8,3x,2g14.6)' ) side_basis(1,iside), & side_indx(1,iside), side_xsim(iside), side_etam(iside) write ( *, '(6x,i8,3x,i8,3x,2g14.6)' ) side_basis(2,iside), & side_indx(2,iside) write ( *, '(6x,i8,3x,i8,3x,2g14.6)' ) side_basis(3,iside), & side_indx(3,iside), side_xsip(iside), side_etap(iside) end do return end subroutine boundary_shape ( bc_tag, normal, region, tangent ) !*****************************************************************************80 ! !! BOUNDARY_SHAPE returns the tangent and normal vectors along the boundary. ! ! Modified: ! ! 08 June 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer BC_TAG, a code which is nonzero at every boundary node, ! and which indicates the type of boundary conditions at that node. ! The code used can vary from problem to problem. ! ! Output, real NORMAL(2), the X and Y components of the normal ! vector to the boundary at the given point. The normal vector should ! point out of the region. ! ! Input, character ( len = 20 ) REGION, the flow problem. ! ! Output, real TANGENT(2), the X and Y components of the ! tangential vector to the boundary, at the given point. The tangent ! vector should point in the counter-clockwise direction, relative to ! a point within the interior of the region. ! implicit none integer bc_tag real ( kind = 8 ) normal(2) character ( len = * ) region real ( kind = 8 ) tangent(2) if ( region == 'CAVITY' ) then call boundary_shape_cavity ( bc_tag, normal, tangent ) else if ( region == 'CHANNEL' ) then call boundary_shape_channel ( bc_tag, normal, tangent ) else if ( region == 'FREESLIP' ) then call boundary_shape_freeslip ( bc_tag, normal, tangent ) else if ( region == 'STEP' ) then call boundary_shape_step ( bc_tag, normal, tangent ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BOUNDARY_SHAPE - Fatal error!' write ( *, '(a)' ) ' Unrecognized region!' stop end if return end subroutine boundary_shape_cavity ( bc_tag, normal, tangent ) !*****************************************************************************80 ! !! BOUNDARY_SHAPE_CAVITY returns the boundary of the cavity problem. ! ! Discussion: ! ! Here is a schematic of the values of BC_TAG along the boundary ! for this problem: ! ! 7 6 6 6 6 6 5 ! 8 0 0 0 0 0 4 ! 8 0 0 0 0 0 4 ! 8 0 0 0 0 0 4 ! 1 2 2 2 2 2 3 ! ! Modified: ! ! 08 June 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer BC_TAG, a code which is nonzero at every boundary node, ! and which indicates the type of boundary conditions at that node. ! The code used can vary from problem to problem. ! ! Output, real NORMAL(2), the X and Y components of the normal ! vector to the boundary at the given point. The normal vector should ! point out of the region. ! ! Output, real TANGENT(2), the X and Y components of the ! tangential vector to the boundary, at the given point. The tangent ! vector should point in the counter-clockwise direction, relative to ! a point within the interior of the region. ! implicit none integer bc_tag real ( kind = 8 ) normal(2) real ( kind = 8 ) tangent(2) if ( bc_tag == 1 ) then normal(1) = - sqrt ( 2.0D+00 ) / 2.0D+00 normal(2) = - sqrt ( 2.0D+00 ) / 2.0D+00 tangent(1) = sqrt ( 2.0D+00 ) / 2.0D+00 tangent(2) = - sqrt ( 2.0D+00 ) / 2.0D+00 else if ( bc_tag == 2 ) then normal(1) = 0.0D+00 normal(2) = - 1.0D+00 tangent(1) = 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 3 ) then normal(1) = sqrt ( 2.0D+00 ) / 2.0D+00 normal(2) = - sqrt ( 2.0D+00 ) / 2.0D+00 tangent(1) = sqrt ( 2.0D+00 ) / 2.0D+00 tangent(2) = sqrt ( 2.0D+00 ) / 2.0D+00 else if ( bc_tag == 4 ) then normal(1) = 1.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = 1.0D+00 else if ( bc_tag == 5 ) then normal(1) = sqrt ( 2.0D+00 ) / 2.0D+00 normal(2) = sqrt ( 2.0D+00 ) / 2.0D+00 tangent(1) = - sqrt ( 2.0D+00 ) / 2.0D+00 tangent(2) = sqrt ( 2.0D+00 ) / 2.0D+00 else if ( bc_tag == 6 ) then normal(1) = 0.0D+00 normal(2) = 1.0D+00 tangent(1) = - 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 7 ) then normal(1) = - sqrt ( 2.0D+00 ) / 2.0D+00 normal(2) = sqrt ( 2.0D+00 ) / 2.0D+00 tangent(1) = - sqrt ( 2.0D+00 ) / 2.0D+00 tangent(2) = - sqrt ( 2.0D+00 ) / 2.0D+00 else if ( bc_tag == 8 ) then normal(1) = - 1.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = - 1.0D+00 else normal(1) = 0.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = 0.0D+00 end if return end subroutine boundary_shape_channel ( bc_tag, normal, tangent ) !*****************************************************************************80 ! !! BOUNDARY_SHAPE_CHANNEL returns the boundary of the channel problem. ! ! Discussion: ! ! Here is a schematic of the values of BC_TAG along the boundary ! for this problem: ! ! 7 6 6 6 6 6 5 ! 8 0 0 0 0 0 4 ! 8 0 0 0 0 0 4 ! 8 0 0 0 0 0 4 ! 1 2 2 2 2 2 3 ! ! Modified: ! ! 08 June 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer BC_TAG, a code which is nonzero at every boundary node, ! and which indicates the type of boundary conditions at that node. ! The code used can vary from problem to problem. ! ! Output, real NORMAL(2), the X and Y components of the normal ! vector to the boundary at the given point. The normal vector should ! point out of the region. ! ! Output, real TANGENT(2), the X and Y components of the ! tangential vector to the boundary, at the given point. The tangent ! vector should point in the counter-clockwise direction, relative to ! a point within the interior of the region. ! implicit none integer bc_tag real ( kind = 8 ) normal(2) real ( kind = 8 ) tangent(2) if ( bc_tag == 1 ) then normal(1) = 0.0D+00 normal(2) = - 1.0D+00 tangent(1) = 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 2 ) then normal(1) = 0.0D+00 normal(2) = - 1.0D+00 tangent(1) = 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 3 ) then normal(1) = 0.0D+00 normal(2) = - 1.0D+00 tangent(1) = 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 4 ) then normal(1) = 1.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = 1.0D+00 else if ( bc_tag == 5 ) then normal(1) = 0.0D+00 normal(2) = 1.0D+00 tangent(1) = - 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 6 ) then normal(1) = 0.0D+00 normal(2) = 1.0D+00 tangent(1) = - 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 7 ) then normal(1) = 0.0D+00 normal(2) = 1.0D+00 tangent(1) = - 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 8 ) then normal(1) = - 1.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = - 1.0D+00 else normal(1) = 0.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = 0.0D+00 end if return end subroutine boundary_shape_freeslip ( bc_tag, normal, tangent ) !*****************************************************************************80 ! !! BOUNDARY_SHAPE_FREESLIP returns the boundary of the freeslip problem. ! ! Discussion: ! ! 12 March 1999: ! I tried using the "diagonal" normals and tangents in the corners, ! but that just caused the code to compute nonzero velocities there ! that didn't make sense. So I'm just extending the flatness. ! ! Here is a schematic of the values of BC_TAG along the boundary ! for this problem: ! ! 7 6 6 6 6 6 5 ! 8 0 0 0 0 0 4 ! 8 0 0 0 0 0 4 ! 8 0 0 0 0 0 4 ! 1 2 2 2 2 2 3 ! ! Modified: ! ! 08 June 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer BC_TAG, a code which is nonzero at every boundary node, ! and which indicates the type of boundary conditions at that node. ! The code used can vary from problem to problem. ! ! Output, real NORMAL(2), the X and Y components of the normal ! vector to the boundary at the given point. The normal vector should ! point out of the region. ! ! Output, real TANGENT(2), the X and Y components of the ! tangential vector to the boundary, at the given point. The tangent ! vector should point in the counter-clockwise direction, relative to ! a point within the interior of the region. ! implicit none integer bc_tag real ( kind = 8 ) normal(2) real ( kind = 8 ) tangent(2) if ( bc_tag == 1 ) then normal(1) = 0.0D+00 normal(2) = - 1.0D+00 tangent(1) = 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 2 ) then normal(1) = 0.0D+00 normal(2) = - 1.0D+00 tangent(1) = 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 3 ) then normal(1) = 0.0D+00 normal(2) = - 1.0D+00 tangent(1) = 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 4 ) then normal(1) = 1.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = 1.0D+00 else if ( bc_tag == 5 ) then normal(1) = 0.0D+00 normal(2) = 1.0D+00 tangent(1) = - 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 6 ) then normal(1) = 0.0D+00 normal(2) = 1.0D+00 tangent(1) = - 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 7 ) then normal(1) = 0.0D+00 normal(2) = 1.0D+00 tangent(1) = - 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 8 ) then normal(1) = - 1.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = - 1.0D+00 else normal(1) = 0.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = 0.0D+00 end if return end subroutine boundary_shape_step ( bc_tag, normal, tangent ) !*****************************************************************************80 ! !! BOUNDARY_SHAPE_STEP returns the boundary of the step problem. ! ! Discussion: ! ! Here is a schematic of the values of BC_TAG along the boundary ! for this problem: ! ! 15 14 14 14 14 14 14 14 14 13 ! 16 0 0 0 0 0 0 0 0 12 ! 16 0 0 0 0 0 0 0 0 12 ! 16 0 0 0 0 0 0 0 0 12 ! 16 0 0 5 6 6 7 0 0 12 ! 16 0 0 4 -1 -1 8 0 0 12 ! 16 0 0 4 -1 -1 8 0 0 12 ! 1 2 2 3 -1 -1 9 10 10 11 ! ! Modified: ! ! 08 June 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer BC_TAG, a code which is nonzero at every boundary node, ! and which indicates the type of boundary conditions at that node. ! The code used can vary from problem to problem. ! ! Output, real NORMAL(2), the X and Y components of the normal ! vector to the boundary at the given point. The normal vector should ! point out of the region. ! ! Output, real TANGENT(2), the X and Y components of the ! tangential vector to the boundary, at the given point. The tangent ! vector should point in the counter-clockwise direction, relative to ! a point within the interior of the region. ! implicit none integer bc_tag real ( kind = 8 ) normal(2) real ( kind = 8 ) tangent(2) if ( bc_tag == 1 ) then normal(1) = 0.0D+00 normal(2) = - 1.0D+00 tangent(1) = 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 2 ) then normal(1) = 0.0D+00 normal(2) = - 1.0D+00 tangent(1) = 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 3 ) then normal(1) = sqrt ( 2.0D+00 ) / 2.0D+00 normal(2) = - sqrt ( 2.0D+00 ) / 2.0D+00 tangent(1) = sqrt ( 2.0D+00 ) / 2.0D+00 tangent(2) = sqrt ( 2.0D+00 ) / 2.0D+00 else if ( bc_tag == 4 ) then normal(1) = 1.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = 1.0D+00 else if ( bc_tag == 5 ) then normal(1) = sqrt ( 2.0D+00 ) / 2.0D+00 normal(2) = - sqrt ( 2.0D+00 ) / 2.0D+00 tangent(1) = sqrt ( 2.0D+00 ) / 2.0D+00 tangent(2) = sqrt ( 2.0D+00 ) / 2.0D+00 else if ( bc_tag == 6 ) then normal(1) = 0.0D+00 normal(2) = - 1.0D+00 tangent(1) = 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 7 ) then normal(1) = - sqrt ( 2.0D+00 ) / 2.0D+00 normal(2) = - sqrt ( 2.0D+00 ) / 2.0D+00 tangent(1) = sqrt ( 2.0D+00 ) / 2.0D+00 tangent(2) = - sqrt ( 2.0D+00 ) / 2.0D+00 else if ( bc_tag == 8 ) then normal(1) = - 1.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = - 1.0D+00 else if ( bc_tag == 9 ) then normal(1) = - sqrt ( 2.0D+00 ) / 2.0D+00 normal(2) = - sqrt ( 2.0D+00 ) / 2.0D+00 tangent(1) = sqrt ( 2.0D+00 ) / 2.0D+00 tangent(2) = - sqrt ( 2.0D+00 ) / 2.0D+00 else if ( bc_tag == 10 ) then normal(1) = 0.0D+00 normal(2) = - 1.0D+00 tangent(1) = 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 11 ) then normal(1) = 0.0D+00 normal(2) = - 1.0D+00 tangent(1) = 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 12 ) then normal(1) = 1.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = 1.0D+00 else if ( bc_tag == 13 ) then normal(1) = 0.0D+00 normal(2) = 1.0D+00 tangent(1) = - 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 14 ) then normal(1) = 0.0D+00 normal(2) = 1.0D+00 tangent(1) = - 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 15 ) then normal(1) = 0.0D+00 normal(2) = 1.0D+00 tangent(1) = - 1.0D+00 tangent(2) = 0.0D+00 else if ( bc_tag == 16 ) then normal(1) = - 1.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = - 1.0D+00 else normal(1) = 0.0D+00 normal(2) = 0.0D+00 tangent(1) = 0.0D+00 tangent(2) = 0.0D+00 end if return end subroutine boundary_value ( bc_tag, pbc, region, region_ymax, ubc, vbc, y ) !*****************************************************************************80 ! !! BOUNDARY_VALUE evaluates boundary condition functions. ! ! Modified: ! ! 11 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer BC_TAG, the value of the boundary condition "tag" ! for the node. The meaning of this tag depends on the problem ! being solved. ! ! Output, real PBC, the value specified for pressure, at ! the point (X,Y), for the problem specified by REGION. If no boundary ! condition is imposed for the particular variable, location and ! problem, simply set the quantity to zero. ! ! Input, character ( len = 20 ) REGION, the flow problem. ! ! Output, real UBC, VBC, the values specified ! for horizontal and vertical velocity, at the point ! (X,Y), for the problem specified by REGION. If no boundary ! condition is imposed for the particular variable, location and ! problem, simply set the quantity to zero. ! ! Input, real Y, the coordinates of the node. ! implicit none integer bc_tag real ( kind = 8 ) pbc character ( len = * ) region real ( kind = 8 ) region_ymax real ( kind = 8 ) ubc real ( kind = 8 ) vbc real ( kind = 8 ) y pbc = 0.0D+00 ubc = 0.0D+00 vbc = 0.0D+00 ! ! CAVITY ! ! 7 6 6 6 6 6 5 ! 8 0 0 0 0 0 4 ! 8 0 0 0 0 0 4 ! 8 0 0 0 0 0 4 ! 1 2 2 2 2 2 3 ! if ( region == 'CAVITY' ) then if ( bc_tag == 6 ) then ubc = 1.0D+00 end if ! ! CHANNEL ! ! 7 6 6 6 6 6 5 ! 8 0 0 0 0 0 4 ! 8 0 0 0 0 0 4 ! 8 0 0 0 0 0 4 ! 1 2 2 2 2 2 3 ! else if ( region == 'CHANNEL' ) then if ( bc_tag == 4 ) then vbc = 0.0D+00 else if ( bc_tag == 5 ) then pbc = 100.0D+00 else if ( bc_tag == 8 ) then ubc = y * ( region_ymax - y ) vbc = 0.0D+00 end if ! ! FREESLIP ! ! 7 6 6 6 6 6 5 ! 8 0 0 0 0 0 4 ! 8 0 0 0 0 0 4 ! 8 0 0 0 0 0 4 ! 1 2 2 2 2 2 3 ! else if ( region == 'FREESLIP' ) then if ( bc_tag == 4 ) then vbc = 0.0D+00 else if ( bc_tag == 8 ) then ubc = y * ( region_ymax - y ) vbc = 0.0D+00 end if ! ! STEP ! ! 15 14 14 14 14 14 14 14 14 13 ! 16 0 0 0 0 0 0 0 0 12 ! 16 0 0 0 0 0 0 0 0 12 ! 16 0 0 0 0 0 0 0 0 12 ! 16 0 0 5 6 6 7 0 0 12 ! 16 0 0 4 -1 -1 8 0 0 12 ! 16 0 0 4 -1 -1 8 0 0 12 ! 1 2 2 3 -1 -1 9 10 10 11 ! else if ( region == 'STEP' ) then if ( bc_tag == 12 ) then vbc = 0.0D+00 else if ( bc_tag == 16 ) then ubc = y * ( region_ymax - y ) vbc = 0.0D+00 end if else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BOUNDARY_VALUE - Fatal error!' write ( *, '(a)' ) ' Unrecognized region!' stop end if return end subroutine ch_cap ( c ) !*****************************************************************************80 ! !! CH_CAP capitalizes a single character. ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! implicit none character c integer itemp itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end subroutine coef_read ( file_name, g, ierror, maxeqn, neqn, nx, ny, region, & nu_inv ) !*****************************************************************************80 ! !! COEF_READ reads the coefficient data from a file. ! ! Modified: ! ! 29 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the file to read. ! ! Output, real G(MAXEQN), the current solution vector. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer MAXEQN, the maximum number of equations. ! ! Output, integer NEQN, the number of equations. ! ! Output, integer NX, NY, the number of elements along a horizontal, or ! vertical line. ! ! Output, character ( len = 20 ) REGION, the flow problem. ! ! Output, real NU_INV, the inverse dynamic viscosity. ! implicit none integer maxeqn character ( len = * ) file_name real ( kind = 8 ) g(maxeqn) integer i integer ierror integer ios integer iunit integer neqn integer nx integer ny character ( len = * ) region real ( kind = 8 ) nu_inv ierror = 0 neqn = 0 nx = 0 ny = 0 region = ' ' nu_inv = 0.0D+00 call get_unit ( iunit ) open ( unit = iunit, file = file_name, status = 'old', & iostat = ios ) if ( ios /= 0 ) then ierror = ios write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COEF_READ - Warning!' write ( *, '(a)' ) ' Could not open the file:' write ( *, '(a)' ) trim ( file_name ) return end if do i = 1, 10 read ( iunit, '(1x)', iostat = ios ) if ( ios /= 0 ) then ierror = ios write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COEF_READ - Warning!' write ( *, '(a)' ) ' Read error!' return end if end do read ( iunit, *, iostat = ios ) neqn if ( ios /= 0 ) then ierror = ios write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COEF_READ - Warning!' write ( *, '(a)' ) ' Read error!' return end if read ( iunit, *, iostat = ios ) nx if ( ios /= 0 ) then ierror = ios write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COEF_READ - Warning!' write ( *, '(a)' ) ' Read error!' return end if read ( iunit, *, iostat = ios ) ny if ( ios /= 0 ) then ierror = ios write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COEF_READ - Warning!' write ( *, '(a)' ) ' Read error!' return end if read ( iunit, '(a)', iostat = ios ) region if ( ios /= 0 ) then ierror = ios write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COEF_READ - Warning!' write ( *, '(a)' ) ' Read error!' return end if read ( iunit, *, iostat = ios ) nu_inv if ( ios /= 0 ) then ierror = ios write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COEF_READ - Warning!' write ( *, '(a)' ) ' Read error!' return end if read ( iunit, '(1x)', iostat = ios ) if ( ios /= 0 ) then ierror = ios write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COEF_READ - Warning!' write ( *, '(a)' ) ' Read error!' return end if if ( maxeqn < neqn ) then ierror = 3 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COEF_READ - Warning!' write ( *, '(a,i8)' ) ' Input NEQN = ', neqn write ( *, '(a,i8)' ) ' Internal MAXEQN = ', maxeqn close ( unit = iunit ) return end if do i = 1, neqn read ( iunit, *, iostat = ios ) g(i) if ( ios /= 0 ) then ierror = ios write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COEF_READ - Warning!' write ( *, '(a)' ) ' Read error!' return end if end do close ( unit = iunit ) return end subroutine coef_write ( file_name, g, ierror, maxeqn, neqn, nx, ny, region, & nu_inv ) !*****************************************************************************80 ! !! COEF_WRITE writes the coefficient data to a file for possible restart. ! ! Modified: ! ! 29 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the file to write. ! ! Input, real G(MAXEQN), the current solution vector. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer MAXEQN, the maximum number of equations. ! ! Input, integer NEQN, the number of equations. ! ! Input, integer NX, NY, the number of elements along a horizontal, or ! vertical line. ! ! Input, character ( len = 20 ) REGION, the flow problem. ! ! Input, real NU_INV, the inverse dynamic viscosity. ! implicit none integer maxeqn character ( len = * ) file_name real ( kind = 8 ) g(maxeqn) integer i integer ierror integer ios integer iunit integer neqn integer nx integer ny character ( len = * ) region real ( kind = 8 ) nu_inv ierror = 0 call get_unit ( iunit ) open ( unit = iunit, file = file_name, status = 'replace', & iostat = ios ) if ( ios /= 0 ) then ierror = ios write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COEF_WRITE - Warning!' write ( *, '(a)' ) ' Could not open the file.' return end if write ( iunit, '(a)' ) '# g.txt FLOW7 solution vector' write ( iunit, '(a)' ) '# ' write ( iunit, '(a)' ) '# NEQN' write ( iunit, '(a)' ) '# NX' write ( iunit, '(a)' ) '# NY' write ( iunit, '(a)' ) '# REGION' write ( iunit, '(a)' ) '# NU_INV' write ( iunit, '(a)' ) '# ' write ( iunit, '(a)' ) '# (G(I),I=1,NEQN)' write ( iunit, '(a)' ) '#' write ( iunit, '(i8)' ) neqn write ( iunit, '(i8)' ) nx write ( iunit, '(i8)' ) ny write ( iunit, '(a)' ) trim ( region ) write ( iunit, '(g14.6)' ) nu_inv write ( iunit, '(1x)' ) do i = 1, neqn write ( iunit, '(g14.6)' ) g(i) end do close ( unit = iunit ) return end subroutine element_print ( maxelm, nelem, node ) !*****************************************************************************80 ! !! ELEMENT_PRINT prints out the elements. ! ! Modified: ! ! 07 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(6,MAXELM), the nodes that make up each element. ! implicit none integer maxelm integer i integer j integer nelem integer node(6,maxelm) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ELEMENT_PRINT -' write ( *, '(a)' ) ' Decomposition of region into elements:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Element Nodes' write ( *, '(a)' ) ' ' do i = 1, nelem write ( *, '(i8,6x,6i8)' ) i, ( node(j,i), j = 1, 6 ) end do return end subroutine element_write ( file_name, ierror, maxelm, nelem, node ) !*****************************************************************************80 ! !! ELEMENT_WRITE writes element information to a file. ! ! Discussion: ! ! The data may be used to plot the finite element mesh and solution. ! ! Modified: ! ! 29 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(6,MAXELM), the nodes that make up each element. ! implicit none integer maxelm character ( len = * ) file_name integer i integer ielem integer ierror integer ios integer iunit integer nelem integer node(6,maxelm) integer npe ierror = 0 call get_unit ( iunit ) open ( unit = iunit, file = file_name, status = 'replace', & iostat = ios ) if ( ios /= 0 ) then ierror = ios write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ELEMENT_WRITE - Warning!' write ( *, '(a)' ) ' Could not open the file!' return end if npe = 6 write ( iunit, * ) nelem write ( iunit, * ) npe do ielem = 1, nelem write ( iunit, * ) node(1:npe,ielem) end do close ( unit = iunit ) return end subroutine eqn_nabor_print ( eqn_nabor, eqn_nabor_max, eqn_nabor_num, neqn ) !*****************************************************************************80 ! !! EQN_NABOR_PRINT prints the equation neighbor array. ! ! Modified: ! ! 18 May 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer EQN_NABOR(EQN_NABOR_MAX,NEQN), the equation neigbors. ! ! Input, integer EQN_NABOR_MAX, the maximum number of equation neighbors. ! ! Input, integer EQN_NABOR_NUM(NEQN), the number of equation neighbors. ! ! Input, integer NEQN, the number of equations and variables. ! implicit none integer neqn integer eqn_nabor_max integer i integer j integer jhi integer jlo integer eqn_nabor(eqn_nabor_max,neqn) integer eqn_nabor_num(neqn) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Equation Neighbor Array' write ( *, '(a)' ) ' ' do i = 1, neqn do jlo = 1, eqn_nabor_num(i), 19 jhi = min ( jlo + 18, eqn_nabor_num(i) ) if ( jlo == 1 ) then write ( *, '(i4,(19i4))' ) i, ( eqn_nabor(j,i), j = jlo, jhi ) else write ( *, '(4x,(19i4))' ) ( eqn_nabor(j,i), j = jlo, jhi ) end if end do end do return end subroutine eqn_nabor_set ( eqn_nabor, eqn_nabor_max, eqn_nabor_num, indx, & neqn, np, node_nabor, node_nabor_max, node_nabor_num ) !*****************************************************************************80 ! !! EQN_NABOR_SET sets up the equation neighbor array. ! ! Discussion: ! ! The equation neighbor array allows us to quickly find out, for each ! equation or variable, the equations or variables that are its "neighbors". ! Equations I and J are neighbors if variable I influences equation J or ! variable J influences equation I. ! ! For our grid, there should be at most 45 equations or variables in ! each neighborhood. ! ! Modified: ! ! 17 May 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NP, the number of nodes. ! ! Input, integer NODE_NABOR(NABOR_MAX,NP), the node neigbors. ! ! Input, integer NODE_NABOR_MAX, the maximum number of node neighbors. ! ! Input, integer NODE_NABOR_NUM(NP), the number of node neighbors. ! implicit none integer eqn_nabor_max integer node_nabor_max integer neqn integer np integer elem integer eqn_nabor(eqn_nabor_max,neqn) integer eqn_nabor_num(neqn) integer i integer idof integer ieqn integer indx(3,np) integer ip integer jp integer k integer nab(3) integer nb integer node_nabor(node_nabor_max,np) integer node_nabor_num(np) ! ! Initialize. ! eqn_nabor(1:eqn_nabor_max,1:neqn) = 0 do ieqn = 1, neqn eqn_nabor(1,ieqn) = ieqn end do eqn_nabor_num(1:neqn) = 1 ! ! For each node... ! do ip = 1, np ! ! For each node neighbor... ! do k = 1, node_nabor_num(ip) jp = node_nabor(k,ip) nb = 1 nab(nb) = indx(1,jp) nb = 2 nab(nb) = indx(2,jp) if ( indx(3,jp) /= 0 ) then nb = 3 nab(nb) = indx(3,jp) end if ! ! For each possible degree of freedom... ! do idof = 1, 3 ieqn = indx(idof,ip) if ( ieqn /= 0 ) then call i4vec_merge_a ( eqn_nabor_num(ieqn), eqn_nabor(1,ieqn), & nb, nab, & eqn_nabor_num(ieqn), eqn_nabor(1,ieqn) ) if ( eqn_nabor_max < eqn_nabor_num(ieqn) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'EQN_NABOR_SET - Fatal error!' write ( *, '(a)' ) ' Number of EQN neighbors exceeded.' write ( *, '(a,i8)' ) ' EQN_NABOR_MAX = ', eqn_nabor_max write ( *, '(a,i8)' ) & ' EQN_NABOR_NUM(IEQN) = ', eqn_nabor_num(ieqn) write ( *, '(a,i8)' ) ' IEQN = ', ieqn stop end if end if end do end do end do return end subroutine equation_print ( eqn, indx, neqn, np ) !*****************************************************************************80 ! !! EQUATION_PRINT prints the equation types. ! ! Modified: ! ! 05 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = 3 ) EQN(NEQN), the type of each equation. ! ! Input, integer INDX(3,NP), mapping from nodes to degrees of freedom. ! ! Input, integer NEQN, the number of finite element equations. ! ! Input, integer NP, the number of nodes. ! implicit none integer neqn integer np character ( len = 3 ) eqn(neqn) integer indx(3,np) integer ip write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'EQUATION_PRINT:' write ( *, '(a,i8)' ) ' Number of unknowns/equations = ', neqn write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Node U V P' write ( *, '(a)' ) ' ' do ip = 1, np if ( indx(3,ip) /= 0 ) then write ( *, '(i8,2x,i8,2x,a3,2x,i8,2x,a3,i8,2x,a3)' ) ip, & indx(1,ip), eqn(indx(1,ip)), & indx(2,ip), eqn(indx(2,ip)), & indx(3,ip), eqn(indx(3,ip)) else write ( *, '(i8,2x,i8,2x,a3,2x,i8,2x,a3,i8,2x,a3)' ) ip, & indx(1,ip), eqn(indx(1,ip)), & indx(2,ip), eqn(indx(2,ip)) end if end do return end subroutine fp ( a, bc_tag, detmap, eqn, g, indx, maxelm, maxeqn, maxnp, & maxquad1, maxquad2, maxside, navier, ncol, nelem, neqn, nlband, node, & np, nquad1, nquad2, nrow, nside, penalty1, penalty2, phi, region, nu_inv, & side_basis, side_elem, side_eqn, side_etam, side_etap, side_indx, & side_xsim, side_xsip, solver, squad1, wquad1, wquad2, xc, yc ) !*********************************************************************** ! !! FP computes the jacobian matrix A of the Navier Stokes residual. ! ! Discussion: ! ! FP differentiates the equations computed in FX. ! ! The differentiated Navier Stokes equations have the form: ! ! d UM-Eqn(I) / d U-Coef(J) = ! ! Integral dW(J)/dx * dW(I)/dx + dW(J)/dy * dW(I)/dy ! + nu_inv * ( W(J) * dUold/dx + Uold * dW(J) /dx + Vold * dW(J)/dy ) ! * W(I) dx dy ! ! d UM-Eqn(I) / d V-Coef(J) = ! ! Integral nu_inv * W(J) * dUold/dy * W(I) dx dy ! ! d UM-Eqn(I) / d P-Coef(J) = ! ! Integral nu_inv * dQ(J)/dx * W(I) dx dy ! ! d VM-Eqn(I) / d U-Coef(J) = ! ! Integral nu_inv * W(J) * dVold/dx * W(I) dx dy ! ! d VM-Eqn(I) / d V-Coef(J) = ! ! Integral dW(J)/dx * dW(I)/dx + dW(J)/dy * dW(I)/dy ! + nu_inv * ( Uold * dW(J)/dx + W(J) * dVold/dy + Vold * dW(J)/dy ) ! * W(I) dx dy ! ! d VM-Eqn(I) / d P-Coef(J) = ! ! Integral nu_inv * dQ(J)/dy * W(I) dx dy ! ! d PC-Eqn(I) / d U-Coef(J) = ! ! Integral dW(J)/dx * Q(I) dx dy ! ! d PC-Eqn(I) / d V-Coef(J) = ! ! Integral dW(J)/dy * Q(I) dx dy ! ! The boundary conditions must also be differentiated. ! ! Modified: ! ! 11 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real A(NROW,MAXEQN), the banded jacobian. ! ! real DETMAP(MAXQUAD2,MAXELM), the reference map determinant. ! ! Input, character ( len = 3 ) EQN(MAXEQN), the type of each equation. ! ! Input, real G(MAXEQN), the current solution vector. ! ! Input, integer INDX(3,MAXNP), mapping from nodes to degrees of freedom. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNP, the maximum number of nodes. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NEQN, the number of finite element equations. ! ! Input, integer NLBAND, the lower bandwidth of the matrix A. ! ! Input, integer NODE(6,MAXELM), the nodes that make up each element. ! ! Input, integer NP, the number of nodes. ! ! Input, integer NROW, the number of rows need to store the matrix A. ! ! Input, real PHI(MAXQUAD2,6,6,MAXELM), basis function values. ! ! Input, character ( len = 20 ) REGION, the flow problem. ! ! Input, real NU_INV, the value of the inverse viscosity. ! implicit none integer maxelm integer maxeqn integer maxnp integer maxquad1 integer maxquad2 integer maxside integer ncol integer nrow real ( kind = 8 ) a(nrow,ncol) integer bc_tag(maxnp) integer col real ( kind = 8 ) detmap(maxquad2,maxelm) real ( kind = 8 ) dpdx real ( kind = 8 ) dpdy real ( kind = 8 ) dqjdr real ( kind = 8 ) dqjdx real ( kind = 8 ) dqjdy real ( kind = 8 ) dtdr real ( kind = 8 ) dtds real ( kind = 8 ) dudr real ( kind = 8 ) duds real ( kind = 8 ) dudx real ( kind = 8 ) dudy real ( kind = 8 ) dvdr real ( kind = 8 ) dvds real ( kind = 8 ) dvdx real ( kind = 8 ) dvdy real ( kind = 8 ) dwideta real ( kind = 8 ) dwidxsi real ( kind = 8 ) dwidr real ( kind = 8 ) dwids real ( kind = 8 ) dwidx real ( kind = 8 ) dwidy real ( kind = 8 ) dwjdeta real ( kind = 8 ) dwjdxsi real ( kind = 8 ) dwjdr real ( kind = 8 ) dwjds real ( kind = 8 ) dwjdx real ( kind = 8 ) dwjdy character ( len = 3 ) eqn(maxeqn) real ( kind = 8 ) eta real ( kind = 8 ) etam real ( kind = 8 ) etap real ( kind = 8 ) g(maxeqn) integer i integer idof integer ielem integer ieqn integer ihor integer indx(3,maxnp) integer ip integer iprs integer iq integer iquad1 integer iquad2 integer iside integer iver integer j integer jhor integer jp integer jprs integer jq integer jver real ( kind = 8 ) length real ( kind = 8 ) n real ( kind = 8 ) navier integer nelem integer neqn integer nlband integer node(6,maxelm) real ( kind = 8 ) normal(2) integer np integer nquad1 integer nquad2 integer nside real ( kind = 8 ) p real ( kind = 8 ) penalty1 real ( kind = 8 ) penalty2 real ( kind = 8 ) phi(maxquad2,6,6,maxelm) real ( kind = 8 ) qi character ( len = 20 ) region real ( kind = 8 ) nu_inv integer row real ( kind = 8 ) s integer side_basis(3,maxside) integer side_elem(maxside) character ( len = 3 ) side_eqn(maxside) real ( kind = 8 ) side_etam(maxside) real ( kind = 8 ) side_etap(maxside) integer side_indx(3,maxside) real ( kind = 8 ) side_xsim(maxside) real ( kind = 8 ) side_xsip(maxside) character ( len = 20 ) solver real ( kind = 8 ) squad1(maxquad1) real ( kind = 8 ) t real ( kind = 8 ) tangent(2) real ( kind = 8 ) u real ( kind = 8 ) v real ( kind = 8 ) weight real ( kind = 8 ) wi real ( kind = 8 ) wj real ( kind = 8 ) wquad1(maxquad1) real ( kind = 8 ) wquad2(maxquad2) real ( kind = 8 ) x real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) xc(maxnp) real ( kind = 8 ) xsi real ( kind = 8 ) xsim real ( kind = 8 ) xsip real ( kind = 8 ) y real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 real ( kind = 8 ) yc(maxnp) ! ! Zero out the matrix. ! a(1:nrow,1:ncol) = 0.0D+00 ! !*****************************************************************************80 ! Handle equations associated with area integrals. !*****************************************************************************80 ! do ielem = 1, nelem ! ! Evaluate the integrand at the quadrature points. ! do iquad2 = 1, nquad2 weight = 0.5D+00 * wquad2(iquad2) * detmap(iquad2,ielem) ! ! For the given quadrature point, evaluate P, U and V. ! call uvp_quad_value ( dpdx, dpdy, dudx, dudy, dvdx, dvdy, & g, ielem, indx, iquad2, maxelm, maxeqn, maxnp, maxquad2, & node, p, phi, u, v ) ! ! Look at all the basis functions in the element IELEM. ! do iq = 1, 6 ip = node(iq,ielem) wi = phi(iquad2,iq,1,ielem) dwidx = phi(iquad2,iq,2,ielem) dwidy = phi(iquad2,iq,3,ielem) qi = phi(iquad2,iq,4,ielem) x = xc(ip) y = yc(ip) call boundary_shape ( bc_tag(ip), normal, region, tangent ) t = tangent(1) * u + tangent(2) * v dudr = tangent(1) * dudx + tangent(2) * dudy dvdr = tangent(1) * dvdx + tangent(2) * dvdy dwidr = tangent(1) * dwidx + tangent(2) * dwidy n = normal(1) * u + normal(2) * v duds = normal(1) * dudx + normal(2) * dudy dvds = normal(1) * dvdx + normal(2) * dvdy dwids = normal(1) * dwidx + normal(2) * dwidy dtdr = tangent(1) * dudr + tangent(2) * dvdr dtds = tangent(1) * duds + tangent(2) * dvds ihor = indx(1,ip) iver = indx(2,ip) iprs = indx(3,ip) ! ! Now compute the derivatives of the functions associated ! with U, V and P, with respect to the coefficients associated ! with basis vectors at each node of the element. ! do jq = 1, 6 jp = node(jq,ielem) wj = phi(iquad2,jq,1,ielem) dwjdx = phi(iquad2,jq,2,ielem) dwjdy = phi(iquad2,jq,3,ielem) dwjdr = tangent(1) * dwjdx + tangent(2) * dwjdy dwjds = normal(1) * dwjdx + normal(2) * dwjdy dqjdx = phi(iquad2,jq,5,ielem) dqjdy = phi(iquad2,jq,6,ielem) dqjdr = tangent(1) * dqjdx + tangent(2) * dqjdy jhor = indx(1,jp) jver = indx(2,jp) jprs = indx(3,jp) ! ! Contribution of the JHOR variable to the U equation. ! if ( eqn(ihor) == 'UM' ) then row = ihor - jhor + 2*nlband + 1 col = jhor a(row,col) = a(row,col) + weight * & ( dwjdx * dwidx + dwjdy * dwidy + nu_inv * navier * & ( wj * dudx + u * dwjdx + v * dwjdy ) * wi ) end if ! ! Contribution of the JHOR variable to the V equation. ! if ( eqn(iver) == 'VM' ) then row = iver - jhor + 2*nlband + 1 col = jhor a(row,col) = a(row,col) + weight * nu_inv * navier * dvdx * wj * wi else if ( eqn(iver) == 'TMP' ) then row = iver - jhor + 2*nlband + 1 col = jhor a(row,col) = a(row,col) + weight * penalty2 * & ( tangent(1) * dwjdr * dwidr + tangent(1) * dwjds * dwids & + nu_inv * navier * ( normal(1) * wj * dtds & + n * tangent(1) * dwjds + tangent(1) * wj * dtdr & + t * tangent(1) * dwjdr ) * wi ) end if ! ! Contribution of the JHOR variable to the P equation. ! if ( 0 < iprs ) then if ( eqn(iprs) == 'PC' ) then row = iprs - jhor + 2*nlband + 1 col = jhor a(row,col) = a(row,col) + weight * dwjdx * qi end if end if ! ! Contribution of the JVER variable to the U equation. ! if ( eqn(ihor) == 'UM' ) then row = ihor - jver + 2*nlband + 1 col = jver a(row,col) = a(row,col) + weight * nu_inv * navier * wj * dudy * wi end if ! ! Contribution of the JVER variable to the V equation. ! if ( eqn(iver) == 'VM' ) then row = iver - jver + 2*nlband + 1 col = jver a(row,col) = a(row,col) + weight * ( dwjdx * dwidx & + dwjdy * dwidy + nu_inv * navier * & ( u * dwjdx + wj * dvdy + v * dwjdy ) * wi ) else if ( eqn(iver) == 'TMP' ) then row = iver - jver + 2*nlband + 1 col = jver a(row,col) = a(row,col) + weight * penalty2 * & ( tangent(2) * dwjdr * dwidr + tangent(2) * dwjds * dwids & + nu_inv * navier * ( normal(2) * wj * dtds & + n * tangent(2) * dwjds + tangent(2) * wj * dtdr & + t * tangent(2) * dwjdr ) * wi ) end if ! ! Contribution of the JVER variable to the P equation. ! if ( 0 < iprs ) then if ( eqn(iprs) == 'PC' ) then row = iprs - jver + 2*nlband + 1 col = jver a(row,col) = a(row,col) + weight * dwjdy * qi end if end if ! ! Contribution of the JPRS variable to the U equation. ! if ( 0 < jprs ) then if ( eqn(ihor) == 'UM' ) then row = ihor - jprs + 2*nlband + 1 col = jprs a(row,col) = a(row,col) + weight * nu_inv * dqjdx * wi end if ! ! Contribution of the JPRS variable to the V equation. ! if ( eqn(iver) == 'VM' ) then row = iver - jprs + 2*nlband + 1 col = jprs a(row,col) = a(row,col) + weight * nu_inv * dqjdy * wi else if ( eqn(iver) == 'TMP' ) then row = iver - jprs + 2*nlband + 1 col = jprs a(row,col) = a(row,col) + weight * penalty2 * nu_inv * dqjdr * wi end if end if end do end do end do end do ! !*****************************************************************************80 ! Handle equations associated with line integrals along the boundary. !*****************************************************************************80 ! ! For each element side that is along the boundary. ! do iside = 1, nside ! ! If the element side condition is 'FS'... ! if ( side_eqn(iside) == 'FS' ) then ielem = side_elem(iside) etam = side_etam(iside) etap = side_etap(iside) xsim = side_xsim(iside) xsip = side_xsip(iside) ! ! Rough estimate of the length of the side. ! iq = side_basis(1,iside) ip = node(iq,ielem) x1 = xc(ip) y1 = yc(ip) iq = side_basis(2,iside) ip = node(iq,ielem) x2 = xc(ip) y2 = yc(ip) iq = side_basis(3,iside) ip = node(iq,ielem) x3 = xc(ip) y3 = yc(ip) length = sqrt ( ( x2 - x1 )**2 + ( y2 - y1 )**2 ) & + sqrt ( ( x3 - x2 )**2 + ( y3 - y2 )**2 ) ! ! KLUDGE. But it's right, I think. ! if ( x3 < x1 ) then length = - length end if do iquad1 = 1, nquad1 s = squad1(iquad1) weight = 0.5D+00 * wquad1(iquad1) * length eta = 0.5D+00 * ( ( 1.0D+00 - s ) * etam + ( 1.0D+00 + s ) * etap ) xsi = 0.5D+00 * ( ( 1.0D+00 - s ) * xsim + ( 1.0D+00 + s ) * xsip ) call x_of_xsi ( eta, ielem, maxelm, maxnp, node, x, xc, xsi, y, yc ) ! ! There are three basis functions on this portion of the boundary. ! do i = 1, 3 iq = side_basis(i,iside) ip = node(iq,ielem) idof = side_indx(i,iside) ieqn = indx(idof,ip) call ref_bf_q6 ( wi, dwideta, dwidxsi, eta, iq, xsi ) call boundary_shape ( bc_tag(ip), normal, region, tangent ) ! ! There are three basis functions that contribute to the values of ! U and V on this portion of the boundary. ! do j = 1, 3 jq = side_basis(j,iside) jp = node(jq,ielem) jhor = indx(1,jp) jver = indx(2,jp) call ref_bf_q6 ( wj, dwjdeta, dwjdxsi, eta, jq, xsi ) if ( eqn(ieqn) == 'TMP' ) then row = ieqn - jhor + 2*nlband+1 col = jhor a(row,col) = a(row,col) + penalty1 * weight & * tangent(1) * ( tangent(1) + tangent(2) ) * wj * wi row = ieqn - jver + 2*nlband+1 col = jver a(row,col) = a(row,col) + penalty1 * weight & * tangent(2) * ( tangent(1) + tangent(2) ) * wj * wi end if end do end do end do end if end do ! !*****************************************************************************80 ! Handle conditions associated with nodal values. !*****************************************************************************80 ! do ip = 1, np ihor = indx(1,ip) iver = indx(2,ip) iprs = indx(3,ip) x = xc(ip) y = yc(ip) call boundary_shape ( bc_tag(ip), normal, region, tangent ) if ( eqn(ihor) == 'UI' .or. eqn(ihor) == 'U0' ) then row = ihor-ihor+2*nlband+1 col = ihor a(row,col) = 1.0D+00 else if ( eqn(ihor) == 'N0' ) then row = ihor-ihor+2*nlband+1 col = ihor a(row,col) = normal(1) row = ihor-iver+2*nlband+1 col = ihor a(row,col) = normal(2) end if if ( eqn(iver) == 'VI' .or. eqn(iver) == 'V0' ) then row = iver-iver+2*nlband+1 col = iver a(row,col) = 1.0D+00 else if ( eqn(iver) == 'T0' ) then row = iver-ihor+2*nlband+1 col = ihor a(row,col) = tangent(1) row = iver-iver+2*nlband+1 col = iver a(row,col) = tangent(2) end if if ( 0 < iprs ) then if ( eqn(iprs) == 'PI' .or. eqn(iprs) == 'P0' ) then row = iprs-iprs+2*nlband+1 col = iprs a(row,col) = 1.0D+00 end if end if end do return end subroutine fp_band_info_print ( a, nrow, ncol, nlband ) !*****************************************************************************80 ! !! FP_BAND_INFO_PRINT prints information about the banded jacobian. ! ! Modified: ! ! 17 May 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(NROW,NCOL), the banded jacobian. ! ! Input, integer NROW, the number of rows need to store the matrix A. ! ! Input, integer NCOL, the number of columns, and the actual number of ! variables and equations. ! ! Input, integer NLBAND, the lower bandwidth of the matrix A. ! implicit none integer ncol integer nrow real ( kind = 8 ) a(nrow,ncol) integer i integer ihi integer ilo integer irow integer j integer jhi integer jlo integer nlband real ( kind = 8 ) norm write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FP_BAND_INFO_PRINT:' ! ! Row sums. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Row sums of absolute values:' write ( *, '(a)' ) ' ' do i = 1, ncol norm = 0.0D+00 jlo = max ( 1, i - nlband ) jhi = min ( ncol, i + nlband ) do j = jlo, jhi norm = norm + abs ( a(i-j+2*nlband+1,j) ) end do write ( *, '(i8,g14.6)' ) i, norm end do ! ! Column sums. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Column sums of absolute values:' write ( *, '(a)' ) ' ' do j = 1, ncol norm = 0.0D+00 ilo = max ( 1, j - nlband ) ihi = min ( ncol, j + nlband ) do i = ilo, ihi norm = norm + abs ( a(i-j+2*nlband+1,j) ) end do write ( *, '(i8,g14.6)' ) j, norm end do ! ! Print a particular row. ! irow = 160 if ( 1 <= irow .and. irow <= ncol ) then write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) 'Entries of row ', irow write ( *, '(a)' ) ' ' i = irow jlo = max ( 1, i - nlband ) jhi = min ( ncol, i + nlband ) do j = jlo, jhi if ( a(i-j+2*nlband+1,j) /= 0.0D+00 ) then write ( *, '(2i8,g14.6)' ) irow, j, a(i-j+2*nlband+1,j) end if end do end if return end subroutine fp_band_width ( indx, maxelm, maxnp, nelem, nlband, node, nrow ) !*****************************************************************************80 ! !! FP_BAND_WIDTH computes the lower band width of the Jacobian matrix. ! ! Discussion: ! ! The routine also finds NROW, the total number of rows required ! to store the matrix in LINPACK general band storage format. ! ! Modified: ! ! 07 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer INDX(3,MAXNP), mapping from nodes to degrees of freedom. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNP, the maximum number of nodes. ! ! Input, integer NELEM, the number of elements. ! ! Output, integer NLBAND, the lower bandwidth of the matrix A. ! ! Input, integer NODE(6,MAXELM), the nodes that make up each element. ! ! Input, integer NP, the number of nodes. ! ! Output, integer NROW, the number of rows needed to store the matrix A. ! implicit none integer maxelm integer maxnp integer i integer ielem integer indx(3,maxnp) integer ip integer ipp integer iq integer iqq integer iuk integer iukk integer j integer nelem integer nlband integer node(6,maxelm) integer nrow nlband = 0 do ielem = 1, nelem do iq = 1, 6 ip = node(iq,ielem) do iuk = 1, 3 i = indx(iuk,ip) if ( 0 < i ) then do iqq = 1, 6 ipp = node(iqq,ielem) do iukk = 1, 3 j = indx(iukk,ipp) if ( 0 < j ) then if ( nlband < j - i ) then nlband = j - i end if end if end do end do end if end do end do end do nrow = 3 * nlband + 1 return end subroutine fp_dif ( a, bc_tag, detmap, eqn, g, indx, maxelm, maxeqn, maxnp, & maxquad1, maxquad2, maxside, navier, ncol, nelem, neqn, nlband, node, np, & nquad1, nquad2, nrow, nside, penalty1, penalty2, phi, region, region_ymax, & res, res2, nu_inv, side_basis, side_elem, side_eqn, side_etam, side_etap, & side_indx, side_xsim, side_xsip, solver, squad1, wquad1, wquad2, xc, yc ) !*****************************************************************************80 ! !! FP_DIF estimates the jacobian matrix by finite differences. ! ! Modified: ! ! 19 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real A(NROW,MAXEQN), the banded jacobian. ! ! Input, real DETMAP(MAXQUAD2,MAXELM), the reference map determinant. ! ! Input, character ( len = 3 ) EQN(MAXEQN), the type of each equation. ! ! Input, real G(MAXEQN), the current solution vector. ! ! Input, integer INDX(3,MAXNP), mapping from nodes to degrees of freedom. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNP, the maximum number of nodes. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NEQN, the number of finite element equations. ! ! Input, integer NLBAND, the lower bandwidth of the matrix A. ! ! Input, integer NODE(6,MAXELM), the nodes that make up each element. ! ! Input, integer NP, the number of nodes. ! ! Input, integer NROW, the number of rows need to store the matrix A. ! ! Input, real PHI(MAXQUAD2,6,6,MAXELM), basis function values. ! ! Input, character ( len = 20 ) REGION, the flow problem. ! ! Input, real NU_INV, the value of the inverse viscosity. ! implicit none integer maxelm integer maxeqn integer maxnp integer maxquad1 integer maxquad2 integer maxside integer ncol integer nrow real ( kind = 8 ) a(nrow,ncol) integer bc_tag(maxnp) real ( kind = 8 ) detmap(maxquad2,maxelm) real ( kind = 8 ) dg real ( kind = 8 ), parameter :: EPS = 0.001D+00 character ( len = 3 ) eqn(maxeqn) real ( kind = 8 ) g(maxeqn) real ( kind = 8 ) gsave integer i integer ieqn integer ihi integer ilo integer indx(3,maxnp) integer irow integer j real ( kind = 8 ) navier integer nelem integer neqn integer nlband integer node(6,maxelm) integer np integer nquad1 integer nquad2 integer nside real ( kind = 8 ) penalty1 real ( kind = 8 ) penalty2 real ( kind = 8 ) phi(maxquad2,6,6,maxelm) character ( len = 20 ) region real ( kind = 8 ) region_ymax real ( kind = 8 ) res(maxeqn) real ( kind = 8 ) res2(maxeqn) real ( kind = 8 ) nu_inv integer side_basis(3,maxside) integer side_elem(maxside) character ( len = 3 ) side_eqn(maxside) real ( kind = 8 ) side_etam(maxside) real ( kind = 8 ) side_etap(maxside) integer side_indx(3,maxside) real ( kind = 8 ) side_xsim(maxside) real ( kind = 8 ) side_xsip(maxside) character ( len = 20 ) solver real ( kind = 8 ) squad1(maxquad1) real ( kind = 8 ) wquad1(maxquad1) real ( kind = 8 ) wquad2(maxquad2) real ( kind = 8 ) xc(maxnp) real ( kind = 8 ) yc(maxnp) ! ! Zero out the matrix. ! a(1:nrow,1:ncol) = 0.0D+00 ! ! Evaluate F(X). ! call fx ( bc_tag, detmap, eqn, g, indx, maxelm, maxeqn, maxnp, maxquad1, & maxquad2, maxside, navier, nelem, neqn, node, np, nquad1, nquad2, nside, & penalty1, penalty2, phi, region, region_ymax, res, nu_inv, side_basis, & side_elem, side_eqn, side_etam, side_etap, side_indx, side_xsim, & side_xsip, squad1, wquad1, wquad2, xc, yc ) ! ! Compare with F(X+deltaX(J)) ! do j = 1, neqn gsave = g(j) g(j) = ( 1.0D+00 + EPS ) * g(j) + sign ( EPS, g(j) ) dg = g(j) - gsave call fx ( bc_tag, detmap, eqn, g, indx, maxelm, maxeqn, maxnp, maxquad1, & maxquad2, maxside, navier, nelem, neqn, node, np, nquad1, nquad2, nside, & penalty1, penalty2, phi, region, region_ymax, res2, nu_inv, side_basis, & side_elem, side_eqn, side_etam, side_etap, side_indx, side_xsim, & side_xsip, squad1, wquad1, wquad2, xc, yc ) ilo = max ( 1, j - nlband ) ihi = min ( neqn, j + nlband ) do ieqn = ilo, ihi irow = ieqn - j + 2 * nlband + 1 a(irow,j) = ( res2(ieqn) - res(ieqn) ) / dg end do g(j) = gsave end do return end subroutine fx ( bc_tag, detmap, eqn, g, indx, maxelm, maxeqn, maxnp, & maxquad1, maxquad2, maxside, navier, nelem, neqn, node, np, nquad1, & nquad2, nside, penalty1, penalty2, phi, region, region_ymax, res, nu_inv, & side_basis, side_elem, side_eqn, side_etam, side_etap, side_indx, & side_xsim, side_xsip, squad1, wquad1, wquad2, xc, yc ) !*****************************************************************************80 ! !! FX computes the residual of the Navier Stokes equations. ! ! Discussion: ! ! By setting the variable NAVIER to 0, you get the residual for ! the Stokes equations! ! ! ! The Navier Stokes equations have the form: ! ! UM-Eqn(I) = ! ! Integral ( ! dU/dx * dW(I)/dx + dU/dy * dW(I)/dy + ! nu_inv * ( U * dU/dx + V * dU/dy + dP/dx ) * W(I) ! ) dx dy = 0 ! ! VM-Eqn(I) = ! ! Integral ( ! dV/dx * dW(I)/dx + dV/dy * dW(I)/dy + ! nu_inv * ( U * dV/dx + V * dV/dy + dP/dy ) * W(I) ! ) dx dy = 0 ! ! PC-Eqn(I) = ! ! Integral ( ! ( dU/dx + dV/dy ) * Q(I) ! ) dx dy = 0 ! ! W is a basis function for U and V, and Q is a basis function for P. ! ! A penalty term may be included with the momentum equations: ! ! UMP-Eqn(I) = ! ! Integral ( ! dU/dx * dW(I)/dx + dU/dy * dW(I)/dy + ! ( ! nu_inv * ( U * dU/dx + V * dU/dy + dP/dx ) + ! ??? / penalty ! ) * W(I) ! ) dx dy = 0 ! ! VMP-Eqn(I) = ! ! Integral ( ! dV/dx * dW(I)/dx + dV/dy * dW(I)/dy + ! ( ! nu_inv * ( U * dV/dx + V * dV/dy + dP/dy ) + ! ??? / penalty ! ) * W(I) ! ) dx dy = 0 ! ! where PENALTY is a user specified nonzero scalar. ! ! Modified: ! ! 11 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real DETMAP(MAXQUAD2,MAXELM), the reference map determinant. ! ! Input, character ( len = 3 ) EQN(MAXEQN), the type of each equation. ! ! Input, real G(MAXEQN), the current solution. ! ! Input, integer INDX(3,MAXNP), mapping from nodes to degrees of freedom. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNP, the maximum number of nodes. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NEQN, the number of finite element equations. ! ! Input, integer NODE(6,MAXELM), the nodes that make up each element. ! ! Input, integer NP, the number of nodes. ! ! Input, real PHI(MAXQUAD2,6,6,MAXELM), basis function values. ! ! Input, character ( len = 20 ) REGION, the flow problem. ! ! Output, real RES(MAXEQN), the residual. ! ! Input, real NU_INV, the inverse viscosity. ! ! Input, real XC(MAXNP), the X coordinates of the nodes. ! ! Input, real YC(MAXNP), the Y coordinates of the nodes. ! implicit none integer maxelm integer maxeqn integer maxnp integer maxquad1 integer maxquad2 integer maxside integer bc_tag(maxnp) real ( kind = 8 ) detmap(maxquad2,maxelm) real ( kind = 8 ) dpdr real ( kind = 8 ) dpdx real ( kind = 8 ) dpdy real ( kind = 8 ) dtdr real ( kind = 8 ) dtds real ( kind = 8 ) dudr real ( kind = 8 ) duds real ( kind = 8 ) dudx real ( kind = 8 ) dudy real ( kind = 8 ) dvdr real ( kind = 8 ) dvds real ( kind = 8 ) dvdx real ( kind = 8 ) dvdy real ( kind = 8 ) dwdeta real ( kind = 8 ) dwdxsi real ( kind = 8 ) dwidr real ( kind = 8 ) dwids real ( kind = 8 ) dwidx real ( kind = 8 ) dwidy character ( len = 3 ) eqn(maxeqn) real ( kind = 8 ) eta real ( kind = 8 ) etam real ( kind = 8 ) etap real ( kind = 8 ) g(maxeqn) integer i integer idof integer ielem integer ieqn integer ihor integer indx(3,maxnp) integer ip integer iprs integer iq integer iquad1 integer iquad2 integer iside integer iver real ( kind = 8 ) length real ( kind = 8 ) n real ( kind = 8 ) navier integer nelem integer neqn integer node(6,maxelm) real ( kind = 8 ) normal(2) integer np integer nquad1 integer nquad2 integer nside real ( kind = 8 ) p real ( kind = 8 ) pbc real ( kind = 8 ) penalty1 real ( kind = 8 ) penalty2 real ( kind = 8 ) phi(maxquad2,6,6,maxelm) real ( kind = 8 ) qi character ( len = 20 ) region real ( kind = 8 ) region_ymax real ( kind = 8 ) res(maxeqn) real ( kind = 8 ) nu_inv real ( kind = 8 ) s integer side_basis(3,maxside) integer side_elem(maxside) character ( len = 3 ) side_eqn(maxside) real ( kind = 8 ) side_etam(maxside) real ( kind = 8 ) side_etap(maxside) integer side_indx(3,maxside) real ( kind = 8 ) side_xsim(maxside) real ( kind = 8 ) side_xsip(maxside) real ( kind = 8 ) squad1(maxquad1) real ( kind = 8 ) t real ( kind = 8 ) tangent(2) real ( kind = 8 ) term real ( kind = 8 ) u real ( kind = 8 ) ubc real ( kind = 8 ) v real ( kind = 8 ) vbc real ( kind = 8 ) weight real ( kind = 8 ) wi real ( kind = 8 ) wquad1(maxquad1) real ( kind = 8 ) wquad2(maxquad2) real ( kind = 8 ) x real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) xc(maxnp) real ( kind = 8 ) xsi real ( kind = 8 ) xsim real ( kind = 8 ) xsip real ( kind = 8 ) y real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 real ( kind = 8 ) yc(maxnp) ! ! Zero out the residual vector. ! res(1:neqn) = 0.0D+00 ! !*****************************************************************************80 ! Handle equations associated with area integrals. !*****************************************************************************80 ! ! Consider an element. ! do ielem = 1, nelem ! ! Evaluate the integrand at each of the quadrature points. ! do iquad2 = 1, nquad2 weight = 0.5D+00 * wquad2(iquad2) * detmap(iquad2,ielem) ! ! Evaluate P, U and V at a quadrature point. ! call uvp_quad_value ( dpdx, dpdy, dudx, dudy, dvdx, dvdy, g, ielem, & indx, iquad2, maxelm, maxeqn, maxnp, maxquad2, node, p, phi, u, v ) ! ! Look at all basis functions in the element IELEM. ! do iq = 1, 6 ip = node(iq,ielem) wi = phi(iquad2,iq,1,ielem) dwidx = phi(iquad2,iq,2,ielem) dwidy = phi(iquad2,iq,3,ielem) qi = phi(iquad2,iq,4,ielem) x = xc(ip) y = yc(ip) call boundary_shape ( bc_tag(ip), normal, region, tangent ) t = tangent(1) * u + tangent(2) * v dudr = tangent(1) * dudx + tangent(2) * dudy dvdr = tangent(1) * dvdx + tangent(2) * dvdy dpdr = tangent(1) * dpdx + tangent(2) * dpdy dwidr = tangent(1) * dwidx + tangent(2) * dwidy n = normal(1) * u + normal(2) * v duds = normal(1) * dudx + normal(2) * dudy dvds = normal(1) * dvdx + normal(2) * dvdy dwids = normal(1) * dwidx + normal(2) * dwidy dtdr = tangent(1) * dudr + tangent(2) * dvdr dtds = tangent(1) * duds + tangent(2) * dvds ihor = indx(1,ip) iver = indx(2,ip) iprs = indx(3,ip) ! ! The horizontal velocity momentum equation. ! if ( eqn(ihor) == 'UM' ) then res(ihor) = res(ihor) + weight * ( dudx * dwidx + dudy * dwidy & + nu_inv * ( navier * u * dudx + navier * v * dudy + dpdx ) * wi ) end if ! ! The vertical velocity momentum equation, or ! the tangential momentum equation, penalized. ! if ( eqn(iver) == 'VM' ) then res(iver) = res(iver) + weight * ( dvdx * dwidx + dvdy * dwidy & + nu_inv * ( navier * u * dvdx + navier * v * dvdy + dpdy ) * wi ) else if ( eqn(iver) == 'TMP' ) then res(iver) = res(iver) + weight * penalty2 * & ( dtdr * dwidr + dtds * dwids + nu_inv * & ( navier * n * dtds + navier * t * dtdr + dpdr ) * wi ) end if ! ! The pressure equations. ! if ( 0 < iprs ) then if ( eqn(iprs) == 'PC' ) then res(iprs) = res(iprs) + weight * ( dudx + dvdy ) * qi end if end if end do end do end do ! !*****************************************************************************80 ! Handle equations associated with line integrals along the boundary. !*****************************************************************************80 ! ! For each element side that is along the boundary. ! do iside = 1, nside ! ! If the element side condition is 'FS'... ! if ( side_eqn(iside) == 'FS' ) then ielem = side_elem(iside) etam = side_etam(iside) etap = side_etap(iside) xsim = side_xsim(iside) xsip = side_xsip(iside) ! ! Rough estimate of the length of the side. ! iq = side_basis(1,iside) ip = node(iq,ielem) x1 = xc(ip) y1 = yc(ip) iq = side_basis(2,iside) ip = node(iq,ielem) x2 = xc(ip) y2 = yc(ip) iq = side_basis(3,iside) ip = node(iq,ielem) x3 = xc(ip) y3 = yc(ip) length = sqrt ( ( x2 - x1 )**2 + ( y2 - y1 )**2 ) & + sqrt ( ( x3 - x2 )**2 + ( y3 - y2 )**2 ) ! ! KLUDGE. But it's right, I think. ! if ( x3 < x1 ) then length = - length end if do iquad1 = 1, nquad1 s = squad1(iquad1) weight = 0.5D+00 * wquad1(iquad1) * length eta = 0.5D+00 * ( ( 1.0D+00 - s ) * etam + ( 1.0D+00 + s ) * etap ) xsi = 0.5D+00 * ( ( 1.0D+00 - s ) * xsim + ( 1.0D+00 + s ) * xsip ) call x_of_xsi ( eta, ielem, maxelm, maxnp, node, x, xc, xsi, y, yc ) call uvp_value ( eta, g, ielem, indx, maxelm, maxeqn, maxnp, node, & p, u, v, xsi ) ! ! There are three basis functions on this portion of the boundary. ! do i = 1, 3 iq = side_basis(i,iside) ip = node(iq,ielem) idof = side_indx(i,iside) ieqn = indx(idof,ip) call ref_bf_q6 ( wi, dwdeta, dwdxsi, eta, iq, xsi ) call boundary_shape ( bc_tag(ip), normal, region, tangent ) ! ! This term is independent of the sign of the tangent... ! term = ( u * tangent(1) + v * tangent(2) ) & * ( tangent(1) + tangent(2) ) * wi if ( eqn(ieqn) == 'TMP' ) then res(ieqn) = res(ieqn) + penalty1 * weight * term end if end do end do end if end do ! !*****************************************************************************80 ! Handle conditions associated with nodal values. !*****************************************************************************80 ! do ip = 1, np ihor = indx(1,ip) iver = indx(2,ip) iprs = indx(3,ip) x = xc(ip) y = yc(ip) call boundary_value ( bc_tag(ip), pbc, region, region_ymax, ubc, vbc, y ) call boundary_shape ( bc_tag(ip), normal, region, tangent ) ! ! Specified value of horizontal velocity at a node. ! if ( eqn(ihor) == 'UI' ) then res(ihor) = g(ihor) - ubc else if ( eqn(ihor) == 'U0' ) then res(ihor) = g(ihor) ! ! Velocity dot local normal vector at a node is zero. ! else if ( eqn(ihor) == 'N0' ) then res(ihor) = g(ihor) * normal(1) + g(iver) * normal(2) end if ! ! Specified value of vertical velocity at a node. ! if ( eqn(iver) == 'VI' ) then res(iver) = g(iver) - vbc else if ( eqn(iver) == 'V0' ) then res(iver) = g(iver) ! ! Velocity dot local tangent vector at a node is zero. ! else if ( eqn(iver) == 'T0' ) then res(iver) = g(ihor) * tangent(1) + g(iver) * tangent(2) end if ! ! Specified value of pressure at a node. ! if ( 0 < iprs ) then if ( eqn(iprs) == 'PI' ) then res(iprs) = g(iprs) - pbc else if ( eqn(iprs) == 'P0' ) then res(iprs) = g(iprs) end if end if end do return end subroutine geometry ( bc_tag, eqn, indx, maxelm, maxeqn, maxnp, nelem, neqn, & node, np, nx, ny, region, region_xmax, region_ymax, xc, yc ) !*****************************************************************************80 ! !! GEOMETRY sets up the geometry for any problem. ! ! Modified: ! ! 11 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(6,MAXELM), the nodes that make up each element. ! ! Input, integer NP, the number of nodes. ! ! Input, integer NX, the number of elements along a horizontal line. ! ! Input, integer NY, the number of elements along a vertical line. ! implicit none integer maxelm integer maxeqn integer maxnp integer bc_tag(maxnp) character ( len = 3 ) eqn(maxeqn) integer indx(3,maxnp) integer nelem integer neqn integer node(6,maxelm) integer np integer nx integer ny character ( len = * ) region real ( kind = 8 ) region_xmax real ( kind = 8 ) region_ymax real ( kind = 8 ) xc(maxnp) real ( kind = 8 ) yc(maxnp) if ( region == 'CAVITY' ) then call geometry_cavity ( bc_tag, eqn, indx, maxelm, maxeqn, maxnp, nelem, & neqn, node, np, nx, ny, region_xmax, region_ymax, xc, yc ) else if ( region == 'CHANNEL' ) then call geometry_channel ( bc_tag, eqn, indx, maxelm, maxeqn, maxnp, nelem, & neqn, node, np, nx, ny, region_xmax, region_ymax, xc, yc ) else if ( region == 'FREESLIP' ) then call geometry_freeslip ( bc_tag, eqn, indx, maxelm, maxeqn, maxnp, nelem, & neqn, node, np, nx, ny, region_xmax, region_ymax, xc, yc ) else if ( region == 'STEP' ) then call geometry_step ( bc_tag, eqn, indx, maxelm, maxeqn, maxnp, nelem, & neqn, node, np, nx, ny, region_xmax, region_ymax, xc, yc ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GEOMETRY - Fatal error!' write ( *, '(a)' ) ' Illegal region.' stop end if return end subroutine geometry_cavity ( bc_tag, eqn, indx, maxelm, maxeqn, maxnp, & nelem, neqn, node, np, nx, ny, region_xmax, region_ymax, xc, yc ) !*****************************************************************************80 ! !! GEOMETRY_CAVITY sets up the geometry for the cavity problem. ! ! Discussion: ! ! The flow region is a square, of width and height 1, with no slip walls on ! the sides and bottom, and open on the top. A tangential force is ! applied uniformly along the top, enforced with a boundary condition ! of the form: ! ! U(X,1) = 1 ! V(X,1) = 0 ! ! Modified: ! ! 11 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(6,MAXELM), the nodes that make up each element. ! ! Input, integer NP, the number of nodes. ! ! Input, integer NX, the number of elements along a horizontal line. ! ! Input, integer NY, the number of elements along a vertical line. ! implicit none integer maxelm integer maxeqn integer maxnp integer bc_tag(maxnp) integer cc integer ee character ( len = 3 ) eqn(maxeqn) integer icol integer indx(3,maxnp) integer irow integer ncol integer ne integer nelem integer neqn integer nn integer node(6,maxelm) integer np integer nrow integer nw integer nx integer ny logical pset real ( kind = 8 ) region_xmax real ( kind = 8 ) region_ymax integer se integer ss integer sw integer ww real ( kind = 8 ) xc(maxnp) real ( kind = 8 ) yc(maxnp) pset = .FALSE. nelem = 0 neqn = 0 np = 0 ncol = 2 * nx + 1 nrow = 2 * ny + 1 do icol = 1, ncol do irow = 1, nrow ! ! 1. Generate the next node. ! np = np + 1 ! ! 2. Assign X, Y coordinates to the node. ! xc(np) = region_xmax * real ( icol - 1, kind = 8 ) & / real ( ncol - 1, kind = 8 ) yc(np) = region_ymax * real ( irow - 1, kind = 8 ) & / real ( nrow - 1, kind = 8 ) ! ! 3. Assign a BC_TAG value. ! if ( irow == 1 ) then if ( icol == 1 ) then bc_tag(np) = 1 else if ( icol < ncol ) then bc_tag(np) = 2 else bc_tag(np) = 3 end if else if ( irow < nrow ) then if ( icol == 1 ) then bc_tag(np) = 8 else if ( icol < ncol ) then bc_tag(np) = 0 else bc_tag(np) = 4 end if else if ( irow == nrow ) then if ( icol == 1 ) then bc_tag(np) = 7 else if ( icol < ncol ) then bc_tag(np) = 6 else bc_tag(np) = 5 end if end if ! ! 4. Assign equation indices to the node. ! ! ! Node along top. ! "Horizontal Velocity Specified" ! "Vertical Velocity 0". ! if ( irow == nrow ) then neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'UI' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'V0' ! ! Node along sides or bottom. ! "Horizontal Velocity 0" ! "Vertical Velocity 0". ! else if ( & icol == 1 .or. & icol == ncol .or. & irow == 1 ) then neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'U0' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'V0' ! ! Interior node. ! else neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'UM' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'VM' end if ! ! Set the pressure condition. ! if ( mod ( irow, 2 ) == 0 .or. mod ( icol, 2 ) == 0 ) then indx(3,np) = 0 else if ( & .not. pset .and. & 2 <= neqn .and. & eqn(neqn-1) == 'UM' .and. & eqn(neqn) == 'VM' ) then neqn = neqn + 1 indx(3,np) = neqn eqn(neqn) = 'P0' pset = .true. else neqn = neqn + 1 indx(3,np) = neqn eqn(neqn) = 'PC' end if ! ! 5. Generate the next pair of elements? ! ! If both the row and the column are odd, and we're not in the first ! column or row, then we can define two new triangular elements ! based at the node. ! ! NW---NN---NE ! |\ | ! | \ | ! | \ | ! WW CC EE ! | \ | ! | \ | ! | \| ! SW---SS---SE ! if ( & 1 < irow .and. mod ( irow, 2 ) == 1 .and. & 1 < icol .and. mod ( icol, 2 ) == 1 ) then sw = np - 2 * nrow - 2 ww = np - 2 * nrow - 1 nw = np - 2 * nrow ss = np - nrow - 2 cc = np - nrow - 1 nn = np - nrow se = np - 2 ee = np - 1 ne = np nelem = nelem + 1 node(1,nelem) = se node(2,nelem) = nw node(3,nelem) = sw node(4,nelem) = ss node(5,nelem) = cc node(6,nelem) = ww nelem = nelem + 1 node(1,nelem) = nw node(2,nelem) = se node(3,nelem) = ne node(4,nelem) = nn node(5,nelem) = cc node(6,nelem) = ee end if end do end do return end subroutine geometry_channel ( bc_tag, eqn, indx, maxelm, maxeqn, maxnp, & nelem, neqn, node, np, nx, ny, region_xmax, region_ymax, xc, yc ) !*****************************************************************************80 ! !! GEOMETRY_CHANNEL sets up the geometry for the channel problem. ! ! Discussion: ! ! Consider a horizontal channel of constant height H, and of length L. ! ! Suppose a parabolic inflow is specified at the left hand opening, ! where X=0, of the form ! ! U(0,Y) = S * Y * (H-Y) ! V(0,Y) = 0 ! P(0,Y) = 0 ! ! where S is any value. ! ! Then the following functions (U,V,P) solve the Navier Stokes ! equations in the region: ! ! U(X,Y) = S * Y * (H-Y) ! V(X,Y) = 0 ! P(X,Y) = -2 * S * X * RE ! ! The standard problem we use has H = 3, L = 10, and chooses ! S = (4/9)*Lambda, so that the maximum value of the parabolic inflow, ! at Y = H/2, is Lambda. Then our formula becomes: ! ! U(X,Y) = (4/9) * Lambda * Y * (3-Y) ! V(X,Y) = 0 ! P(X,Y) = -2 * (4/9) * Lambda * X * RE ! ! Modified: ! ! 22 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(6,MAXELM), the nodes that make up each element. ! ! Input, integer NP, the number of nodes. ! ! Input, integer NX, the number of elements along a horizontal line. ! ! Input, integer NY, the number of elements along a vertical line. ! implicit none integer maxelm integer maxeqn integer maxnp integer bc_tag(maxnp) integer cc integer ee character ( len = 3 ) eqn(maxeqn) integer icol integer indx(3,maxnp) integer irow integer ncol integer ne integer nelem integer neqn integer nn integer node(6,maxelm) integer np integer nrow integer nw integer nx integer ny real ( kind = 8 ) region_xmax real ( kind = 8 ) region_ymax integer se integer ss integer sw integer ww real ( kind = 8 ) xc(maxnp) real ( kind = 8 ) yc(maxnp) nelem = 0 neqn = 0 np = 0 ncol = 2 * nx + 1 nrow = 2 * ny + 1 do icol = 1, ncol do irow = 1, nrow ! ! 1. Generate the next node. ! np = np + 1 ! ! 2. Assign X, Y coordinates to the node. ! xc(np) = region_xmax * real ( icol - 1 ) & / real ( ncol - 1, kind = 8 ) yc(np) = region_ymax * real ( irow - 1 ) & / real ( nrow - 1, kind = 8 ) ! ! 3. Assign a BC_TAG value. ! if ( irow == 1 ) then if ( icol == 1 ) then bc_tag(np) = 1 else if ( icol < ncol ) then bc_tag(np) = 2 else bc_tag(np) = 3 end if else if ( irow < nrow ) then if ( icol == 1 ) then bc_tag(np) = 8 else if ( icol < ncol ) then bc_tag(np) = 0 else bc_tag(np) = 4 end if else if ( irow == nrow ) then if ( icol == 1 ) then bc_tag(np) = 7 else if ( icol < ncol ) then bc_tag(np) = 6 else bc_tag(np) = 5 end if end if ! ! 4. Assign equation indices to the node. ! ! Corners. ! if ( & ( icol == 1 .or. icol == ncol ) .and. & ( irow == 1 .or. irow == nrow ) ) then neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'U0' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'V0' ! ! Node within left hand inflow boundary. ! H specified. ! V 0. ! else if ( icol == 1 ) then neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'UI' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'V0' ! ! Node within right hand outflow boundary. ! V = 0. ! else if ( icol == ncol ) then neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'UM' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'V0' ! ! Node within fixed wall. ! "Normal Velocity Zero" ! "Tangential Velocity Zero". ! else if ( irow == 1 .or. irow == nrow ) then neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'N0' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'T0' ! ! Interior node. ! else neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'UM' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'VM' end if ! ! Set the pressure condition. ! if ( mod ( irow, 2 ) == 0 .or. mod ( icol, 2 ) == 0 ) then indx(3,np) = 0 else if ( irow == nrow .and. icol == ncol ) then neqn = neqn + 1 indx(3,np) = neqn eqn(neqn) = 'PI' else neqn = neqn + 1 indx(3,np) = neqn eqn(neqn) = 'PC' end if ! ! 5. Generate the next pair of elements? ! ! If both the row and the column are odd, and we're not in the first ! column or row, then we can define two new triangular elements ! based at the node. ! ! NW---NN---NE ! |\ | ! | \ | ! | \ | ! WW CC EE ! | \ | ! | \ | ! | \| ! SW---SS---SE ! if ( 1 < irow .and. mod ( irow, 2 ) == 1 .and. & 1 < icol .and. mod ( icol, 2 ) == 1 ) then sw = np - 2 * nrow - 2 ww = np - 2 * nrow - 1 nw = np - 2 * nrow ss = np - nrow - 2 cc = np - nrow - 1 nn = np - nrow se = np - 2 ee = np - 1 ne = np nelem = nelem + 1 node(1,nelem) = se node(2,nelem) = nw node(3,nelem) = sw node(4,nelem) = ss node(5,nelem) = cc node(6,nelem) = ww nelem = nelem + 1 node(1,nelem) = nw node(2,nelem) = se node(3,nelem) = ne node(4,nelem) = nn node(5,nelem) = cc node(6,nelem) = ee end if end do end do return end subroutine geometry_freeslip ( bc_tag, eqn, indx, maxelm, maxeqn, maxnp, & nelem, neqn, node, np, nx, ny, region_xmax, region_ymax, xc, yc ) !*****************************************************************************80 ! !! GEOMETRY_FREESLIP sets up the geometry for the freeslip problem. ! ! Modified: ! ! 11 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(6,MAXELM), the nodes that make up each element. ! ! Input, integer NP, the number of nodes. ! ! Input, integer NX, the number of elements along a horizontal line. ! ! Input, integer NY, the number of elements along a vertical line. ! implicit none integer maxelm integer maxeqn integer maxnp integer bc_tag(maxnp) integer cc integer ee character ( len = 3 ) eqn(maxeqn) integer icol integer indx(3,maxnp) integer irow integer ncol integer ne integer nelem integer neqn integer nn integer node(6,maxelm) integer np integer nrow integer nw integer nx integer ny real ( kind = 8 ) region_xmax real ( kind = 8 ) region_ymax integer se integer ss integer sw integer ww real ( kind = 8 ) xc(maxnp) real ( kind = 8 ) yc(maxnp) nelem = 0 neqn = 0 np = 0 ncol = 2 * nx + 1 nrow = 2 * ny + 1 do icol = 1, ncol do irow = 1, nrow ! ! 1. Generate the next node. ! np = np + 1 ! ! 2. Assign X, Y coordinates to the node. ! xc(np) = region_xmax * real ( icol - 1, kind = 8 ) & / real ( ncol - 1, kind = 8 ) yc(np) = region_ymax * real ( irow - 1, kind = 8 ) & / real ( nrow - 1, kind = 8 ) ! ! 3. Assign a BC_TAG value. ! if ( irow == 1 ) then if ( icol == 1 ) then bc_tag(np) = 1 else if ( icol < ncol ) then bc_tag(np) = 2 else bc_tag(np) = 3 end if else if ( irow < nrow ) then if ( icol == 1 ) then bc_tag(np) = 8 else if ( icol < ncol ) then bc_tag(np) = 0 else bc_tag(np) = 4 end if else if ( irow == nrow ) then if ( icol == 1 ) then bc_tag(np) = 7 else if ( icol < ncol ) then bc_tag(np) = 6 else bc_tag(np) = 5 end if end if ! ! 4. Assign equation indices to the node. ! ! Corners. ! if ( icol == 1 .and. & ( irow == 1 .or. irow == nrow ) ) then neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'N0' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'TMP' else if ( icol == ncol .and. & ( irow == 1 .or. irow == nrow ) ) then neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'N0' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'TMP' ! ! Node within left hand inflow boundary. ! H specified. ! V = 0. ! else if ( icol == 1 ) then neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'UI' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'V0' ! ! Node within right hand outflow boundary. ! V = 0. ! else if ( icol == ncol ) then neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'UM' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'V0' ! ! Node within fixed wall. ! "Normal Velocity Zero" ! "Tangential Slip / Momentum / Penalty term". ! else if ( irow == 1 .or. irow == nrow ) then neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'N0' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'TMP' ! ! Interior node. ! else neqn = neqn + 1 indx(1,np) = neqn eqn(neqn) = 'UM' neqn = neqn + 1 indx(2,np) = neqn eqn(neqn) = 'VM' end if ! ! Set the pressure condition. ! if ( mod ( irow, 2 ) == 0 .or. mod ( icol, 2 ) == 0 ) then indx(3,np) = 0 else if ( irow == nrow .and. icol == ncol ) then neqn = neqn + 1 indx(3,np) = neqn eqn(neqn) = 'P0' else neqn = neqn + 1 indx(3,np) = neqn eqn(neqn) = 'PC' end if ! ! 5