Help with MMOC (Modified method of characteristics) code please

[woodward]

I am having trouble with a bug in FEMWATER (not in GMS) and was hoping someone here might be able to help me.

The bug is in the modified method of characteristics lagrangian concentration backtracking code for hexahedral elements only.

The FEMWATER code fragment below is supposed to figure out whether cartesian point XXQ(3) lies within hexahedron XX(8,3). It does this by transforming XX to a reference hexahedron (cube) SI(8,3), and then seeing if the transformed location of XXQ, which is called SNEW(3), lies inside SI. It is using some kind of iterative method to calculate SNEW. But it doesn't always work. Sometimes on the first pass through SUM2=0 and this causes a DIV/0 error. Does anyone know this method, or can you see a bug? Thanks so much!

SUBROUTINE BRKINS (N, XX, XXQ, CC, CS, FLAG, OUTEL)

C

C

C THIS SUBROUTINE COMPUTES THE NEW ISOPARAMETRIC COORDINATES

C INSIDE THE GIVEN HEXAHEDRON FOR THE GIVEN POINT (ZQ, YQ, ZQ).

C IF THEY ARE BETWEEN -1 AND 1, THE POINT IS INSIDE, AND A NEW

C CONCENTRATION IS COMPUTED.

C

C

IMPLICIT REAL * 8 (A-H, O-Z)

C

INCLUDE 'gwpara.inc'

C

DIMENSION CS(MAXNPK), FLAG(MAXNPK)

DIMENSION OUTEL(MAXNPK)

C

DIMENSION XX(8, 3), XXQ(3), CC(8)

DIMENSION SI(8, 3), SOLD(3), SNEW(3), IORDER(2, 3)

DATA SI /

& -1.0D0, 1.0D0, 1.0D0, -1.0D0, -1.0D0, 1.0D0, 1.0D0, -1.0D0,

& -1.0D0, -1.0D0, 1.0D0, 1.0D0, -1.0D0, -1.0D0, 1.0D0, 1.0D0,

& -1.0D0, -1.0D0, -1.0D0, -1.0D0, 1.0D0, 1.0D0, 1.0D0, 1.0D0 /

DATA IORDER / 2, 3, 1, 3, 1, 2 /

C

DO J = 1, 3

SNEW(J) = 0.0D0

END DO

K = 0

ISTOP = 0

C

DO WHILE ((ISTOP .EQ. 0) .AND. (K .LT. 100))

C

K = K + 1

RESS = 0.0D0

C

C COMPUTE NEW ISOPARAMETRIC COORDINATES.

C

DO J = 1, 3

C

SOLD(J) = SNEW(J)

JJ1 = IORDER(1, J)

JJ2 = IORDER(2, J)

C

SUM1 = 0.0D0

SUM2 = 0.0D0

DO I = 1, 8

T1 = (1.0D0 + SI(I, JJ1) * SNEW(JJ1)) * (1.0D0 +

& SI(I, JJ2) * SNEW(JJ2)) * XX(I, J)

SUM1 = SUM1 + T1

SUM2 = SI(I, J) * T1 + SUM2

END DO

C

SNEW(J) = (XXQ(J) * 8.0D0 - SUM1) / SUM2

RESJ = DABS (SNEW(J) - SOLD(J))

DABSJ = DABS (SNEW(J))

IF (DABSJ .GT. 1.0D0) RESJ = RESJ / DABSJ

RESS = RESS + RESJ

C

END DO

C

C CHECK FOR CONVERGENCE.

C

IF (RESS .LT. 1.0D-4) THEN

ISTOP = 1

END IF

C

END DO

C

C IF THERE IS A VALID SOLUTION, COMPUTE THE NEW

C CONCENTRATION.

C

IF ((DABS (SNEW(1)) .LE. 1.001D0) .AND.

& (DABS (SNEW(2)) .LE. 1.001D0) .AND.

& (DABS (SNEW(3)) .LE. 1.001D0)) THEN

C

FLAG(N) = 4.0D0