Jump to content
GMS, SMS, and WMS User Forum

Help with MMOC (Modified method of characteristics) code please


woodward
 Share

Recommended Posts

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

Link to comment
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Loading...
 Share

×
×
  • Create New...