!/*****************************************************************************/ ! * ! * Elmer/Ice, a glaciological add-on to Elmer ! * http://elmerice.elmerfem.org ! * ! * ! * This program is free software; you can redistribute it and/or ! * modify it under the terms of the GNU General Public License ! * as published by the Free Software Foundation; either version 2 ! * of the License, or (at your option) any later version. ! * ! * This program is distributed in the hope that it will be useful, ! * but WITHOUT ANY WARRANTY; without even the implied warranty of ! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! * GNU General Public License for more details. ! * ! * You should have received a copy of the GNU General Public License ! * along with this program (in file fem/GPL-2); if not, write to the ! * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, ! * Boston, MA 02110-1301, USA. ! * ! *****************************************************************************/ ! ****************************************************************************** ! * ! * Authors: Joe Todd ! * Email: ! * Web: http://elmerice.elmerfem.org ! * ! * ! ***************************************************************************** !A routine for getting frontal advance in 3D, given velocity and melt SUBROUTINE FrontAdvance3D ( Model, Solver, dt, TransientSimulation ) USE CalvingGeometry USE DefUtils IMPLICIT NONE !----------------------------------------------- TYPE(Model_t) :: Model TYPE(Solver_t), TARGET :: Solver REAL(KIND=dp) :: dt LOGICAL :: TransientSimulation !----------------------------------------------- TYPE(Mesh_t), POINTER :: Mesh TYPE(Nodes_t), TARGET :: FrontNodes, ColumnNodes TYPE(ValueList_t), POINTER :: Params TYPE(Variable_t), POINTER :: Var, VeloVar, MeltVar, NormalVar, TangledVar INTEGER :: i, j, k, m, n, DOFs, TotalNodes, NodesPerLevel, ExtrudedLevels,& FaceNodeCount, DummyInt, col, hits, ierr, FrontLineCount, county,& ShiftIdx, ShiftToIdx, NoTangledGroups, PivotIdx, CornerIdx, & SecondIdx, FirstTangleIdx, LastTangleIdx INTEGER, POINTER :: Perm(:), FrontPerm(:)=>NULL(), TopPerm(:)=>NULL(), & FrontNodeNums(:)=>NULL() INTEGER, ALLOCATABLE :: FNColumns(:), FrontColumnList(:), FrontLocalNodeNumbers(:), & NodeNumbers(:), TangledGroup(:), TangledPivotIdx(:), UpdatedDirection(:) REAL(KIND=dp) :: FrontOrientation(3),RotationMatrix(3,3),& NodeVelo(3), NodeMelt(3), NodeNormal(3), ColumnNormal(3), CornerDirection(2),& MeltRate, Displace(3), NodeHolder(3), DangerGrad, ShiftTo, direction, & ShiftDist, y_coord(2), epsShift, ShiftToY, LongRangeLimit, MaxDisplacement, LimitZ, & p1(2),p2(2),q1(2),q2(2),intersect(2), LeftY, RightY, EpsTangle,thisEps,Shift, thisY REAL(KIND=dp), POINTER :: Advance(:) REAL(KIND=dp), ALLOCATABLE :: Rot_y_coords(:,:), Rot_z_coords(:,:), ColumnNormals(:,:), & TangledShiftTo(:) LOGICAL :: Found, Debug, Parallel, Boss, ShiftLeft, LeftToRight, MovedOne, ShiftSecond, & Protrusion, SqueezeLeft, SqueezeRight, FirstTime=.TRUE., intersect_flag, FrontMelting, & IgnoreVelo LOGICAL, ALLOCATABLE :: DangerZone(:), WorkLogical(:), UpdatedColumn(:),& Tangled(:), DontMove(:) CHARACTER(LEN=MAX_NAME_LEN) :: SolverName, VeloVarName, MeltVarName, & NormalVarName, FrontMaskName, TopMaskName, TangledVarName SAVE :: FirstTime !----------------------------------------------- ! Initialisation !----------------------------------------------- Debug = .FALSE. Parallel = (ParEnv % PEs > 1) Boss = ParEnv % MyPE == 0 SolverName = "FrontAdvance3D" Params => Solver % Values Mesh => Solver % Mesh !The main solver var contains the magnitude of front advance Var => Solver % Variable Advance => Var % Values Perm => Var % Perm DOFs = Var % DOFs IF(Var % DOFs /= 3) CALL Fatal(SolverName, "Variable should have 3 DOFs...") IgnoreVelo = ListGetLogical(Params, "Ignore Velocity", Found) IF(.NOT. Found) THEN IgnoreVelo = .FALSE. ELSE CALL Info(SolverName, "Ignoring velocity (melt undercutting only)") END IF IF(.NOT. IgnoreVelo) THEN !Get the flow solution VeloVarName = ListGetString(Params, "Flow Solution Variable Name", Found) IF(.NOT. Found) THEN CALL Info(SolverName, "Flow Solution Variable Name not found, assuming 'Flow Solution'") VeloVarName = "Flow Solution" END IF VeloVar => VariableGet(Mesh % Variables, VeloVarName, .TRUE., UnfoundFatal=.TRUE.) END IF !Get melt rate MeltVarName = ListGetString(Params, "Melt Variable Name", Found) IF(.NOT. Found) THEN CALL Info(SolverName, "Melt Variable Name not found, assuming no frontal melting") FrontMelting = .FALSE. ELSE FrontMelting = .TRUE. MeltVar => VariableGet(Mesh % Variables, MeltVarName, .TRUE., UnfoundFatal=.TRUE.) END IF TangledVarName = "Tangled" TangledVar => VariableGet(Mesh % Variables, TangledVarName, .TRUE., UnfoundFatal=.TRUE.) TangledVar % Values = 0.0_dp !Get front normal vector NormalVarName = ListGetString(Params, "Normal Vector Variable Name", UnfoundFatal=.TRUE.) NormalVar => VariableGet(Mesh % Variables, NormalVarName, .TRUE., UnfoundFatal=.TRUE.) !Get the orientation of the calving front, compute rotation matrix !TODO: generalize and link FrontOrientation = GetFrontOrientation(Model) RotationMatrix = ComputeRotationMatrix(FrontOrientation) DangerGrad = ListGetConstReal(Params, "Front Gradient Threshold", Found) IF(.NOT.Found) THEN CALL Info(SolverName, "'Front Gradient Threshold' not found, setting to 5.0") DangerGrad = 5.0 END IF !When correcting for unprojectability, to prevent interp problems, !ensure that it is SLIGHTLY beyond beyond parallel, by adding this !value (metres) to the displacement epsShift = ListGetConstReal(Params, "Front Projectability Epsilon", Found) IF(.NOT.Found) THEN CALL Info(SolverName, "'Front Projectability Epsilon' not found, setting to 1.0") epsShift = 1.0 END IF epsTangle = ListGetConstReal(Params, "Front Tangle Epsilon", Found) IF(.NOT.Found) THEN CALL Info(SolverName, "'Front Tangle Epsilon' not found, setting to 1.0") epsTangle = 1.0 END IF LongRangeLimit = ListGetConstReal(Params, "Column Max Longitudinal Range", Found) IF(.NOT.Found) THEN CALL Info(SolverName, "'Column Max Longitudinal Range' not found, setting to 300.0") LongRangeLimit = 300.0_dp END IF MaxDisplacement = ListGetConstReal(Params, "Maximum Node Displacement", Found) IF(.NOT.Found) THEN CALL Info(SolverName, "'Maximum Node Displacement' not found, setting to 1.0E4.") MaxDisplacement = 1.0E4_dp END IF !------------------------------------------- ! Find FNColumns CALL MPI_AllReduce(MAXVAL(Mesh % ParallelInfo % GlobalDOFs), TotalNodes, & 1, MPI_INTEGER, MPI_MAX, ELMER_COMM_WORLD,ierr) ExtrudedLevels = ListGetInteger(CurrentModel % Simulation,'Extruded Mesh Levels',Found) IF(.NOT. Found) ExtrudedLevels = & ListGetInteger(CurrentModel % Simulation,'Remesh Extruded Mesh Levels',Found) IF(.NOT. Found) CALL Fatal(SolverName,& "Unable to find 'Extruded Mesh Levels' or 'Remesh Extruded Mesh Levels'") NodesPerLevel = TotalNodes / ExtrudedLevels ALLOCATE(FNColumns(Mesh % NumberOfNodes)) FNColumns = MOD(Mesh % ParallelInfo % GlobalDOFs, NodesPerLevel) !Get the front line FrontMaskName = "Calving Front Mask" TopMaskName = "Top Surface Mask" CALL MakePermUsingMask( Model, Solver, Mesh, TopMaskName, & .FALSE., TopPerm, dummyint) CALL MakePermUsingMask( Model, Solver, Mesh, FrontMaskName, & .FALSE., FrontPerm, FaceNodeCount) CALL GetDomainEdge(Model, Mesh, TopPerm, FrontNodes, & FrontNodeNums, Parallel, FrontMaskName, Simplify=.FALSE.) !Pass FrontNodeNums to all CPUs IF(Boss) FrontLineCount = SIZE(FrontNodeNums) CALL MPI_BCAST(FrontLineCount , 1, MPI_INTEGER, 0, ELMER_COMM_WORLD, ierr) !Determine whether front columns are arranged !left to right, and reorder if not IF(Boss) THEN NodeHolder(1) = FrontNodes % x(1) NodeHolder(2) = FrontNodes % y(1) NodeHolder(3) = FrontNodes % z(1) NodeHolder = MATMUL(RotationMatrix, NodeHolder) y_coord(1) = NodeHolder(2) NodeHolder(1) = FrontNodes % x(FrontLineCount) NodeHolder(2) = FrontNodes % y(FrontLineCount) NodeHolder(3) = FrontNodes % z(FrontLineCount) NodeHolder = MATMUL(RotationMatrix, NodeHolder) y_coord(2) = NodeHolder(2) LeftToRight = y_coord(2) > y_coord(1) IF(.NOT. LeftToRight) THEN IF(Debug) PRINT *,'Debug, switching to LeftToRight' FrontNodeNums = FrontNodeNums(FrontLineCount:1:-1) FrontNodes % x = FrontNodes % x(FrontLineCount:1:-1) FrontNodes % y = FrontNodes % y(FrontLineCount:1:-1) FrontNodes % z = FrontNodes % z(FrontLineCount:1:-1) END IF END IF IF(.NOT. Boss) ALLOCATE(FrontNodeNums(FrontLineCount)) CALL MPI_BCAST(FrontNodeNums , FrontLineCount, MPI_INTEGER, 0, ELMER_COMM_WORLD, ierr) !FrontNodeNums is the GlobalDOFs of the nodes on the front !So this can be easily converted into a FNColumn like list ALLOCATE(FrontColumnList(FrontLineCount), & FrontLocalNodeNumbers(FrontLineCount), & DangerZone(FrontLineCount), & ColumnNormals(FrontLineCount,3)) FrontColumnList = MOD(FrontNodeNums, NodesPerLevel) FrontLocalNodeNumbers = -1 DangerZone = .FALSE. ColumnNormals = 0.0_dp DO i=1,FrontLineCount n = FrontNodeNums(i) DO j=1,Mesh % NumberOfNodes IF(Mesh % ParallelInfo % GlobalDOFs(j) == n) THEN FrontLocalNodeNumbers(i) = j END IF END DO END DO !-------------------------------------- ! Action: Compute lagrangian displacement for all nodes ! This is the main function of the Solver. ! Everything following this DO loop is taking care of ! unprojectability/high gradient etc !-------------------------------------- DO i=1,Mesh % NumberOfNodes IF(Perm(i) <= 0) CYCLE IF(FrontMelting) THEN IF(MeltVar % Perm(i) <= 0) & CALL Fatal(SolverName, "Permutation error on front node!") !Scalar melt value from Plume solver MeltRate = MeltVar % Values(MeltVar % Perm(i)) NodeNormal(1) = NormalVar % Values(((NormalVar % Perm(i)-1)*NormalVar % DOFs) + 1) NodeNormal(2) = NormalVar % Values(((NormalVar % Perm(i)-1)*NormalVar % DOFs) + 2) NodeNormal(3) = NormalVar % Values(((NormalVar % Perm(i)-1)*NormalVar % DOFs) + 3) NodeMelt = NodeNormal * MeltRate ELSE NodeMelt = 0.0_dp END IF IF(IgnoreVelo) THEN NodeVelo = 0.0 ELSE !Compute front normal component of velocity NodeVelo(1) = VeloVar % Values(((VeloVar % Perm(i)-1)*VeloVar % DOFs) + 1) NodeVelo(2) = VeloVar % Values(((VeloVar % Perm(i)-1)*VeloVar % DOFs) + 2) NodeVelo(3) = VeloVar % Values(((VeloVar % Perm(i)-1)*VeloVar % DOFs) + 3) END IF Displace = 0.0 Displace(1) = NodeVelo(1) - NodeMelt(1) Displace(2) = NodeVelo(2) - NodeMelt(2) Displace(3) = NodeVelo(3) - NodeMelt(3) Displace = Displace * dt IF(MAXVAL(Displace) > MaxDisplacement) THEN WRITE(Message,'(A,i0,A)') "Maximum allowable front displacement exceeded for node ",i,". Limiting..." CALL Warn(SolverName, Message) Displace = Displace * (MaxDisplacement/MAXVAL(Displace)) END IF Advance((Perm(i)-1)*DOFs + 1) = Displace(1) Advance((Perm(i)-1)*DOFs + 2) = Displace(2) Advance((Perm(i)-1)*DOFs + 3) = Displace(3) !First and last (i.e. lateral margin) columns don't move IF((FNColumns(i) == FrontColumnList(1)) .OR. & (FNColumns(i) == FrontColumnList(FrontLineCount))) THEN Advance((Perm(i)-1)*DOFs + 1) = 0.0_dp Advance((Perm(i)-1)*DOFs + 2) = 0.0_dp Advance((Perm(i)-1)*DOFs + 3) = 0.0_dp END IF END DO !---------------------------------------------------------- ! Now we need to look for and limit regions on the ! where high gradient + 'straightness' makes ! remeshing into vertical columns troublesome. !---------------------------------------------------------- ! ! Five things: ! - Limit horizontal range (otherwise, Remesh will smear a column of nodes ! all over the place) ! - Limit undercut range (longitudinal range) to prevent mesh degeneracy ! - Prevent retreat at node near lateral margins ! - Check for 'tangled' regions, where fixing unprojectability one way ! causes unprojectability the other way. ! - Prevent Unprojectability - requires knowledge of neighbouring columns !---------------------------------------------------------- !------------------------------------------- ! Cycle columns, looking at average normal for the column !------------------------------------------- ! Normal vector is limited by projectability, so if average is high, ! all are high ! ! Note that we are cycling the parallel gathered line, so each part ! will only have a few columns (if any) ! So, mark troublesome columns where present, then communicate !------------------------------------------- DO i=1,FrontLineCount col = FrontColumnList(i) hits = COUNT(FNColumns == col) IF(hits == 0) CYCLE ColumnNormal = 0.0_dp !Gather normal vector DO j=1,Mesh % NumberOfNodes IF(FNColumns(j) == col) THEN ColumnNormal(1) = ColumnNormal(1) + & NormalVar % Values(((NormalVar % Perm(j)-1)*NormalVar % DOFs) + 1) ColumnNormal(2) = ColumnNormal(2) + & NormalVar % Values(((NormalVar % Perm(j)-1)*NormalVar % DOFs) + 2) ColumnNormal(3) = ColumnNormal(3) + & NormalVar % Values(((NormalVar % Perm(j)-1)*NormalVar % DOFs) + 3) END IF END DO !unit vector ColumnNormal = ColumnNormal / hits !Convert to front coordinate system ColumnNormal = MATMUL( RotationMatrix, ColumnNormal ) !Save for later ColumnNormals(i,1:3) = ColumnNormal(1:3) !We're concerned with the lateral (rather than vertical) component IF( ABS(ColumnNormal(2) / ColumnNormal(3)) > DangerGrad ) THEN DangerZone(i) = .TRUE. END IF END DO !Communicate between partitions !This is an 'OR' condition insofar as only one partition !will have actually computed the gradient IF(Parallel) THEN DO i=1,FrontLineCount CALL SParIterAllReduceOR(DangerZone(i)) END DO END IF !Mark before and after dangerzone too ALLOCATE(WorkLogical(FrontLineCount)) WorkLogical = .FALSE. DO i=1,FrontLineCount IF(DangerZone(i)) WorkLogical(i) = .TRUE. IF(i > 1) THEN IF(DangerZone(i-1)) WorkLogical(i) = .TRUE. END IF IF(i < FrontLineCount) THEN IF(DangerZone(i+1)) WorkLogical(i) = .TRUE. END IF END DO DangerZone = WorkLogical DEALLOCATE(WorkLogical) !Find and store leftmost and rightmost (rotated) coordinate of !each column for checking projectability later !Rot_y_coords(:,1) is leftmost, and (:,2) is right most y coord of !each column's nodes, and Rot_z_coords(:,:) is the min(1)/max(2) z ALLOCATE(Rot_y_coords(FrontLineCount,2),& Rot_z_coords(FrontLineCount,2)) Rot_y_coords(:,1) = HUGE(0.0_dp) Rot_y_coords(:,2) = -HUGE(0.0_dp) Rot_z_coords(:,1) = HUGE(0.0_dp) Rot_z_coords(:,2) = -HUGE(0.0_dp) DO i=1,FrontLineCount col = FrontColumnList(i) IF(COUNT(FNColumns == col) == 0) CYCLE DO j=1,Mesh % NumberOfNodes IF(FNColumns(j) /= col) CYCLE NodeHolder(1) = Mesh % Nodes % x(j) + Advance((Perm(j)-1)*DOFs + 1) NodeHolder(2) = Mesh % Nodes % y(j) + Advance((Perm(j)-1)*DOFs + 2) NodeHolder(3) = Mesh % Nodes % z(j) + Advance((Perm(j)-1)*DOFs + 3) NodeHolder = MATMUL(RotationMatrix, NodeHolder) Rot_y_coords(i,1) = MIN(Rot_y_coords(i,1), NodeHolder(2)) Rot_y_coords(i,2) = MAX(Rot_y_coords(i,2), NodeHolder(2)) Rot_z_coords(i,1) = MIN(Rot_z_coords(i,1), NodeHolder(3)) Rot_z_coords(i,2) = MAX(Rot_z_coords(i,2), NodeHolder(3)) END DO END DO !----------------------------------- ! Limit lateral range (to 0) in DangerZone !----------------------------------- ! The physical interpretation of this is that, when melt undercutting occurs ! in a region of high gradient, the unmelted parts of that column also shift ! with the melt. ! This is a minor limitation, but better than Free Surface Equation! DO i=1,FrontLineCount IF(.NOT. DangerZone(i)) CYCLE col = FrontColumnList(i) hits = COUNT(FNColumns == col) IF(hits == 0) CYCLE ALLOCATE(NodeNumbers(hits), & ColumnNodes % x(hits),& ColumnNodes % y(hits),& ColumnNodes % z(hits)) ColumnNodes % NumberOfNodes = hits !Gather nodenumbers in column county = 0 DO j=1,Mesh % NumberOfNodes IF(FNColumns(j) /= col) CYCLE county = county + 1 NodeNumbers(county) = j ColumnNodes % x(county) = Mesh % Nodes % x(j) + Advance((Perm(j)-1)*DOFs + 1) ColumnNodes % y(county) = Mesh % Nodes % y(j) + Advance((Perm(j)-1)*DOFs + 2) ColumnNodes % z(county) = Mesh % Nodes % z(j) + Advance((Perm(j)-1)*DOFs + 3) END DO !direction - which way is this part of the front pointing? ShiftLeft = ColumnNormals(i,2) > 0 IF(ShiftLeft) THEN ShiftTo = HUGE(ShiftTo) ELSE ShiftTo = -HUGE(ShiftTo) END IF !Rotate points and find furthest left (or right) DO j=1,ColumnNodes % NumberOfNodes NodeHolder(1) = ColumnNodes % x(j) NodeHolder(2) = ColumnNodes % y(j) NodeHolder(3) = ColumnNodes % z(j) NodeHolder = MATMUL(RotationMatrix, NodeHolder) IF(ShiftLeft) THEN IF(NodeHolder(2) < ShiftTo) ShiftTo = NodeHolder(2) ELSE IF(NodeHolder(2) > ShiftTo) ShiftTo = NodeHolder(2) END IF END DO !Now, for each node in column, compute the displacement (in rotated y coordinate) DO j=1,ColumnNodes % NumberOfNodes NodeHolder(1) = ColumnNodes % x(j) NodeHolder(2) = ColumnNodes % y(j) NodeHolder(3) = ColumnNodes % z(j) NodeHolder = MATMUL(RotationMatrix, NodeHolder) ShiftDist = ShiftTo - NodeHolder(2) Displace = 0.0_dp Displace(2) = ShiftDist Displace = MATMUL(TRANSPOSE(RotationMatrix), Displace) !Adjust the variable values to shift nodes in line. Advance((Perm(NodeNumbers(j))-1)*DOFs + 1) = & Advance((Perm(NodeNumbers(j))-1)*DOFs + 1) + Displace(1) Advance((Perm(NodeNumbers(j))-1)*DOFs + 2) = & Advance((Perm(NodeNumbers(j))-1)*DOFs + 2) + Displace(2) IF(Debug) PRINT *,'node: ',NodeNumbers(j),' shifting: ',ShiftDist,' to ',ShiftTo,' point: ', & ColumnNodes % x(j),ColumnNodes % y(j),ColumnNodes % z(j) END DO DEALLOCATE(NodeNumbers, & ColumnNodes % x, & ColumnNodes % y, & ColumnNodes % z) END DO !----------------------------------- ! Limit longitudinal range everywhere !----------------------------------- ! Two options for how to do this: ! Drag nodes back with melt, or effectively ! limit melt (keeping nodes forward) ! Opt for the latter, or else we're ! forcing melting to have an effect !----------------------------------- DO i=1,FrontLineCount col = FrontColumnList(i) hits = COUNT(FNColumns == col) IF(hits == 0) CYCLE IF((Rot_z_coords(i,2) - Rot_z_coords(i,1)) > LongRangeLimit) THEN LimitZ = Rot_z_coords(i,2) - LongRangeLimit Rot_z_coords(i,1) = LimitZ DO j=1,Mesh % NumberOfNodes IF(FNColumns(j) /= col) CYCLE NodeHolder(1) = Mesh % Nodes % x(j) + Advance((Perm(j)-1)*DOFs + 1) NodeHolder(2) = Mesh % Nodes % y(j) + Advance((Perm(j)-1)*DOFs + 2) NodeHolder(3) = Mesh % Nodes % z(j) + Advance((Perm(j)-1)*DOFs + 3) NodeHolder = MATMUL(RotationMatrix, NodeHolder) IF(NodeHolder(3) >= LimitZ) CYCLE Displace = 0.0_dp Displace(3) = LimitZ - NodeHolder(3) Displace = MATMUL(TRANSPOSE(RotationMatrix), Displace) !Adjust the variable values to shift nodes in line. Advance((Perm(j)-1)*DOFs + 1) = & Advance((Perm(j)-1)*DOFs + 1) + Displace(1) Advance((Perm(j)-1)*DOFs + 2) = & Advance((Perm(j)-1)*DOFs + 2) + Displace(2) IF(Debug) PRINT *,ParEnv % MyPE,' Debug, limiting node range: col ',& i,' node ',j,' limit: ',LimitZ,' displace: ',Displace END DO END IF END DO !----------------------------------- ! Now check projectability !----------------------------------- ! If an unprojectable portion of the front is found, ! the inland nodes remain in place, while those further ! advanced shift sideways to maintain projectability ! ! Typical case is melt plume at end of recently calved ! straight edge. ! ! requires MPI comms !Gather rotated y coord of all columns DO i=1,FrontLineCount CALL MPI_AllReduce(MPI_IN_PLACE, Rot_y_coords(i,1), & 1, MPI_DOUBLE_PRECISION, MPI_MIN, ELMER_COMM_WORLD,ierr) CALL MPI_AllReduce(MPI_IN_PLACE, Rot_y_coords(i,2), & 1, MPI_DOUBLE_PRECISION, MPI_MAX, ELMER_COMM_WORLD,ierr) CALL MPI_AllReduce(MPI_IN_PLACE, Rot_z_coords(i,1), & 1, MPI_DOUBLE_PRECISION, MPI_MIN, ELMER_COMM_WORLD,ierr) CALL MPI_AllReduce(MPI_IN_PLACE, Rot_z_coords(i,2), & 1, MPI_DOUBLE_PRECISION, MPI_MAX, ELMER_COMM_WORLD,ierr) IF(Boss .AND. Debug) PRINT *,'Debug, rot_y_coords: ',i,rot_y_coords(i,:) IF(Boss .AND. Debug) PRINT *,'Debug, rot_z_coords: ',i,rot_z_coords(i,:) END DO ALLOCATE(UpdatedColumn(FrontLineCount), & UpdatedDirection(FrontLineCount), & Tangled(FrontLineCount),& TangledGroup(FrontLineCount),& TangledPivotIdx(FrontLineCount),& TangledShiftTo(FrontLineCount),& DontMove(FrontLineCount)) UpdatedColumn = .FALSE. UpdatedDirection = 0 Tangled = .FALSE. TangledGroup = 0 TangledPivotIdx = 0 TangledShiftTo = 0.0_dp DontMove = .FALSE. DontMove(1) = .TRUE. DontMove(FrontLineCount) = .TRUE. !TODO: Deal with DontMove and tangling properly... ! Ensure that comms is not necessary for DontMove !Test for thin sliver at lateral margins - potential to get !stuck like this, so impose an angle here !We require that the 2nd node from the end not retreat beyond !a 45 degree angle inland. PYTHAGORAS !Shift the 2nd column in the lateral direction only, !to the point where it makes a 45 degree angle w/ !the corner. ! * ! / | ! / | ! * <- * DO i=1,2 IF(i==1) THEN CornerIdx = 1 SecondIdx = 2 direction = 1.0_dp ELSE CornerIdx = FrontLineCount SecondIdx = FrontLineCount-1 direction = -1.0_dp END IF !Conveniently, first time round we want the left most coord (Rot_y_coords(:,1)) !whereas second time we want the right (Rot_y_cords(:,2)) CornerDirection(1) = Rot_y_coords(SecondIdx,i) - Rot_y_coords(CornerIdx,1) CornerDirection(2) = Rot_z_coords(SecondIdx,i) - Rot_z_coords(CornerIdx,1) !Unit vector CornerDirection = CornerDirection / (SUM(CornerDirection ** 2.0)**0.5) IF(Debug) PRINT *,ParEnv % MyPE, 'CornerDirection ',i,': ',CornerDirection IF( CornerDirection(2) < (-1.0_dp / (2.0_dp ** 0.5)) ) THEN !Shift 2nd node and mark DontMove ShiftToY = Rot_y_coords(CornerIdx,1) + & direction * ((Rot_z_coords(CornerIdx,1) - Rot_z_coords(SecondIdx,1))) Rot_y_coords(SecondIdx,1) = ShiftToY Rot_y_coords(SecondIdx,2) = ShiftToY DontMove(SecondIdx) = .TRUE. IF(Debug) PRINT *,ParEnv % MyPE,' debug, side ',i,' of front is too inwards' DO k=1,Mesh % NumberOfNodes IF(FNColumns(k) /= FrontColumnList(SecondIdx)) CYCLE NodeHolder(1) = Mesh % Nodes % x(k) + Advance((Perm(k)-1)*DOFs + 1) NodeHolder(2) = Mesh % Nodes % y(k) + Advance((Perm(k)-1)*DOFs + 2) NodeHolder(3) = Mesh % Nodes % z(k) + Advance((Perm(k)-1)*DOFs + 3) NodeHolder = MATMUL(RotationMatrix, NodeHolder) Displace = 0.0_dp Displace(2) = ShiftToY - NodeHolder(2) Displace = MATMUL(TRANSPOSE(RotationMatrix), Displace) IF(Debug) PRINT *,ParEnv % MyPE, ' Debug, shifting next to corner node ',k,& NodeHolder,' by ',Displace !Adjust the variable values to shift nodes in line. Advance((Perm(k)-1)*DOFs + 1) = Advance((Perm(k)-1)*DOFs + 1) + Displace(1) Advance((Perm(k)-1)*DOFs + 2) = Advance((Perm(k)-1)*DOFs + 2) + Displace(2) END DO END IF END DO !----------------------------------------------- ! Cycle the line, looking for unprojectable ! (but not tangled) regions. ! e.g. plume melting behind a vertical portion !----------------------------------------------- !In case shifting nodes affects another partition, we loop so long as !at least one partition has MovedOne MovedOne = .TRUE. county = 0 DO WHILE(MovedOne) county = county + 1 IF(county > 100) CALL Fatal(SolverName, "Infinite loop!") IF(Debug) PRINT *,ParEnv % MyPE, 'Debug, iterating projectability ', county MovedOne = .FALSE. UpdatedColumn = .FALSE. DO i=2,FrontLineCount !If already detangled, skip !IF(Tangled(i)) CYCLE !epsShift here IF((Rot_y_coords(i,1) - Rot_y_coords(i-1,2)) < 0.9*epsShift) THEN IF(Debug) PRINT *, 'Debug diff ',i,i-1, & (Rot_y_coords(i,1) - Rot_y_coords(i-1,2)), epsShift IF(DontMove(i) .AND. DontMove(i-1)) THEN CALL Warn(SolverName,& "Both nodes are marked Dont Move, stopping shifting.") EXIT END IF !Work out which to shift !If one is marked DontMove (corner node, or 2nd inland and weird shape) ! then move the other !Otherwise, shift whichever is further forward IF(DontMove(i)) THEN ShiftSecond = .FALSE. DontMove(i-1) = .TRUE. IF(Debug) PRINT *,ParEnv % MyPE,'Debug, do not move second: ',i ELSE IF(DontMove(i-1)) THEN ShiftSecond = .TRUE. DontMove(i) = .TRUE. IF(Debug) PRINT *,ParEnv % MyPE,'Debug, do not move first: ',i-1 ELSE IF(SUM(Rot_z_coords(i,:)) > SUM(Rot_z_coords(i-1,:))) THEN ShiftSecond = .TRUE. ELSE ShiftSecond = .FALSE. END IF END IF IF(ShiftSecond) THEN ShiftIdx = i ShiftToIdx = i-1 ELSE ShiftIdx = i-1 ShiftToIdx = i END IF !Calculate ShiftTo, depends on ShiftSecond !Also update shifted Rot_y_coords IF(ShiftSecond) THEN !If this node has already been moved the other way, leave it alone IF((UpdatedDirection(ShiftIdx) < 0) .AND. .NOT. DontMove(ShiftToIdx)) CYCLE ShiftTo = Rot_y_coords(i-1,2) + epsShift IF(Rot_y_coords(i,2) < ShiftTo) UpdatedDirection(ShiftIdx) = 1 Rot_y_coords(i,1) = MAX(Rot_y_coords(i,1), ShiftTo) Rot_y_coords(i,2) = MAX(Rot_y_coords(i,2), ShiftTo) ELSE !If this node has already been moved the other way, leave it alone IF((UpdatedDirection(ShiftIdx) > 0) .AND. .NOT. DontMove(ShiftToIdx)) CYCLE ShiftTo = Rot_y_coords(i,1) - epsShift IF(Rot_y_coords(i,1) > ShiftTo) UpdatedDirection(ShiftIdx) = -1 Rot_y_coords(i-1,1) = MIN(Rot_y_coords(i-1,1), ShiftTo) Rot_y_coords(i-1,2) = MIN(Rot_y_coords(i-1,2), ShiftTo) END IF IF(Debug) PRINT *,ParEnv % MyPE, ' Debug, updating col ',ShiftIdx,& ' rot_y_coords: ',Rot_y_coords(ShiftIdx,:) UpdatedColumn(ShiftIdx) = .TRUE. !Shift all the nodes in this column DO j=1,Mesh % NumberOfNodes IF(FNColumns(j) /= FrontColumnList(ShiftIdx)) CYCLE NodeHolder(1) = Mesh % Nodes % x(j) + Advance((Perm(j)-1)*DOFs + 1) NodeHolder(2) = Mesh % Nodes % y(j) + Advance((Perm(j)-1)*DOFs + 2) NodeHolder(3) = Mesh % Nodes % z(j) + Advance((Perm(j)-1)*DOFs + 3) NodeHolder = MATMUL(RotationMatrix, NodeHolder) Displace = 0.0_dp Displace(2) = ShiftTo - NodeHolder(2) !Check if node already projectable IF(ShiftSecond) THEN IF(Displace(2) < 0 ) CYCLE ELSE IF(Displace(2) > 0) CYCLE END IF Displace = MATMUL(TRANSPOSE(RotationMatrix), Displace) IF(Debug) PRINT *,ParEnv % MyPE, 'Debug, shifting node ',j,' col ',ShiftIdx,'xyz: ',& Mesh % Nodes % x(j)+Advance((Perm(j)-1)*DOFs + 1),& Mesh % Nodes % y(j)+Advance((Perm(j)-1)*DOFs + 2),& Mesh % Nodes % z(j)+Advance((Perm(j)-1)*DOFs + 3),& ' by ',Displace,' to ensure projectability. 1' !Adjust the variable values to shift nodes in line. Advance((Perm(j)-1)*DOFs + 1) = Advance((Perm(j)-1)*DOFs + 1) + Displace(1) Advance((Perm(j)-1)*DOFs + 2) = Advance((Perm(j)-1)*DOFs + 2) + Displace(2) END DO END IF END DO IF(ANY(UpdatedColumn)) MovedOne = .TRUE. IF(Debug) PRINT *,ParEnv % MyPE, ' Debug, MovedOne: ',MovedOne IF(Debug) PRINT *,ParEnv % MyPE, ' Debug, UpdatedColumn: ',UpdatedColumn END DO !----------------------------------------------- ! Check for pinnacles or rifts which can't be made ! projectable. Set the columns in a line and mark them ! for removal by Remesh.F90 !----------------------------------------------- !TODO: At least 3 separate issues here: ! ! *This one should be fixed now, at the end we check for groups disappearing* ! 1 : If there's a tangled group, then another, superset tangled ! group, this causes "programming error, tangled group has zero nodes" ! potential solution here is to just remove this error message, because ! its not really a problem *I DON'T THINK*. However, if so, need to ! check for it, and shift TangledGroups down one, so Remesh behaves ! ! *This one should now be fixed, because we do all the unprojectable shifting first:* ! 2 : Here we iteratively: check tangled, then check unprojectable ! If correcting unprojectability causes further tangling, this ! isn't caught. Potential solutions: ! Improve parallel unprojectability checker (dirty fix, won't catch every poss case) ! Better: recheck previously tangled regions ! ! 3 : Long straight sides which also happen to tangle with the end of their headlands ! are not well handled. All the nodes along these sides are shifted out the way ! in the new mesh NoTangledGroups = 0 DO i=1,FrontLineCount-2 IF(Tangled(i)) CYCLE j = i+2 !If the column two columns away is less than 2*epsShift away, and there !is a change of direction, problems... IF((Rot_y_coords(j,1) - Rot_y_coords(i,2)) < 2*epsShift) THEN IF(((Rot_z_coords(i,1) < Rot_z_coords(i+1,1)) .NEQV. & (Rot_z_coords(i+1,1) < Rot_z_coords(j,1))) .OR. & ((Rot_z_coords(i,2) < Rot_z_coords(i+1,2)) .NEQV. & (Rot_z_coords(i+1,2) < Rot_z_coords(j,2)))) THEN !Either a pinnacle or a slit (i.e protrusion or rift) Protrusion = SUM(Rot_z_coords(i+1,:)) > SUM(Rot_z_coords(i,:)) IF(Debug) PRINT *,'Debug, protrusion: ',Protrusion NoTangledGroups = NoTangledGroups + 1 PivotIdx = i+1 DO k=2,PivotIdx p1(1) = Rot_y_coords(k,2) p1(2) = Rot_z_coords(k,2) p2(1) = Rot_y_coords(k-1,2) p2(2) = Rot_z_coords(k-1,2) DO m=FrontLineCount-1,PivotIdx,-1 IF(k==m) CYCLE !first two will always intersect by definition, not what we want q1(1) = Rot_y_coords(m,1) q1(2) = Rot_z_coords(m,1) q2(1) = Rot_y_coords(m+1,1) q2(2) = Rot_z_coords(m+1,1) CALL LineSegmentsIntersect ( p1, p2, q1, q2, intersect, intersect_flag ) IF(intersect_flag) EXIT END DO IF(intersect_flag) EXIT END DO IF(intersect_flag) THEN !Found an intersection point (intersect) PRINT *,'Debug, found tangle intersection ',intersect,' leaving last tangled nodes: ',k,m Tangled(k:m) = .TRUE. TangledPivotIdx(k:m) = PivotIdx TangledShiftTo(k:m) = intersect(1) IF(Protrusion) THEN TangledGroup(k:m) = NoTangledGroups ELSE TangledGroup(k:m) = -NoTangledGroups END IF ELSE PRINT *,'Debug, found no intersection, so nodes are not QUITE tangled: ',i,j Tangled(i:j) = .TRUE. TangledPivotIdx(i:j) = PivotIdx TangledShiftTo(i:j) = (SUM(rot_y_coords(i,:)) + SUM(rot_y_coords(j,:))) / 4.0_dp IF(Protrusion) THEN TangledGroup(i:j) = NoTangledGroups ELSE TangledGroup(i:j) = -NoTangledGroups END IF END IF END IF END IF END DO !Check for a tangled group 'swallowing' another county = 0 DO i=1,NoTangledGroups IF(COUNT(TangledGroup == i) > 0) CYCLE county = county + 1 !Shift group numbers down one DO j=1,SIZE(TangledGroup) IF(TangledGroup(j) > i) TangledGroup(j) = TangledGroup(j) - 1 END DO END DO NoTangledGroups = NoTangledGroups - county !Strategy: ! Shift all nodes to near (offset) the y-coordinate ! of the tangle pivot (i+1, above) DO i=1,NoTangledGroups n = COUNT(TangledGroup == i) IF(n==0) THEN n = COUNT(TangledGroup == -i) Protrusion = .FALSE. ELSE Protrusion = .TRUE. END IF IF(n==0) CALL Fatal(SolverName, "Programming error: tangled group has 0 nodes?") !Get the pivot node index, ShiftToY, and tangled node range FirstTangleIdx = 0 LastTangleIdx = 0 DO j=1,FrontLineCount IF(TangledGroup(j) == i .OR. TangledGroup(j) == -i) THEN IF(FirstTangleIdx == 0) FirstTangleIdx = j LastTangleIdx = j PivotIdx = TangledPivotIdx(j) ShiftToY = TangledShiftTo(j) !EXIT END IF END DO IF(LastTangleIdx - FirstTangleIdx /= n-1) CALL Fatal(SolverName, & "Programming error: wrong number of nodes in TangledGroup") IF(Debug) PRINT *,'Debug, tangled pivot index, ShiftToY: ', PivotIdx, ShiftToY LeftY = ShiftToY - (epsTangle * ((n-1) / 2.0_dp)) RightY = ShiftToY + (epsTangle * ((n-1) / 2.0_dp)) !Potentially 'squeezed' on either side by neighbour columns SqueezeLeft = .FALSE. SqueezeRight = .FALSE. IF(LeftY < Rot_y_coords(FirstTangleIdx-1,2) + epsTangle) THEN SqueezeLeft = .TRUE. IF( ( (Rot_y_coords(LastTangleIdx+1,1) - epsTangle) - (epsTangle * (n-1)) ) < & (Rot_y_coords(FirstTangleIdx-1,2) + epsTangle)) THEN SqueezeRight = .TRUE. END IF END IF IF(RightY > (Rot_y_coords(LastTangleIdx+1,1) - epsTangle)) THEN SqueezeRight = .TRUE. IF( ( (Rot_y_coords(LastTangleIdx+1,1) - epsTangle) - (epsTangle * (n-1)) ) < & (Rot_y_coords(FirstTangleIdx-1,2) + epsTangle)) THEN SqueezeLeft = .TRUE. END IF END IF IF(Debug) PRINT *,'Debug, LeftY: ',LeftY,' RightY: ',RightY,' SqueezeL: ',& SqueezeLeft,' SqueezeR: ',SqueezeRight,' prev: ',Rot_y_coords(FirstTangleIdx-1,2),& ' next: ',Rot_y_coords(LastTangleIdx+1,1) !If squeezed, adjust the new y coord of the tangled nodes IF(SqueezeLeft .AND. SqueezeRight) THEN thisEps = (Rot_y_coords(LastTangleIdx+1,1) - Rot_y_coords(FirstTangleIdx-1,2)) / (n+1) LeftY = Rot_y_coords(FirstTangleIdx-1,2) + thisEps RightY = Rot_y_coords(LastTangleIdx+1,1) - thisEps ShiftToY = (RightY + LeftY) / 2.0_dp !midpoint ELSE IF(SqueezeLeft) THEN Shift = (Rot_y_coords(FirstTangleIdx-1,2) + epsTangle) - LeftY LeftY = LeftY + Shift RightY = RightY + Shift ShiftToY = ShiftToY + Shift ELSE IF(SqueezeRight) THEN Shift = (Rot_y_coords(LastTangleIdx+1,1) - epsTangle) - RightY LeftY = LeftY + Shift RightY = RightY + Shift ShiftToY = ShiftToY + Shift END IF !Where tangling occurs, we shift the tangled columns to be !1m apart in a series. Then Remesh gets rid of them DO j=FirstTangleIdx,LastTangleIdx IF(.NOT. (TangledGroup(j) == i .OR. TangledGroup(j) == -i)) & CALL Fatal(SolverName, "Programming error: node in specified idx range not tangled?") thisY = LeftY + ((REAL(j - FirstTangleIdx)/(n-1)) * (RightY - LeftY)) IF(Debug) PRINT *,'Debug, thisY: ',thisY, j DO k=1,Mesh % NumberOfNodes IF(FNColumns(k) /= FrontColumnList(j)) CYCLE NodeHolder(1) = Mesh % Nodes % x(k) + Advance((Perm(k)-1)*DOFs + 1) NodeHolder(2) = Mesh % Nodes % y(k) + Advance((Perm(k)-1)*DOFs + 2) NodeHolder(3) = Mesh % Nodes % z(k) + Advance((Perm(k)-1)*DOFs + 3) IF(Debug) PRINT *, ParEnv % MyPE, ' Debug, pre shift tangled node ',k,': ',NodeHolder NodeHolder = MATMUL(RotationMatrix, NodeHolder) Displace = 0.0_dp Displace(2) = thisY - NodeHolder(2) Displace = MATMUL(TRANSPOSE(RotationMatrix), Displace) IF(Debug) PRINT *,ParEnv % MyPE, ' Debug, shifting node ',k, ' col ',j,& ' by ',Displace,' to detangle.' !Adjust the variable values to shift nodes in line. Advance((Perm(k)-1)*DOFs + 1) = Advance((Perm(k)-1)*DOFs + 1) + Displace(1) Advance((Perm(k)-1)*DOFs + 2) = Advance((Perm(k)-1)*DOFs + 2) + Displace(2) !Set this variable so that Remesh knows TangledVar % Values(TangledVar % Perm(k)) = 1.0_dp * ABS(TangledGroup(j)) END DO END DO END DO !--------------------------------------- !Done, just deallocations FirstTime = .FALSE. DEALLOCATE(FrontPerm, & TopPerm, & FrontNodeNums, & Rot_y_coords, & UpdatedColumn, & UpdatedDirection, & Tangled, & TangledGroup, & TangledPivotIdx, & TangledShiftTo, & FrontColumnList, & FrontLocalNodeNumbers) END SUBROUTINE FrontAdvance3D