kNaN: E:/turkuaz/krylov/codes/newton/1.2b/physics/bratu/anaJ/jf.f90 Source File

00001 !bratu problem
00002 subroutine jf(X,JAC,n)
00003 use localParameters
00004 implicit none
00005 integer i,j,k,n,z
00006 integer index
00007 real(8) X(n),JAC(n,n)
00008 real(8) dx,dy
00009 
00010 ! call index from the index matrix and mesh info!!!!
00011 
00012 
00013 dx=(Lx)/(nx-1.d0)
00014 dy=(Ly)/(ny-1.d0)
00015 
00016 
00017 
00018 
00019 JAC=0.d0
00020 
00021 
00022 
00023 
00024 !bottom boundary 
00025 
00026 j=1     
00027 z=(j-1)*nx
00028 do i=1,nx
00029         index=(z+i-1)*unk+1
00030         JAC(index,index)=1.d0
00031 end do
00032 
00033 !top boundary 
00034 
00035 j=ny     
00036 z=(j-1)*nx
00037 do i=1,nx
00038         index=(z+i-1)*unk+1
00039         JAC(index,index)=1.d0
00040 end do
00041 
00042 
00043 !left boundary 
00044 
00045 i=1     
00046 
00047 do j=2,ny-1
00048         index=((j-1)*nx+i-1)*unk+1
00049         JAC(index,index)=1.d0
00050 end do
00051 
00052 !right boundary 
00053 
00054 i=nx     
00055 
00056 do j=2,ny-1
00057         index=((j-1)*nx+i-1)*unk+1
00058         JAC(index,index)=1.d0
00059 end do
00060 
00061 
00062 ! interior
00063 
00064 
00065 do j=2,ny-1
00066         
00067         z=(j-1)*nx
00068 
00069         do i=2,ny-1
00070                 
00071                 index=(z+i-1)*unk+1
00072 
00073 
00074 
00075                 JAC(index,index)=-2.d0/dx**2.d0-2.d0/dy**2.d0+lambda*dexp(X(index))
00076                 JAC(index,index-unk)=1.d0/dx**2.d0
00077                 JAC(index,index+unk)=1.d0/dx**2.d0
00078                 JAC(index,index+unk*nx)=1.d0/dy**2.d0
00079                 JAC(index,index-unk*nx)=1.d0/dy**2.d0
00080 
00081 
00082         end do
00083 
00084 
00085 
00086 end do
00087 
00088 
00089 end subroutine jf