I did make the following changes. 1: I saved the vector f as a 1-dimensional array instead of a 2D one as given in the example 2: In the example code, the following lines are used. xi=hx*(ix-1)/lx and yi=hy*(iy-1)/ly Since ix is an integer, xi will always be rounded of to the nearest integer although it is defined as a double precision. Hence instead I used: xi=hx*real(ix-1)/lx and yi=hy*real(iy-1)/ly 3: In the example program q is set to 1. However to solve the Poisson equation q should be 0.0d0
The core program can be found below When applying Dirichlet boundary conditions everywhere, I manage to obtain the correct answer, however when applying Von Neumann boundary conditions, even on one side, the error is very large (maximum value for the error is 0.89. The core program can be found below. Does anyone have an idea where I went wrong?
The core program is now ax=0.0d0 bx=1.0d0 ay=0.0d0 by=1.0d0 nx=gridx !gridx=10 ny=gridz !gridz=10 q=0.0d0 ipar(:)=0 dpar(:)=0.0d0 spar(:)=0.0d0
do iy=1,ny+1 do ix=1,nx+1 xi=hx*real(ix-1)/lx yi=hy*real(iy-1)/ly cx=sin(2.0d0*pi*xi) cy=sin(2.0d0*pi*yi) u(ix+(iy-1)*(nx+1))=cx*cy f(ix+(iy-1)*(nx+1))=((8.0d0*pi*pi))*u(ix+(iy-1)*(nx+1)) u(ix+(iy-1)*(nx+1))=u(ix+(iy-1)*(nx+1))+1.0d0 end do end do
!do iy=1,nx+1 !bd_ax=-2.0d0*pi*sin(2.0d0*pi*real(iy-1)/ny) !bd_bx=2.0d0*pi*sin(2.0d0*pi*real(iy-1)/ny) !end do do ix=1,nx+1 bd_ay=-2.0d0*pi*sin(2.0d0*pi*real(ix-1)/nx) bd_by=2.0d0*pi*sin(2.0d0*pi*real(ix-1)/nx) end do
! Initializing simple data structures of Poisson Library for 2D Helmholtz Solver call d_init_Helmholtz_2D(ax, bx, ay, by, nx, ny, BCtype, q, ipar, dpar, stat) if(stat.ne.0)then write(*,*)'Double precision 2D Helmholtz example FAILED to compute the solution...' end if
! Computing the approximate solution of 2D Helmholtz problem ! between the Commit step and the subsequent call to the Solver routine! ! Otherwise the results may be wrong. call d_commit_Helmholtz_2D(f, bd_ax, bd_bx, bd_ay, bd_by, xhandle, ipar, dpar, stat) if(stat.ne.0)then write(*,*)'Double precision 2D Helmholtz example FAILED to compute the solution...' end if
! NOTE: Boundary data stored in the arrays bd_ax, bd_bx, bd_ay, bd_by should not be changed ! between the Commit step and the subsequent call to the Solver routine! call d_Helmholtz_2D(f, bd_ax, bd_bx, bd_ay, bd_by, xhandle, ipar, dpar, stat) if(stat.ne.0)then write(*,*)'Double precision 2D Helmholtz example FAILED to compute the solution...' end if
c1=0.0d0 do iy=1,ny+1 write(*,*) (ix-1)*hx,(iy-1)*hy,f(ix+(iy-1)*(ny+1))-u(ix+(iy-1)*(ny+1)) end do
do i=1,gridx+1 do j=1,gridz+1 write(*,*)i,j, 'ans',f((i)+(j-1)*(gridx+1)),u((i)+(j-1)*(gridx+1)) end do end do
I do now still have a problem when running the code with all Von Neumann boundary conditions. When I attempt this I get the following MKL Poisson Library warning: The problem is degenerate up to rounding errors! the approximate solution that provides the mimimal Euclidian norm or the solution will be computed...
An error is expected since an infinite number of solutions is possible when using 4 Von Neumann boundary conditions, however the gradients of the numerical solution should match the gradients found with the analytical solution. This seems not to be the case at the bottom left and bottom right corner of the computational domain. Everywhere else, the solution is OK. Do you have any idea how this larger error in the lower corner points (Error value of 0.18 on a 10x10 grid) is initiated and is there a way for me to prevent it.
When you obtain a finite-difference solution for a 2-D problem with discontinuities in the boundary conditions, you should go in with reasonable expectations. Typically, discontinuities are smoothed over a small region. As the mesh is made finer, this smoothed region becomes smaller. Just as an experiment cannot be set up with discontinuous boundary values, a finite-difference approximation will be inaccurate in regions of high (or infinite) gradients.
If you want accurate solutions near a discontinuity, you should consider alternative methods of solution.
It is often troublesome when contrived problems are solved numerically whose virtue is that a simple analytical solution is possible. The numerical solution should not be expected to have the same properties as a Fourier series solution (for example, convergence to the mean at a point of discontinuity).
In case when you use only zero Neumann boundary condition the kernel space is not empty. Indifferential equation it is mean that gradients of all solutions are equal but it is not clear how to calculate gradient of grid solution. You can calculate numeriacal aproximation of gradients in each internal point of solution but what about boundary point? Could you explain how you calculate gradient of solution in corner points?