1 |
C $Id: ini_rhs.F,v 1.1 2006/05/12 21:58:06 ce107 Exp $ |
2 |
CStartOfInterface |
3 |
SUBROUTINE INI_RHS |
4 |
C /==========================================================\ |
5 |
C | SUBROUTINE INI_RHS | |
6 |
C | o Initialise 2d conjugate gradient solver right-hand side| |
7 |
C |==========================================================| |
8 |
C | Set a source term b in Ax = b (1) | |
9 |
C | We solve (1) with neumann bc's whic means that b must | |
10 |
C | integrate out to zero over the whole domain. If b does | |
11 |
C | not integrate out to zero then the solution will | |
12 |
C | converge, but to a non-zero final residual. | |
13 |
C \==========================================================/ |
14 |
IMPLICIT NONE |
15 |
|
16 |
C === Global variables === |
17 |
#include "SIZE.h" |
18 |
#include "EEPARAMS.h" |
19 |
#include "PARAMS.h" |
20 |
#include "CG2D.h" |
21 |
|
22 |
C === Routine arguments === |
23 |
CEndOFInterface |
24 |
|
25 |
C === Local variables === |
26 |
C iG, jG - Global coordinate index |
27 |
C faceArea - Temporary used to hold cell face areas. |
28 |
C I,J,K - Loop counters |
29 |
C iMid, jMid - Global coords of mid-point of model domain |
30 |
INTEGER I, J |
31 |
INTEGER iG, jG |
32 |
INTEGER iMid, jMid |
33 |
INTEGER iQ, i3Q |
34 |
|
35 |
C-- Get model global domain mid-point |
36 |
iMid = Nx/2 |
37 |
iQ = Nx/4 |
38 |
i3Q = 3*Nx/4 |
39 |
jMid = Ny/2 |
40 |
|
41 |
C-- Set dummy source term for elliptic equation |
42 |
DO J=1-OLy,sNy+OLy |
43 |
DO I=1-OLx,sNx+OLx |
44 |
cg2d_b(I,J) = 0. _d 0 |
45 |
cg2d_x(I,J) = 0. _d 0 |
46 |
ENDDO |
47 |
ENDDO |
48 |
DO J=1,sNy |
49 |
DO I=1,sNx |
50 |
cg2d_x(I,J) = 0. _d 0 |
51 |
cg2d_b(I,J) = 0. _d 0 |
52 |
C-- Set +/-1 source function in middle of domain. |
53 |
iG = myXGlobalLo + I - 1 |
54 |
jG = myYGlobalLo + J - 1 |
55 |
IF ( iG .EQ. 1 .AND. jG .EQ. 1 ) cg2d_b(I,J) = 1. _d 0 |
56 |
IF ( iG .EQ. Nx .AND. jG .EQ. Ny ) cg2d_b(I,J) = -1. _d 0 |
57 |
ENDDO |
58 |
ENDDO |
59 |
C-- Update overlap regions synchronously |
60 |
CALL EXCH_XY_R8(cg2d_b) |
61 |
CALL EXCH_XY_R8(cg2d_x) |
62 |
C |
63 |
RETURN |
64 |
END |
65 |
|