woodward Posted August 25, 2010 Report Share Posted August 25, 2010 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)CCC THIS SUBROUTINE COMPUTES THE NEW ISOPARAMETRIC COORDINATESC 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 NEWC CONCENTRATION IS COMPUTED.CC 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 = 0C DO WHILE ((ISTOP .EQ. 0) .AND. (K .LT. 100))C K = K + 1 RESS = 0.0D0CC COMPUTE NEW ISOPARAMETRIC COORDINATES.C DO J = 1, 3C 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 DOC 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 + RESJC END DOCC CHECK FOR CONVERGENCE.C IF (RESS .LT. 1.0D-4) THEN ISTOP = 1 END IFC END DOCC IF THERE IS A VALID SOLUTION, COMPUTE THE NEWC CONCENTRATION.C IF ((DABS (SNEW(1)) .LE. 1.001D0) .AND. & (DABS (SNEW(2)) .LE. 1.001D0) .AND. & (DABS (SNEW(3)) .LE. 1.001D0)) THENC FLAG(N) = 4.0D0 Quote Link to comment Share on other sites More sharing options...
Recommended Posts
Join the conversation
You can post now and register later. If you have an account, sign in now to post with your account.