c----------------------------------------------------------------------- c implicit.f c Written by Mijan Huq c Physics 387n c ******************T-shaped stencil****************** c c The following subroutine takes in as its arguments u at n and n-1 and c computes u at n+1. The routine requires a temporary array of size NxN c where N is the size of the 1D grid on which we are solving the wave c equation. c To form the matrix we have to take into account the boundary conditions. c Here we take zero boundary conditions at u_1=0 and u_N=0 c nx := grid spacing c dt := time step c dx := space step c unp:= u at n+1 c unm:= u at n-1 c un := u at n c mat:= Matrix to be formed c c ipvt := Integer array of dimension N. Pivot array for LINPACK c rhs := array that will contain rhs of matrix equation and solution on return c from LINPACK. Could recycle unp and use as rhs and solution... subroutine impl_update1(nx, dt, dx, unp, unm, un, mat, ipvt, rhs) implicit none integer nx, ipvt(nx) real*8 dt, dx c Grid function declarations real*8 unp(nx), unm(nx), un(nx) c Matrix and LINPACK arrays real*8 mat(nx,nx), rhs(nx) c Local variables integer i, j, job real*8 m, m2, fact, rcond logical trace trace = .TRUE. m = dt / dx m2= m**2 fact = 1.0 + 2.0 * m2 c write(0,*)m,m2,fact c First zero out matrix do i = 1, nx do j = 1, nx mat(i,j) = 0.0 end do end do c Form the matrix taking into account the boundary conditions c Here we take fixed boundaries. The fixed boundaries yield that c unp(1) = 0.0 and unp(nx) = 0.0 c We are left with a nx-2 by nx-2 matrix to solve for the interior points c First initialize row 1 and row nx mat(1,1) = 1.0 mat(2,2) = fact mat(2,3) = -m2 mat(nx-1,nx-1) = fact mat(nx-1,nx-2) = -m2 mat(nx,nx) = 1.0 c Set up rest of matrix do i = 3, nx-2 mat(i,i-1) = -m2 mat(i,i) = fact mat(i,i+1) = -m2 end do c Form right hand side of matrix equation c Don't forget that we require the interior only c Zero boundaries rhs(1) = 0.0 do i = 2, nx-1 rhs(i) = 2.0 * un(i) - unm(i) end do c Zero boundaries rhs(nx) = 0.0 c Now use DGECO to factorize matrix and get condition number c Could use DGEFA to get just factorization... c unp passed in as workspace call DGECO(mat,nx, nx, ipvt, rcond, unp) if(trace)then write(0,*)'Reciprocal condition number = ',rcond end if if(1.0 + rcond .eq. 1.0)then write(0,*)'Matrix formed is singular!!! Exiting...' stop end if c Use DGESL to solve the linear system job = 0 call DGESL(mat, nx, nx, ipvt, rhs, job) c Copy the solution from rhs to unp do i = 1, nx unp(i) = rhs(i) end do return end c----------------------------------------------------------------------- c subroutine to form matrix for T-shaped stencil - reflecting boundaries subroutine form_mat1(nx, dt, dx, unm, un, mat, ipvt, temp) implicit none integer nx, ipvt(nx) real*8 dt, dx c Grid function declarations real*8 unm(nx), un(nx) c Matrix and LINPACK arrays real*8 mat(nx,nx), temp(nx) c Local variables integer i, j, job real*8 m, m2, fact, rcond logical trace trace = .TRUE. m = dt / dx m2= m**2 fact = 1.0 + 2.0 * m2 c write(0,*)m,m2,fact c First zero out matrix do i = 1, nx do j = 1, nx mat(i,j) = 0.0 end do end do c Form the matrix taking into account the boundary conditions c Here we take fixed boundaries. The fixed boundaries yield that c unp(1) = 0.0 and unp(nx) = 0.0 c We are left with a nx-2 by nx-2 matrix to solve for the interior points c First initialize row 1 and row nx mat(1,1) = 1.0 mat(2,2) = fact mat(2,3) = -m2 mat(nx-1,nx-1) = fact mat(nx-1,nx-2) = -m2 mat(nx,nx) = 1.0 c Set up rest of matrix do i = 3, nx-2 mat(i,i-1) = -m2 mat(i,i) = fact mat(i,i+1) = -m2 end do c Form right hand side of matrix equation c Don't forget that we require the interior only c Zero boundaries temp(1) = 0.0 do i = 2, nx-1 temp(i) = 2.0 * un(i) - unm(i) end do c Zero boundaries temp(nx) = 0.0 c Now use DGECO to factorize matrix and get condition number c Could use DGEFA to get just factorization... c unp passed in as workspace call DGECO(mat,nx, nx, ipvt, rcond, temp) if(trace)then write(0,*)'Reciprocal condition number = ',rcond end if if(1.0 + rcond .eq. 1.0)then write(0,*)'Matrix formed is singular!!! Exiting...' stop end if return end c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c subroutine to set up right-hand-side and solve for n+1 c for T stencil subroutine solve1(nx, dt, dx, unp, unm, un, mat, ipvt, rhs) implicit none integer nx, ipvt(nx) real*8 dt, dx c Grid function declarations real*8 unp(nx), unm(nx), un(nx) c Matrix and LINPACK arrays real*8 mat(nx,nx), rhs(nx) integer i, j, job real*8 m, m2, fact, rcond logical trace c Form right hand side of matrix equation c Don't forget that we require the interior only c boundaries rhs(1) = 0.0 do i = 2, nx-1 rhs(i) = 2.0 * un(i) - unm(i) end do c boundaries rhs(nx) = 0.0 c Use DGESL to solve the linear system job = 0 call DGESL(mat, nx, nx, ipvt, rhs, job) c Copy the solution from rhs to unp do i = 1, nx unp(i) = rhs(i) end do return end c----------------------------------------------------------------------- c subroutine to form matrix for 9 point stencil - reflecting boundaries subroutine form_mat2(nx, dt, dx, unm, un, mat, ipvt, temp) implicit none integer nx, ipvt(nx) real*8 dt, dx c Grid function declarations real*8 unm(nx), un(nx) c Matrix and LINPACK arrays real*8 mat(nx,nx), temp(nx) c Local variables integer i, j, job real*8 m, m2, a, rcond, thetax, thetat logical trace trace = .TRUE. thetax = 0.5 thetat = 0.5 m = dt / dx m2= m**2 a = 0.5 * ( thetat - m2 * thetax) c First zero out matrix do i = 1, nx do j = 1, nx mat(i,j) = 0.0 end do end do c Form the matrix taking into account the boundary conditions c Here we take fixed boundaries. The fixed boundaries yield that c unp(1) = 0.0 and unp(nx) = 0.0 c We are left with a nx-2 by nx-2 matrix to solve for the interior points c First initialize row 1 and row nx mat(1,1) = 1.0 mat(2,2) = 1.0 - 2.0 * a mat(2,3) = a mat(nx-1,nx-1) = 1.0 - 2.0 * a mat(nx-1,nx-2) = a mat(nx,nx) = 1.0 c Set up rest of matrix do i = 3, nx-2 mat(i,i-1) = a mat(i,i) = 1.0 - 2.0 * a mat(i,i+1) = a end do c Form right hand side of matrix equation c Don't forget that we require the interior only c Now use DGECO to factorize matrix and get condition number c Could use DGEFA to get just factorization... c temp passed in as workspace call DGECO(mat,nx, nx, ipvt, rcond, temp) if(trace)then write(0,*)'Reciprocal condition number = ',rcond end if if(1.0 + rcond .eq. 1.0)then write(0,*)'Matrix formed is singular!!! Exiting...' stop end if return end c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c subroutine to set up right-hand-side and solve for n+1 c for nine point stencil subroutine solve2(nx, dt, dx, unp, unm, un, mat, ipvt, rhs) implicit none integer nx, ipvt(nx) real*8 dt, dx c Grid function declarations real*8 unp(nx), unm(nx), un(nx) c Matrix and LINPACK arrays real*8 mat(nx,nx), rhs(nx) integer i, j, job real*8 m, m2, a, rcond, thetax, thetat logical trace trace = .TRUE. thetax = 0.5 thetat = 0.5 m = dt / dx m2= m**2 a = 0.5 * ( thetat - m2 * thetax) c Form right hand side of matrix equation c Don't forget that we require the interior only c boundaries rhs(1) = 0.0 do i = 2, nx-1 rhs(i) = -a*unm(i+1) + (2.0 * a - 1.0) * unm(i) . -a*unm(i-1) + (2.0 * a + m2 ) * un(i+1) . +2.0 * ( - 2.0 * a - m2 + 1.0)*un(i) . + ( 2.0 * a + m2) * un(i-1) end do c boundaries rhs(nx) = 0.0 c Use DGESL to solve the linear system job = 0 call DGESL(mat, nx, nx, ipvt, rhs, job) c Copy the solution from rhs to unp do i = 1, nx unp(i) = rhs(i) end do return end c-----------------------------------------------------------------------