Dear Alexander,
I checked the accuracy of d_Helmholtz_3D and it was much much larger than 1E-16. Since the solver uses a standard seven-point discretization, I verified the accuracy by the following code. (I have a uniform mesh for my problem.)
[cpp]for ( k=1; kGetTuple( i + j *NX + k *NX*NY);
phi_im1jk = Phi->GetTuple( (i-1) + j *NX + k *NX*NY);
phi_ip1jk = Phi->GetTuple( (i+1) + j *NX + k *NX*NY);
phi_ijm1k = Phi->GetTuple( i + (j-1)*NX + k *NX*NY);
phi_ijp1k = Phi->GetTuple( i + (j+1)*NX + k *NX*NY);
phi_ijkm1 = Phi->GetTuple( i + j *NX + (k-1)*NX*NY);
phi_ijkp1 = Phi->GetTuple( i + j *NX + (k+1)*NX*NY);
X_i = X->GetTuple( i );
X_im1 = X->GetTuple( i-1 );
Y_j = Y->GetTuple( j );
Y_jm1 = Y->GetTuple( j-1 );
Z_k = Z->GetTuple( k );
Z_km1 = Z->GetTuple( k-1 );
rhs= f->GetTuple( i + j*NX+ k*NX*NY);
res = ( phi_im1jk + phi_ip1jk - 2*phi_ijk)/pow( X_i - X_im1, 2 ) +
( phi_ijm1k + phi_ijp1k - 2*phi_ijk)/pow( Y_j - Y_jm1, 2 ) +
( phi_ijkm1 + phi_ijkp1 - 2*phi_ijk)/pow( Z_k - Z_km1, 2 ) +
rhs;
}
}
}[/cpp]
When I print res, the residual is about 1E-1. Is there something that I have to be careful when using the function? I need to have an accuracy about 1E-16. Please advise.
Thanks,
Ahmad
