C $Id: ini_rhs.F,v 1.2 2006/05/12 22:34:55 ce107 Exp $ CStartOfInterface SUBROUTINE INI_RHS C /==========================================================\ C | SUBROUTINE INI_RHS | C | o Initialise 2d conjugate gradient solver right-hand side| C |==========================================================| C | Set a source term b in Ax = b (1) | C | We solve (1) with neumann bc's whic means that b must | C | integrate out to zero over the whole domain. If b does | C | not integrate out to zero then the solution will | C | converge, but to a non-zero final residual. | C \==========================================================/ IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "CG2D.h" C === Routine arguments === CEndOFInterface C === Local variables === C iG, jG - Global coordinate index C faceArea - Temporary used to hold cell face areas. C I,J,K - Loop counters C iMid, jMid - Global coords of mid-point of model domain INTEGER I, J INTEGER iG, jG INTEGER iMid, jMid INTEGER iQ, i3Q C-- Get model global domain mid-point iMid = Nx/2 iQ = Nx/4 i3Q = 3*Nx/4 jMid = Ny/2 C-- Set dummy source term for elliptic equation DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx cg2d_b(I,J) = 0. _d 0 cg2d_x(I,J) = 0. _d 0 ENDDO ENDDO DO J=1,sNy DO I=1,sNx cg2d_x(I,J) = 0. _d 0 cg2d_b(I,J) = 0. _d 0 C-- Set +/-1 source function in middle of domain. iG = myXGlobalLo + I - 1 jG = myYGlobalLo + J - 1 IF ( iG .EQ. 1 .AND. jG .EQ. 1 ) cg2d_b(I,J) = 1. _d 0 IF ( iG .EQ. Nx .AND. jG .EQ. Ny ) cg2d_b(I,J) = -1. _d 0 ENDDO ENDDO C-- Update overlap regions synchronously CALL EXCH_XY_R8(cg2d_b) CALL EXCH_XY_R8(cg2d_x) C RETURN END