diff --git a/elmerice/Solvers/Calving3D_lset.F90 b/elmerice/Solvers/Calving3D_lset.F90 index fe1fa91a74..5090bc6223 100644 --- a/elmerice/Solvers/Calving3D_lset.F90 +++ b/elmerice/Solvers/Calving3D_lset.F90 @@ -78,6 +78,7 @@ SUBROUTINE Find_Calving3D_LSet ( Model, Solver, dt, TransientSimulation ) angle,angle0,a1(2),a2(2),b1(2),b2(2),a2a1(2),isect(2),front_extent(4), & buffer, gridmesh_dx, FrontLeft(2), FrontRight(2), ElemEdge(2,5), & CalvingLimit, CrevPenetration, PrevValue, PartMinDist, & + Calv_delta_t, Calv_dmax, Calv_k, RandomNumber, x, maxprob, Mu, Prob, calv_f, & #ifdef USE_ISO_C_BINDINGS rt0, rt #else @@ -90,10 +91,10 @@ SUBROUTINE Find_Calving3D_LSet ( Model, Solver, dt, TransientSimulation ) EdgeX(:), EdgeY(:), EdgePoly(:,:), CrevOrient(:,:) CHARACTER(LEN=MAX_NAME_LEN) :: SolverName, DistVarname, & FrontMaskName,TopMaskName,BotMaskName,LeftMaskName,RightMaskName,InflowMaskName, & - PC_EqName, Iso_EqName, VTUSolverName, EqName + PC_EqName, Iso_EqName, VTUSolverName, EqName, FileName LOGICAL :: Found, Parallel, Boss, Debug, FirstTime = .TRUE., CalvingOccurs=.FALSE., & - SaveParallelActive, LeftToRight, inside, Complete,& - LatCalvMargins, FullThickness, UnfoundConstraint + SaveParallelActive, LeftToRight, inside, Complete, SaveToFile, & + LatCalvMargins, FullThickness, UnfoundConstraint, RandomCalving, FileCreated LOGICAL, ALLOCATABLE :: RemoveNode(:), IMNOnFront(:), IMOnMargin(:), & IMNOnLeft(:), IMNOnRight(:), IMElmONFront(:), IMElmOnLeft(:), IMElmOnRight(:), FoundNode(:), & IMElemOnMargin(:), DeleteMe(:), IsCalvingNode(:), PlaneEdgeElem(:), EdgeNode(:), UsedElem(:), & @@ -104,7 +105,7 @@ SUBROUTINE Find_Calving3D_LSet ( Model, Solver, dt, TransientSimulation ) SAVE :: FirstTime, SolverName, Params, Parallel, Boss, dim, Debug, & DistVarName, PC_EqName, Iso_EqName, LeftConstraint, & RightConstraint, FrontConstraint,TopMaskName, BotMaskName, & - LeftMaskName, RightMaskName, FrontMaskName, InflowMaskName + LeftMaskName, RightMaskName, FrontMaskName, InflowMaskName, FileCreated !---------------Get Variables and Parameters for Solution----------- @@ -146,6 +147,82 @@ SUBROUTINE Find_Calving3D_LSet ( Model, Solver, dt, TransientSimulation ) IF(Debug) PRINT *,'Front, Left, Right constraints: ',FrontConstraint,LeftConstraint,RightConstraint END IF !FirstTime + RandomCalving = ListGetLogical(Params,"Random Calving") + IF(RandomCalving) THEN + IF(Boss) THEN + + Calv_k = ListGetConstReal(Params, "Calving k",Found) + IF(.NOT. Found) CALL FATAL(SolverName, "Random calving specified but no 'Calving k' value") + Calv_delta_t = ListGetConstReal(Params, "Calving delta t",Found) + IF(.NOT. Found) CALL FATAL(SolverName, "Random calving specified but no 'Calving delta t' value") + Calv_dmax = ListGetConstReal(Params, "Calving dmax",Found) + IF(.NOT. Found) CALL FATAL(SolverName, "Random calving specified but no 'Calving dmax' value") + Calv_f = ListGetConstReal(Params, "Calving f",Found) + IF(.NOT. Found) CALL FATAL(SolverName, "Random calving specified but no 'Calving f' value") + + SaveToFile = ListGetLogical(Params,"Save Probability to File") + + ! calculate dcrev from calving probability + !x = (dcrev - dmin) / (1 - dmin) + !x[dcrev < dmin] = 0 + !mu = (x**k) / ( (x**k) + (1-x)**k ) + !prob = 1 - np.exp(-mu * delta_t) + + maxprob = 1.0_dp - EXP(-1.0_dp * Calv_delta_t * Calv_f) + + ! one random number per timestep + CALL Random_Seed() + CALL Random_Number(RandomNumber) + + IF(RandomNumber >= maxprob + AEPS) THEN + CALL INFO(SolverName, 'shortening random number') + Prob = maxprob + ELSE + Prob = RandomNumber + END IF + + Mu = -LOG(1.0_dp-Prob)/(Calv_delta_t*calv_f) + + CrevPenetration = Calv_dmax - (LOG(1/mu) / Calv_k) + + IF(SaveToFile) THEN + Filename = ListGetString(Params,"Probability File Name", Found) + IF(.NOT. Found) THEN + CALL WARN(SolverName, 'Output file name not given so using Probability.txt') + Filename = "Probability.txt" + END IF + + ! write to file + IF(FileCreated) THEN + OPEN( 47, FILE=filename, STATUS='UNKNOWN', ACCESS='APPEND') + ELSE + OPEN( 47, FILE=filename, STATUS='UNKNOWN') + WRITE(47, '(A)') "Calving Probability Output File" + WRITE(47, '(A)') "TimeStep, Time, RandomNumber, Probability, Mu, DCrev" + END IF + + WRITE(47, *) GetTimestep(), GetTime(), RandomNumber, Prob, Mu, CrevPenetration + CLOSE(47) + FileCreated = .TRUE. + END IF + ELSE + CrevPenetration = 0.0_dp + END IF + + CALL MPI_ALLREDUCE(MPI_IN_PLACE, CrevPenetration, 1, MPI_DOUBLE_PRECISION, MPI_MAX, ELMER_COMM_WORLD, ierr) + + PRINT*, 'CrevPenetration: ', CrevPenetration + + WRITE (Message,'(A,F1.2,A)') 'Calving occuring using crevasse penetration:', CrevPenetration, & + 'due to probability.' + CALL INFO(SolverName, Message) + ELSE + CrevPenetration = ListGetConstReal(Params, "Crevasse Penetration",Found, DefValue = 1.0_dp) + IF(.NOT. Found) CALL Info(SolverName, "No Crevasse Penetration specified so assuming full thickness") + PRINT*, 'CrevPenetration: ', CrevPenetration + IF(CrevPenetration > 1 .OR. CrevPenetration < 0) CALL FATAL(SolverName, "Crevasse Penetraion must be between 0-1") + END IF + CalvingOccurs = .FALSE. Mesh => Model % Mesh @@ -1495,6 +1572,7 @@ SUBROUTINE Find_Calving3D_LSet ( Model, Solver, dt, TransientSimulation ) ! EXIT ! END IF !END DO + CalvingValues(CalvingPerm) = DistValues(DistPerm) + 1 CALL WARN(SolverName, 'No crevasses so not calculating signed distance') !RETURN ELSE diff --git a/elmerice/Solvers/CalvingGeometry.F90 b/elmerice/Solvers/CalvingGeometry.F90 index 1eb5c49b7e..3ea3e308b7 100644 --- a/elmerice/Solvers/CalvingGeometry.F90 +++ b/elmerice/Solvers/CalvingGeometry.F90 @@ -1451,15 +1451,15 @@ FUNCTION PointInPolygon2D(Polygon, Point, buffer) RESULT(inside) windingnumber=100 DO i=1, n-1 ! polygon y i <= point y - IF(ZPolygon(2,i) <= ZPoint(2)) THEN !start with y<=P.y - IF(ZPolygon(2, i+1) > ZPoint(2)) THEN !upward crossing + IF(ZPolygon(2,i) <= ZPoint(2) + buf) THEN !start with y<=P.y + IF(ZPolygon(2, i+1) > ZPoint(2) - buf) THEN !upward crossing left=IsLeft(ZPolygon(:, i), ZPolygon(:, i+1), ZPoint(:)) IF(left > buf) THEN !p is to left of intersect windingnumber=windingnumber+1 !valid up intersect END IF END IF ELSE !start at y> point y - IF(ZPolygon(2, i+1) <= ZPoint(2)) THEN ! downward crossing + IF(ZPolygon(2, i+1) <= ZPoint(2) + buf) THEN ! downward crossing Left = IsLeft(ZPolygon(:, i), ZPolygon(:, i+1), ZPoint(:)) IF(left < buf) THEN ! p right of edge windingnumber=windingnumber-1 @@ -2830,7 +2830,7 @@ SUBROUTINE GetDomainEdge(Model, Mesh, TopPerm, OrderedNodes, OrderedNodeNums, Pa !------------ DEALLOCATIONS ------------------ - DEALLOCATE(OnEdge, UnorderedNodeNums, GlobalCorners, CornerParts, PCornerCounts) + DEALLOCATE(OnEdge, UnorderedNodeNums, GlobalCorners, CornerParts, PCornerCounts, OrderedNodeNums) IF(Boss .AND. Parallel) THEN !Deallocations DEALLOCATE(UnorderedNodes % x, & @@ -2843,6 +2843,8 @@ SUBROUTINE GetDomainEdge(Model, Mesh, TopPerm, OrderedNodes, OrderedNodeNums, Pa OrderedGlobalNodeNums) END IF + IF(.NOT. Boss) DEALLOCATE(UnorderedNodes % x, UnorderedNodes % y, UnorderedNodes % z) + END SUBROUTINE GetDomainEdge ! Copies over time variables and creates coordinate vars. Basically pinched @@ -3489,14 +3491,15 @@ SUBROUTINE SwitchMesh(Model, Solver, OldMesh, NewMesh) TYPE(Matrix_t), POINTER :: WorkMatrix=>NULL() LOGICAL :: Found, Global, GlobalBubbles, Debug, DoPrevValues, & NoMatrix, DoOptimizeBandwidth, PrimaryVar, HasValuesInPartition, & - PrimarySolver + PrimarySolver,CreatedParMatrix LOGICAL, POINTER :: UnfoundNodes(:)=>NULL(), BulkUnfoundNodes(:)=>NULL() - INTEGER :: i,j,k,DOFs, nrows,n, dummyint + INTEGER :: i,j,k,DOFs, nrows,n, dummyint, ierr INTEGER, POINTER :: WorkPerm(:)=>NULL(), SolversToIgnore(:)=>NULL(), & SurfaceMaskPerm(:)=>NULL(), BottomMaskPerm(:)=>NULL() REAL(KIND=dp), POINTER :: WorkReal(:)=>NULL(), WorkReal2(:)=>NULL(), PArray(:,:) => NULL() REAL(KIND=dp) :: FrontOrientation(3), RotationMatrix(3,3), UnRotationMatrix(3,3), & globaleps, localeps + LOGICAL, ALLOCATABLE :: PartActive(:) CHARACTER(LEN=MAX_NAME_LEN) :: SolverName, WorkName INTERFACE @@ -3555,6 +3558,10 @@ END SUBROUTINE InterpolateMeshToMesh !---------------------------------------------- Var => OldMesh % Variables + + ALLOCATE(PartActive(ParEnv % PEs)) + CreatedParMatrix = .FALSE. + DO WHILE( ASSOCIATED(Var) ) DoPrevValues = ASSOCIATED(Var % PrevValues) @@ -3633,6 +3640,9 @@ END SUBROUTINE InterpolateMeshToMesh END IF END IF + IF(.NOT. CreatedParMatrix) & + CALL MPI_AllGather(.NOT. NoMatrix, 1, MPI_LOGICAL, PartActive, 1, MPI_LOGICAL, ELMER_COMM_WORLD, ierr) + IF ( ASSOCIATED(Var % EigenValues) ) THEN n = SIZE(Var % EigenValues) @@ -3676,6 +3686,20 @@ END SUBROUTINE InterpolateMeshToMesh IF(ASSOCIATED(WorkSolver % Matrix)) CALL FreeMatrix(WorkSolver % Matrix) WorkSolver % Matrix => WorkMatrix + ! bit of a hack + ! since ParEnv become a pointer to ParMatrix we need to ensure one ParMatrix is formed + ! it needs to be from a solver present on all parts hence the all gather further up. + ! it seems we only need to this once per timestep/interpolation as ParEnv will have some thing + ! to point to. If we don't do this ParEnv % PEs, % MyPE etc. all become nans mucking eveything up! + IF ( ASSOCIATED(WorkSolver % Matrix) .and. ALL(PartActive) .and. .NOT. CreatedParMatrix) THEN + IF (.NOT. ASSOCIATED(WorkSolver % Matrix % ParMatrix) ) THEN + WorkSolver % Mesh => NewMesh + + CALL ParallelInitMatrix( WorkSolver, WorkSolver % Matrix, WorkPerm) + CreatedParMatrix = .TRUE. + END IF + END IF + NULLIFY(WorkMatrix) !NOTE: We don't switch Solver % Variable here, because @@ -4648,27 +4672,28 @@ SUBROUTINE GetCalvingEdgeNodes(Mesh, Parallel, Shared, TotalCount) END SUBROUTINE GetCalvingEdgeNodes - SUBROUTINE MeshVolume(Mesh, Parallel, Volume, ElemMask) + SUBROUTINE MeshVolume(Mesh, Parallel, Volume, ElemMask, Centroid) TYPE(Mesh_t), POINTER :: Mesh LOGICAL :: Parallel REAL(kind=dp) :: Volume LOGICAL, OPTIONAL :: ElemMask(:) + REAL(kind=dp), OPTIONAL :: Centroid(3) !----------------------------- TYPE(Element_t), POINTER :: Element INTEGER :: i, j, NBdry, NBulk, n, ierr INTEGER, ALLOCATABLE :: ElementNodes(:) REAL(kind=dp), ALLOCATABLE :: Vertices(:,:), Vectors(:,:), PartVolume(:) - REAL(kind=dp) :: det, det1, det2, det3 + REAL(kind=dp) :: det, Centre(3) NBdry = Mesh % NumberOfBoundaryElements NBulk = Mesh % NumberOfBulkElements ALLOCATE(Vertices(4,3), Vectors(3,3)) - ! calculate volume of each bulk tetra. Add these together to get mesh volume Volume = 0.0_dp + IF(PRESENT(Centroid)) Centroid = 0.0_dp DO, i=1, NBulk IF(PRESENT(ElemMask)) THEN IF(.NOT. ElemMask(i)) CYCLE @@ -4697,11 +4722,17 @@ SUBROUTINE MeshVolume(Mesh, Parallel, Volume, ElemMask) - Vectors(1,2) * (Vectors(2,1)*Vectors(3,3) - Vectors(2,3)*Vectors(3,1)) & + Vectors(1,3) * (Vectors(2,1)*Vectors(3,2) - Vectors(2,2)*Vectors(3,1))) + Centre(1) = SUM(Vertices(:,1))/4 + Centre(2) = SUM(Vertices(:,2))/4 + Centre(3) = SUM(Vertices(:,3))/4 + ! tetra volume = det/6 Volume = Volume + Det/6 - + IF(PRESENT(Centroid)) Centroid = Centroid + Det/6 * Centre END DO + IF(PRESENT(Centroid)) Centroid = Centroid / Volume + ! if parallel calculate total mesh volume over all parts IF(Parallel) THEN ALLOCATE(PartVolume(ParEnv % PEs)) @@ -5319,7 +5350,7 @@ SUBROUTINE InterpolateUnfoundSharedPoint3D( NodeNumber, Mesh, Variables, Unfound ! all parallel communication changed to use NoUsedNeighbours so neighbouring procs ! of those with zero suppnodes (no info) do not over allocate (eg allocate nans) !share SuppNodeMask - ALLOCATE(PartSuppNodeMask(NoUsedNeighbours+1, 25, MaskCount)) + ALLOCATE(PartSuppNodeMask(NoUsedNeighbours+1, 50, MaskCount)) PartSuppNodeMask = .FALSE. PartSuppNodeMask(1,:NoSuppNodes,:) = SuppNodeMask counter=0 @@ -5338,7 +5369,7 @@ SUBROUTINE InterpolateUnfoundSharedPoint3D( NodeNumber, Mesh, Variables, Unfound END DO !share SuppNodePMask for prevvalues - ALLOCATE(PartSuppNodePMask(NoUsedNeighbours+1, 25, PMaskCount)) + ALLOCATE(PartSuppNodePMask(NoUsedNeighbours+1, 50, PMaskCount)) PartSuppNodePMask = .FALSE. PartSuppNodePMask(1,:NoSuppNodes,:) = SuppNodePMask counter=0 @@ -6217,6 +6248,8 @@ SUBROUTINE GetFrontCorners(Model, Solver, FrontLeft, FrontRight) NULLIFY(SidePerm) END DO + DEALLOCATE(FrontPerm, TopPerm, LeftPerm, RightPerm) + END SUBROUTINE GetFrontCorners SUBROUTINE ValidateNPCrevassePaths(Mesh, CrevassePaths, OnLeft, OnRight, FrontLeft, FrontRight, & @@ -6233,7 +6266,8 @@ SUBROUTINE ValidateNPCrevassePaths(Mesh, CrevassePaths, OnLeft, OnRight, FrontLe REAL(KIND=dp) :: RotationMatrix(3,3), UnRotationMatrix(3,3), FrontDist, MaxDist, & ShiftTo, Dir1(2), Dir2(2), CCW_value,a1(2),a2(2),b1(2),b2(2),intersect(2), & StartX, StartY, EndX, EndY, Orientation(3), temp, NodeHolder(3), err_buffer,& - yy, zz, gradient, c, intersect_z, SideCorner(3), MinDist, TempDist, IsBelowMean + yy, zz, gradient, c, intersect_z, SideCorner(3), MinDist, TempDist, IsBelowMean,& + PolyMin, PolyMax REAL(KIND=dp), ALLOCATABLE :: ConstrictDirection(:,:), REdge(:,:), Polygons(:,:) REAL(KIND=dp), POINTER :: WorkReal(:) TYPE(CrevassePath_t), POINTER :: CurrentPath, OtherPath, WorkPath, LeftPath, RightPath @@ -6255,6 +6289,13 @@ SUBROUTINE ValidateNPCrevassePaths(Mesh, CrevassePaths, OnLeft, OnRight, FrontLe Debug = .FALSE. Snakey = .TRUE. + IF(PRESENT(GridSize)) THEN + err_buffer = GridSize/1000.0 + ELSE + err_buffer = AEPS + END IF + IF( err_buffer < AEPS) err_buffer = AEPS + ! if on lateral margin need to make sure that glacier corner is within crev. ! if it lies outside the crev then the crev isn't really on front but on the lateral corner ! first and last both actually on same lateral margin @@ -6295,6 +6336,8 @@ SUBROUTINE ValidateNPCrevassePaths(Mesh, CrevassePaths, OnLeft, OnRight, FrontLe CurrentPath => WorkPath END DO + DEALLOCATE(Polygons, PolyStart, PolyEnd) + CALL GetCalvingPolygons(Mesh, CrevassePaths, EdgeX, EdgeY, Polygons, PolyStart, PolyEnd, GridSize) ! invalid lateral crevs must first be removed before this subroutine CurrentPath => CrevassePaths path=0 @@ -6338,7 +6381,7 @@ SUBROUTINE ValidateNPCrevassePaths(Mesh, CrevassePaths, OnLeft, OnRight, FrontLe IF(Onside /= 0 .AND. LatCalvMargins) AddLateralMargins = .TRUE. orientation(3) = 0.0_dp - IF( ABS(StartX-EndX) < AEPS) THEN + IF( ABS(StartX-EndX) < err_buffer) THEN ! front orientation is aligned with y-axis Orientation(2) = 0.0_dp IF(EndY > StartY) THEN @@ -6346,7 +6389,7 @@ SUBROUTINE ValidateNPCrevassePaths(Mesh, CrevassePaths, OnLeft, OnRight, FrontLe ELSE Orientation(1)=-1.0_dp END IF - ELSE IF (ABS(StartY-EndY) StartX) THEN @@ -6358,11 +6401,15 @@ SUBROUTINE ValidateNPCrevassePaths(Mesh, CrevassePaths, OnLeft, OnRight, FrontLe CALL ComputePathExtent(CrevassePaths, Mesh % Nodes, .TRUE.) ! endx always greater than startx ! check if yextent min smaller than starty - IF(CurrentPath % Right == StartY .OR. & - CurrentPath % Right == EndY) THEN - Orientation(2)=1.0_dp - ELSE + + PolyMin = MINVAL(Polygons(2,PolyStart(path):PolyEnd(path))) + PolyMax = MAXVAL(Polygons(2,PolyStart(path):PolyEnd(path))) + + IF(ABS(CurrentPath % Right - PolyMax) > & + CurrentPath % Left - PolyMin) THEN Orientation(2)=-1.0_dp + ELSE + Orientation(2)=1.0_dp END IF Orientation(1)=Orientation(2)*(EndY-StartY)/(StartX-EndX) END IF @@ -6396,12 +6443,6 @@ SUBROUTINE ValidateNPCrevassePaths(Mesh, CrevassePaths, OnLeft, OnRight, FrontLe REdge(3,i) = NodeHolder(3) END DO - IF(PRESENT(GridSize)) THEN - err_buffer = GridSize/10 - ELSE - err_buffer = 0.0_dp - END IF - ! crop edge around crev ends crop=0 DO i=1, EdgeLength @@ -7615,9 +7656,6 @@ SUBROUTINE EnforceLateralMargins(Model, SolverParams) xx = Mesh % Nodes % x(i) yy = Mesh % Nodes % y(i) - inside = PointInPolygon2D(RailPoly, (/xx,yy/)) - IF(inside) CYCLE - IF(LeftPerm(i) > 0) THEN ! check if on left side mindist = HUGE(1.0_dp) DO j=1, Nl-1 @@ -7653,6 +7691,9 @@ SUBROUTINE EnforceLateralMargins(Model, SolverParams) END IF IF(FrontPerm(i) > 0) THEN ! check if front is on rail eg advance on narrowing rails + inside = PointInPolygon2D(RailPoly, (/xx,yy/)) + IF(inside) CYCLE + mindist = HUGE(1.0_dp) DO j=1, Nr-1 tempdist = PointLineSegmDist2D((/xR(j), yR(j)/),(/xR(j+1), yR(j+1)/), (/xx, yy/)) @@ -7851,8 +7892,8 @@ SUBROUTINE CalvingStatsMMG(Params, Mesh, Mask, ElemMask, FileCreated, MaxBergVol LOGICAL, ALLOCATABLE :: FoundNode(:), UsedElem(:), IcebergElem(:), GotNode(:), & NodeCount(:) CHARACTER(LEN=MAX_NAME_LEN) :: Filename - REAL(kind=dp), ALLOCATABLE :: BergVolumes(:), BergExtents(:) - REAL(kind=dp) :: BergVolume, extent(4) + REAL(kind=dp), ALLOCATABLE :: BergVolumes(:), BergExtents(:), BergCentroids(:) + REAL(kind=dp) :: BergVolume, extent(4), Centroid(3) Filename = ListGetString(Params,"Calving Stats File Name", Found) IF(.NOT. Found) THEN @@ -7867,7 +7908,7 @@ SUBROUTINE CalvingStatsMMG(Params, Mesh, Mask, ElemMask, FileCreated, MaxBergVol !limit here of 10 possible mesh 'islands' ALLOCATE(FoundNode(NNodes), NodeCount(NNodes), ElNodes(4), & UsedElem(NBulk), IceBergElem(NBulk), BergVolumes(100), & - BergExtents(100 * 4)) + BergExtents(100 * 4), BergCentroids(100*3)) FoundNode = .FALSE. NodeCount = .NOT. Mask UsedElem = .FALSE. !count of elems used @@ -7907,7 +7948,7 @@ SUBROUTINE CalvingStatsMMG(Params, Mesh, Mask, ElemMask, FileCreated, MaxBergVol IcebergElem(i) = .TRUE. END DO iceberg = iceberg + 1 - CALL MeshVolume(Mesh, .FALSE., BergVolume, IcebergElem) + CALL MeshVolume(Mesh, .FALSE., BergVolume, IcebergElem, Centroid) CALL IcebergExtent(Mesh, IcebergElem, Extent) IF(SIZE(BergVolumes) < Iceberg) CALL DoubleDPVectorSize(BergVolumes) @@ -7916,8 +7957,11 @@ SUBROUTINE CalvingStatsMMG(Params, Mesh, Mask, ElemMask, FileCreated, MaxBergVol IF(SIZE(BergExtents) < Iceberg*4) CALL DoubleDPVectorSize(BergExtents) BergExtents(iceberg*4-3:iceberg*4) = Extent + IF(SIZE(BergCentroids) < Iceberg*3) CALL DoubleDPVectorSize(BergCentroids) + BergCentroids(iceberg*3-2:iceberg*3) = Centroid + IF(Iceberg > 0) THEN ! not first time - PRINT*, 'Iceberg no.', Iceberg, BergVolume, 'extent', extent + PRINT*, 'Iceberg no.', Iceberg, BergVolume, 'extent', extent, 'centroid', centroid END IF END IF END DO @@ -7940,8 +7984,9 @@ SUBROUTINE CalvingStatsMMG(Params, Mesh, Mask, ElemMask, FileCreated, MaxBergVol DO i=1,iceberg - WRITE(36, '(A,i0,A,F20.0,A,F12.4,F12.4,F12.4,F12.4)') 'Iceberg ',i, ' Volume ', BergVolumes(i),& - ' Extent ', BergExtents(i*4-3:i*4) + WRITE(36, '(A,i0,A,F20.0,A,F20.4,F20.4,F20.4,F20.4,A,F20.4,F20.4,F20.4)') & + 'Iceberg ',i, ' Volume ', BergVolumes(i),& + ' Extent ', BergExtents(i*4-3:i*4), ' Centroid ', BergCentroids(i*3-2:i*3) END DO @@ -8086,6 +8131,8 @@ SUBROUTINE CheckFrontBoundary(Model, FrontConstraint, RightConstraint, LeftConst DO j=1, SIZE(ElNodes) IF(TopPerm(ElNodes(j)) /= 0) CYCLE IF(BottomPerm(ElNodes(j)) /= 0) CYCLE + IF(LeftPerm(ElNodes(j)) /= 0) CYCLE + IF(RightPerm(ElNodes(j)) /= 0) CYCLE Neighbours => Mesh % ParallelInfo % NeighbourList(ElNodes(j)) % Neighbours DO k=1, SIZE(Neighbours) IF(Neighbours(k) == ParEnv % MyPE) CYCLE @@ -8806,6 +8853,7 @@ SUBROUTINE SaveTerminusPosition(Model, Solver, Mesh, Boss) FirstTime=.TRUE. GotNode = .FALSE. counter = 0 + LastNode = 0 DO WHILE(LastNode /= FrontRight(1)) Found = .FALSE. IF(FirstTime) THEN diff --git a/elmerice/Solvers/CalvingGlacierAdvance3D.F90 b/elmerice/Solvers/CalvingGlacierAdvance3D.F90 index 019fe64cc3..cfe552d143 100644 --- a/elmerice/Solvers/CalvingGlacierAdvance3D.F90 +++ b/elmerice/Solvers/CalvingGlacierAdvance3D.F90 @@ -51,7 +51,7 @@ SUBROUTINE GlacierAdvance3D ( Model, Solver, dt, TransientSimulation ) !----------------------------------------------- TYPE(Mesh_t), POINTER :: Mesh TYPE(Nodes_t), TARGET :: FrontNodes - TYPE(ValueList_t), POINTER :: Params + TYPE(ValueList_t), POINTER :: Params, Material TYPE(Variable_t), POINTER :: Var, VeloVar, MeltVar, NormalVar, TangledVar, & DistVar ! TO DO clean all unused variables @@ -70,13 +70,14 @@ SUBROUTINE GlacierAdvance3D ( Model, Solver, dt, TransientSimulation ) MeltRate, Displace(3), y_coord(2), epsShift, LongRangeLimit, MaxDisplacement, & EpsTangle,thisEps,Shift, thisY,xx,yy,TempDist,MinDist,xt,yt,t, & a1(2), a2(2), b1(2), b2(2), b3(2), intersect(2), DistRailNode, RDisplace(3),& - buffer, VeloFactor + buffer, VeloFactor, zb REAL(KIND=dp), POINTER :: Advance(:) REAL(KIND=dp), ALLOCATABLE :: Rot_y_coords(:,:), Rot_z_coords(:,:), & xL(:),yL(:),xR(:),yR(:),xRail(:),yRail(:) LOGICAL :: Found, Debug, Parallel, Boss, ShiftLeft, LeftToRight, MovedOne, ShiftSecond, & Protrusion, SqueezeLeft, SqueezeRight, FirstTime=.TRUE., intersect_flag, FrontMelting, & - MovePastRailNode, HitsRails, Reverse, ThisBC, MoveBulk + MovePastRailNode, HitsRails, Reverse, ThisBC, MoveBulk, MoveBase,& + UnFoundFatal LOGICAL, ALLOCATABLE :: DangerZone(:), WorkLogical(:), & DontMove(:), FrontToRight(:), FrontToLeft(:) CHARACTER(LEN=MAX_NAME_LEN) :: SolverName, VeloVarName, MeltVarName, & @@ -206,6 +207,10 @@ SUBROUTINE GlacierAdvance3D ( Model, Solver, dt, TransientSimulation ) MoveBulk = ListGetLogical(Params,"MoveBulk", Found, DefValue=.FALSE.) IF(.NOT. Found) CALL Info(SolverName, "Not moving bulk as default") + MoveBase = ListGetLogical(Params,"Account for bedrock", Found, DefValue=.FALSE.) + ! still in testing + !IF(.NOT. Found) CALL Fatal(SolverName, "'Account for bedrock' not found") + !Get the front line FrontMaskName = "Calving Front Mask" TopMaskName = "Top Surface Mask" @@ -585,6 +590,35 @@ SUBROUTINE GlacierAdvance3D ( Model, Solver, dt, TransientSimulation ) END DO + !do we need to check base base of new advance location? + IF(MoveBase) THEN + + i = ListGetInteger( CurrentModel % Bodies(Mesh % Elements(1) % BodyId) % Values, & + 'Material') + Material => CurrentModel % Materials(i) % Values !TODO, this is not generalised + + DO i=1, Mesh % NumberOfNodes + IF(Perm(i) == 0) CYCLE + IF(FrontPerm(i) == 0) CYCLE + Mesh % Nodes % x(i) = Mesh % Nodes % x(i) + Advance((Perm(i)-1)*DOFs + 1) + Mesh % Nodes % y(i) = Mesh % Nodes % y(i) + Advance((Perm(i)-1)*DOFs + 2) + Mesh % Nodes % z(i) = Mesh % Nodes % z(i) + Advance((Perm(i)-1)*DOFs + 3) + zb = ListGetRealAtNode(Material, "Min Zs Bottom",i,UnfoundFatal=.TRUE.) + + Mesh % Nodes % x(i) = Mesh % Nodes % x(i) - Advance((Perm(i)-1)*DOFs + 1) + Mesh % Nodes % y(i) = Mesh % Nodes % y(i) - Advance((Perm(i)-1)*DOFs + 2) + Mesh % Nodes % z(i) = Mesh % Nodes % z(i) - Advance((Perm(i)-1)*DOFs + 3) + + IF(Mesh % Nodes % z(i) + Advance((Perm(i)-1)*DOFs + 3) < zb) THEN + Advance((Perm(i)-1)*DOFs + 3) = zb - Mesh % Nodes % z(i) + IF(Mesh % Nodes % z(i) < zb) THEN + Advance((Perm(i)-1)*DOFs + 3) = 0.0_dp + END IF + END IF + END DO + END IF + + ! find lateral boundary tags DO i=1,Model % NumberOfBCs ThisBC = ListGetLogical(Model % BCs(i) % Values,LeftMaskName,Found) diff --git a/elmerice/Solvers/CalvingRemeshMMG.F90 b/elmerice/Solvers/CalvingRemeshMMG.F90 index f34f933d96..0e65093407 100644 --- a/elmerice/Solvers/CalvingRemeshMMG.F90 +++ b/elmerice/Solvers/CalvingRemeshMMG.F90 @@ -79,8 +79,8 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) front_BC_ID, front_body_id, my_calv_front,calv_front, ncalv_parts, & group_calve, comm_calve, group_world,ecode, NElNodes, target_bodyid,gdofs(4), & PairCount,RPairCount, NCalvNodes, croot, nonCalvBoss,& - NVerts, NTetras, NPrisms, NTris, NQuads, NEdges - INTEGER, POINTER :: NodeIndexes(:), geom_id + NVerts, NTetras, NPrisms, NTris, NQuads, NEdges, dummyint + INTEGER, POINTER :: NodeIndexes(:), geom_id, TopPerm(:)=>NULL(), FrontPerm(:)=>NULL() INTEGER, ALLOCATABLE :: Prnode_count(:), cgroup_membs(:),disps(:), & PGDOFs_send(:),pcalv_front(:),GtoLNN(:),EdgePairs(:,:),REdgePairs(:,:), ElNodes(:),& Nodes(:), IslandCounts(:), pNCalvNodes(:,:), TetraQuality(:) @@ -91,15 +91,15 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) REAL(KIND=dp), POINTER :: WorkArray(:,:) => NULL() LOGICAL, ALLOCATABLE :: calved_node(:), remeshed_node(:), fixed_node(:), fixed_elem(:), & elem_send(:), RmElem(:), RmNode(:),new_fixed_node(:), new_fixed_elem(:), FoundNode(:,:), & - UsedElem(:), NewNodes(:), RmIslandNode(:), RmIslandElem(:), PartGotNodes(:) + UsedElem(:), NewNodes(:), RmIslandNode(:), RmIslandElem(:), PartGotNodes(:), SendNode(:) LOGICAL :: ImBoss, Found, Isolated, Debug,DoAniso,NSFail,CalvingOccurs,& RemeshOccurs,CheckFlowConvergence, NoNewNodes, RSuccess, Success,& SaveMMGMeshes, SaveMMGSols, PauseSolvers, PauseAfterCalving, FixNodesOnRails, & SolversPaused, NewIceberg, GotNodes(4), CalvingFileCreated=.FALSE., SuppressCalv,& - DistributedMesh,SaveTerminus + DistributedMesh,SaveTerminus,RemeshFront,RemeshIfNoCalving CHARACTER(LEN=MAX_NAME_LEN) :: SolverName, CalvingVarName, MeshName, SolName, & premmgls_meshfile, mmgls_meshfile, premmgls_solfile, mmgls_solfile,& - RepartMethod, Filename + RepartMethod, Filename, DistVarName TYPE(Variable_t), POINTER :: TimeVar INTEGER :: Time, remeshtimestep, proc, idx, island, node, MaxLSetIter, mmgloops REAL(KIND=dp) :: TimeReal, PreCalveVolume, PostCalveVolume, CalveVolume, LsetMinQuality @@ -183,7 +183,16 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) remesh_thresh = ListGetConstReal(SolverParams,"Remeshing Distance", DefValue=1000.0_dp) LsetMinQuality = ListGetConstReal(SolverParams,"Mesh Min Quality", DefValue=0.00001_dp) RmcValue = ListGetConstReal(SolverParams,"Mesh Rmc Value", DefValue=1e-15_dp) - CalvingVarName = ListGetString(SolverParams,"Calving Variable Name", DefValue="Calving Lset") + CalvingVarName = ListGetString(SolverParams, "Calving Variable Name", Found) + IF(.NOT. Found) THEN + CALL WARN(SolverName, "'Levelset Variable Name' not set so assuming 'Calving Lset'.") + CalvingVarName = "Calving LSet" + END IF + DistVarName = ListGetString(SolverParams, "Distance Variable Name", Found) + IF(.NOT. Found) THEN + CALL WARN(SolverName, "'Distance Variable Name' not set so assuming 'Distance'.") + DistVarName = "Distance" + END IF SaveMMGMeshes = ListGetLogical(SolverParams,"Save MMGLS Meshes", DefValue=.FALSE.) SaveMMGSols = ListGetLogical(SolverParams,"Save MMGLS Sols", DefValue=.FALSE.) IF(SaveMMGMeshes) THEN @@ -203,6 +212,16 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) FixNodesOnRails = ListGetLogical(SolverParams,"Fix Nodes On Rails", DefValue=.TRUE.) SuppressCalv = ListGetLogical(SolverParams,"Suppress Calving", DefValue=.FALSE.) SaveTerminus = ListGetLogical(SolverParams,"Save Terminus", DefValue=.TRUE.) + RemeshFront = ListGetLogical(SolverParams,"Remesh Full Calving Front", Found) + IF(.NOT. Found) THEN + CALL Info(SolverName, "Not Found 'Remesh Full Calving Front' asssuming TRUE") + RemeshFront = .TRUE. + END IF + RemeshIfNoCalving = ListGetLogical(SolverParams,"Remesh if no calving", Found) + IF(.NOT. Found) THEN + CALL Info(SolverName, "Not Found 'Remesh if no calving' asssuming TRUE") + RemeshIfNoCalving = .TRUE. + END IF ! calving algo passes through 202 elems to mmg i = ListGetInteger( Model % Bodies(Mesh % Elements(1) % BodyId) % Values, 'Material') @@ -238,7 +257,7 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) fixed_node = .FALSE. ! TO DO some other wah to define remeshed nodes - DistanceVar => VariableGet(Mesh % Variables, "Distance", .TRUE., UnfoundFatal=.TRUE.) + DistanceVar => VariableGet(Mesh % Variables, DistVarName, .TRUE., UnfoundFatal=.TRUE.) ALLOCATE(test_dist(NNodes)) test_dist = DistanceVar % Values(DistanceVar % Perm(:)) remeshed_node = test_dist < remesh_thresh @@ -246,7 +265,7 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) !Get the calving levelset function (-ve inside calving event, +ve in intact ice) !------------------- IF (CalvingOccurs) THEN - CalvingVar => VariableGet(Mesh % Variables, "Calving Lset", .TRUE., UnfoundFatal=.TRUE.) + CalvingVar => VariableGet(Mesh % Variables, CalvingVarName, .TRUE., UnfoundFatal=.TRUE.) ALLOCATE(test_lset(NNodes),& calved_node(NNodes)& @@ -257,7 +276,11 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) calved_node = test_lset < 0.0 ! calving front boundary nodes may have lset value greater than remesh dist DO i=1, NNodes - newdist = MINVAL((/test_lset(i), test_dist(i)/)) + IF(RemeshFront) THEN + newdist = MINVAL((/test_lset(i), test_dist(i)/)) + ELSE + newdist = test_lset(i) + END IF remeshed_node(i) = newdist < remesh_thresh END DO END IF @@ -308,9 +331,8 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) IF(elem_send(ParentElem % ElementIndex)) elem_send(i) = .TRUE. END DO - DEALLOCATE(fixed_node) !reuse later END IF - + DEALLOCATE(fixed_node) !reuse later !This does nothing yet but it will be important - determine !the discrete calving zones, each of which will be separately remeshed @@ -502,7 +524,11 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) ! so use dist so front boundary nodes aren't fixed IF(CalvingOccurs) THEN DO i=1,GNNode - newdist = MINVAL((/Gtest_lset(i), Gtest_dist(i)/)) + IF(RemeshFront) THEN + newdist = MINVAL((/Gtest_lset(i), Gtest_dist(i)/)) + ELSE + newdist = Gtest_lset(i) + END IF fixed_node(i) = newdist > remesh_thresh END DO ELSE @@ -518,6 +544,14 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) END IF END IF ! My calving front > 1 + ! if no calving and we don't want to remesh then go to end + ! need gathered mesh for terminus output + ! also need to update mesh deformation variables + IF(.NOT. RemeshIfNoCalving .AND. .NOT. CalvingOccurs .AND. .NOT. NSFail) THEN + RSuccess =.FALSE. + GO TO 30 + END IF + !Nominated partition does the remeshing IF(ImBoss) THEN IF (CalvingOccurs) THEN @@ -864,94 +898,97 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) NBulk = NewMeshR % NumberOfBulkElements NBdry = NewMeshR % NumberOfBoundaryElements - !sometimes some icebergs are not removed fully from the MMG levelset - !find all connected nodes in groups and remove in groups not in the - !main icebody - !limit here of 10 possible mesh 'islands' - ALLOCATE(FoundNode(10, NNodes), ElNodes(4), & - UsedElem(NBulk)) - FoundNode = .FALSE.! count of nodes found - UsedElem = .FALSE. !count of elems used - island=0 ! count of different mesh islands - NoNewNodes=.TRUE. ! whether node has neighour - DO WHILE(COUNT(FoundNode) < NNodes) - IF(NoNewNodes) THEN - island = island + 1 - NewIceberg = .TRUE. - END IF - NoNewNodes = .TRUE. - DO i=1, NBulk - IF(UsedElem(i)) CYCLE - Element => NewMeshR % Elements(i) - ElNodes = Element % NodeIndexes - ! if there are not any matching nodes and its not a new iceberg - DO j=1, 4 - GotNodes(j) = ANY(FoundNode(:, ElNodes(j))) + IF(RemeshFront) THEN ! if just remesh where calving occurs cannot remove + !islands as mesh already split + !sometimes some icebergs are not removed fully from the MMG levelset + !find all connected nodes in groups and remove in groups not in the + !main icebody + !limit here of 10 possible mesh 'islands' + ALLOCATE(FoundNode(10, NNodes), ElNodes(4), & + UsedElem(NBulk)) + FoundNode = .FALSE.! count of nodes found + UsedElem = .FALSE. !count of elems used + island=0 ! count of different mesh islands + NoNewNodes=.TRUE. ! whether node has neighour + DO WHILE(COUNT(FoundNode) < NNodes) + IF(NoNewNodes) THEN + island = island + 1 + NewIceberg = .TRUE. + END IF + NoNewNodes = .TRUE. + DO i=1, NBulk + IF(UsedElem(i)) CYCLE + Element => NewMeshR % Elements(i) + ElNodes = Element % NodeIndexes + ! if there are not any matching nodes and its not a new iceberg + DO j=1, 4 + GotNodes(j) = ANY(FoundNode(:, ElNodes(j))) + END DO + IF(ALL(.NOT. GotNodes) .AND. .NOT. NewIceberg) CYCLE + NewIceberg = .FALSE. + UsedElem(i) = .TRUE. + FoundNode(island, ElNodes) = .TRUE. + NoNewNodes = .FALSE. END DO - IF(ALL(.NOT. GotNodes) .AND. .NOT. NewIceberg) CYCLE - NewIceberg = .FALSE. - UsedElem(i) = .TRUE. - FoundNode(island, ElNodes) = .TRUE. - NoNewNodes = .FALSE. END DO - END DO - ALLOCATE(IslandCounts(10)) - DO i=1,10 - IslandCounts(i) = COUNT(FoundNode(i, :)) - END DO - Island = COUNT(IslandCounts > 0) - DEALLOCATE(IslandCounts) ! reused - - ! if iceberg not removed mark nodes and elems - IF(Island > 1) THEN - CALL WARN(SolverName, 'Mesh island found after levelset- removing...') - ALLOCATE(RmIslandNode(NNodes), RmIslandElem(Nbulk+NBdry),& - IslandCounts(Island)) - RmIslandNode = .FALSE. - RmIslandElem = .FALSE. - counter=0 - DO i=1, 10 - IF(COUNT(FoundNode(i, :)) == 0) CYCLE - counter=counter+1 - IslandCounts(counter) = COUNT(FoundNode(i, :)) - END DO - DO i=1, Island - IF(IslandCounts(i) < MAXVAL(IslandCounts)) THEN - DO j=1, NNodes - IF(.NOT. FoundNode(i,j)) CYCLE! if not part of island - RmIslandNode(j) = .TRUE. - END DO - END IF - END DO - ! mark all elems with rm nodes as need removing - nodes = PACK((/ (i,i=1,SIZE(RmIslandNode)) /), RmIslandNode .eqv. .TRUE.) - DO i=1, Nbulk+Nbdry - Element => NewMeshR % Elements(i) - ElNodes = Element % NodeIndexes - !any([(any(A(i) == B), i = 1, size(A))]) - IF( .NOT. ANY([(ANY(ElNodes(i)==Nodes), i = 1, SIZE(ElNodes))])) CYCLE - RmIslandElem(i) = .TRUE. + ALLOCATE(IslandCounts(10)) + DO i=1,10 + IslandCounts(i) = COUNT(FoundNode(i, :)) END DO + Island = COUNT(IslandCounts > 0) + DEALLOCATE(IslandCounts) ! reused + + ! if iceberg not removed mark nodes and elems + IF(Island > 1) THEN + CALL WARN(SolverName, 'Mesh island found after levelset- removing...') + ALLOCATE(RmIslandNode(NNodes), RmIslandElem(Nbulk+NBdry),& + IslandCounts(Island)) + RmIslandNode = .FALSE. + RmIslandElem = .FALSE. + counter=0 + DO i=1, 10 + IF(COUNT(FoundNode(i, :)) == 0) CYCLE + counter=counter+1 + IslandCounts(counter) = COUNT(FoundNode(i, :)) + END DO + DO i=1, Island + IF(IslandCounts(i) < MAXVAL(IslandCounts)) THEN + DO j=1, NNodes + IF(.NOT. FoundNode(i,j)) CYCLE! if not part of island + RmIslandNode(j) = .TRUE. + END DO + END IF + END DO + ! mark all elems with rm nodes as need removing + nodes = PACK((/ (i,i=1,SIZE(RmIslandNode)) /), RmIslandNode .eqv. .TRUE.) + DO i=1, Nbulk+Nbdry + Element => NewMeshR % Elements(i) + ElNodes = Element % NodeIndexes + !any([(any(A(i) == B), i = 1, size(A))]) + IF( .NOT. ANY([(ANY(ElNodes(i)==Nodes), i = 1, SIZE(ElNodes))])) CYCLE + RmIslandElem(i) = .TRUE. + END DO - !Chop out missed iceberg if they exist - CALL CutMesh(NewMeshR, RmIslandNode, RmIslandElem) + !Chop out missed iceberg if they exist + CALL CutMesh(NewMeshR, RmIslandNode, RmIslandElem) - !modify rmelem and rmnode to include island removal - counter=0 - DO i=1, SIZE(RmElem) - IF(RmElem(i)) CYCLE - counter=counter+1 - IF(RmIslandElem(Counter)) RmElem(i) =.TRUE. - END DO - counter=0 - DO i=1, SIZE(RmNode) - IF(RmNode(i)) CYCLE - counter=counter+1 - IF(RmIslandNode(Counter)) RmNode(i) =.TRUE. - END DO - CALL INFO(SolverName, 'Mesh island removed', level=10) - END IF + !modify rmelem and rmnode to include island removal + counter=0 + DO i=1, SIZE(RmElem) + IF(RmElem(i)) CYCLE + counter=counter+1 + IF(RmIslandElem(Counter)) RmElem(i) =.TRUE. + END DO + counter=0 + DO i=1, SIZE(RmNode) + IF(RmNode(i)) CYCLE + counter=counter+1 + IF(RmIslandNode(Counter)) RmNode(i) =.TRUE. + END DO + CALL INFO(SolverName, 'Mesh island removed', level=10) + END IF + END IF !remeshfront END IF ! CalvingOccurs @@ -1047,11 +1084,12 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) 30 CONTINUE !Wait for all partitions to finish - IF(My_Calv_Front>0) THEN - CALL MPI_BARRIER(COMM_CALVE,ierr) - CALL MPI_COMM_FREE(COMM_CALVE,ierr) - CALL MPI_GROUP_FREE(group_world,ierr) - END IF + + CALL MPI_GROUP_FREE(group_world,ierr) + CALL MPI_GROUP_FREE(group_calve,ierr) + IF(My_Calv_Front > 0) & + CALL MPI_COMM_FREE(COMM_CALVE,ierr) + CALL MPI_BARRIER(ELMER_COMM_WORLD,ierr) CALL MPI_BCAST(CalvingFileCreated, 1, MPI_LOGICAL, my_cboss, ELMER_COMM_WORLD, ierr) @@ -1122,7 +1160,77 @@ SUBROUTINE CalvingRemeshMMG( Model, Solver, dt, Transient ) END IF END DO - IF(SaveTerminus) CALL SaveTerminusPosition(Model, Solver, GatheredMesh, ImBoss) + IF(SaveTerminus) THEN + IF(RemeshFront) THEN ! got entire front + CALL SaveTerminusPosition(Model, Solver, GatheredMesh, ImBoss) + ELSE ! only remeshing where bergs calve + CALL MakePermUsingMask( Model, Solver, GatheredMesh, "Calving Front Mask", & + .FALSE., FrontPerm, dummyint) + CALL MakePermUsingMask( Model, Solver, GatheredMesh, "Top Surface Mask", & + .FALSE., TopPerm, dummyint) + + IF(.NOT. ASSOCIATED(GatheredMesh % Repartition)) THEN + ALLOCATE(GatheredMesh % Repartition(GatheredMesh % NumberOfBulkElements+& + GatheredMesh % NumberOfBoundaryElements), STAT=ierr) + IF(ierr /= 0) PRINT*, ParEnv % MyPE, 'Unable to allocate mesh % repartition' + END IF + + GatheredMesh % Repartition = ParEnv % MyPE + 1 + IF(ImBoss) DEALLOCATE(fixed_node) + DEALLOCATE(elem_send) + ALLOCATE(SendNode(GatheredMesh % NumberOfNodes), & + fixed_node(GatheredMesh % NumberOfNodes)) + SendNode = .FALSE.; fixed_node = .FALSE. + DO i=1, GatheredMesh % NumberOfNodes + IF( (TopPerm(i)>0) .AND. (FrontPerm(i)>0)) THEN + SendNode(i) = .TRUE. + END IF + END DO + + ALLOCATE(elem_send(GatheredMesh % NumberOfBulkElements+& + GatheredMesh % NumberOfBoundaryElements)) + elem_send = .FALSE. + DO i=1,GatheredMesh % NumberOfBulkElements + Element => GatheredMesh % Elements(i) + CALL Assert(Element % ElementIndex == i, SolverName,"Misnumbered bulk element.") + NodeIndexes => Element % NodeIndexes + NElNodes = Element % TYPE % NumberOfNodes + IF(ANY(SendNode(NodeIndexes(1:NElNodes)))) THEN + elem_send(i) = .TRUE. + IF(.NOT. ALL(SendNode(NodeIndexes(1:NElNodes)))) THEN + fixed_node(NodeIndexes(1:NElNodes)) = .TRUE. + END IF + END IF + END DO + + SendNode = SendNode .OR. fixed_node + + !Cycle boundary elements, checking parent elems + !BC elements follow parents + DO i=GatheredMesh % NumberOfBulkElements+1, GatheredMesh % NumberOfBulkElements+& + GatheredMesh % NumberOfBoundaryElements + Element => GatheredMesh % Elements(i) + ParentElem => Element % BoundaryInfo % Left + IF(.NOT. ASSOCIATED(ParentElem)) THEN + ParentElem => Element % BoundaryInfo % Right + END IF + CALL Assert(ASSOCIATED(ParentElem),SolverName,"Boundary element has no parent!") + IF(elem_send(ParentElem % ElementIndex)) elem_send(i) = .TRUE. + END DO + + DO i=1, GatheredMesh % NumberOfBulkElements + GatheredMesh % NumberOfBoundaryElements + IF(Elem_send(i)) GatheredMesh % RePartition(i) = my_cboss + 1 + END DO + + ParMetisMesh => RedistributeMesh(Model, GatheredMesh, .TRUE., .FALSE.) + + CALL SaveTerminusPosition(Model, Solver, ParMetisMesh, ImBoss) + + CALL ReleaseMesh(ParMetisMesh) + ParMetisMesh => NULL() + DEALLOCATE(FrontPerm,TopPerm) + END IF ! remeshfront + END IF ! saveterminus !Call zoltan to determine redistribution of mesh ! then do the redistribution diff --git a/elmerice/Solvers/GroundedSolver.F90 b/elmerice/Solvers/GroundedSolver.F90 index e33b093806..4f60de929e 100644 --- a/elmerice/Solvers/GroundedSolver.F90 +++ b/elmerice/Solvers/GroundedSolver.F90 @@ -80,7 +80,7 @@ SUBROUTINE GroundedSolver( Model,Solver,dt,TransientSimulation ) TYPE(Nodes_t), SAVE :: Nodes LOGICAL :: AllocationsDone = .FALSE., GotIt, stat,UnFoundFatal=.TRUE.,& - AllGrounded = .FALSE., useLSvar = .FALSE. + AllGrounded = .FALSE., useLSvar = .FALSE., Active INTEGER :: i, mn, n, t, Nn, istat, DIM, MSum, ZSum, bedrockSource INTEGER, POINTER :: Permutation(:), bedrockPerm(:), LSvarPerm(:) @@ -103,6 +103,8 @@ SUBROUTINE GroundedSolver( Model,Solver,dt,TransientSimulation ) Permutation => PointerToVariable % Perm VariableValues => PointerToVariable % Values + Active = ANY(Permutation > 0) + CALL INFO(SolverName, 'Computing grounded mask from geometry', level=3) !-------------------------------------------------------------- @@ -281,7 +283,7 @@ SUBROUTINE GroundedSolver( Model,Solver,dt,TransientSimulation ) END IF END DO - IF ( ParEnv % PEs>1 ) CALL ParallelSumVector( Solver % Matrix, VariableValues, 1 ) + IF ( ParEnv % PEs>1 .AND. Active) CALL ParallelSumVector( Solver % Matrix, VariableValues, 1 ) CALL INFO( SolverName , 'Done') diff --git a/elmerice/Solvers/ProjectCalving.F90 b/elmerice/Solvers/ProjectCalving.F90 index ffa283ec33..a81a23f9b7 100644 --- a/elmerice/Solvers/ProjectCalving.F90 +++ b/elmerice/Solvers/ProjectCalving.F90 @@ -619,7 +619,7 @@ SUBROUTINE ProjectCalving( Model,Solver,dt,TransientSimulation ) myNodes = myNodes + (PlaneNodes-Parenv % Pes*myNodes) IF ( Parenv % PEs>1 ) THEN - CALL CheckBuffer( (2 * (2 * (DOFs_3D+1)*SUM(PointStore(:) % NoInt)) + PlaneNodes )) + CALL CheckBuffer( (4 * (2 * (DOFs_3D+1)*SUM(PointStore(:) % NoInt)) + PlaneNodes )) CALL SendPoints() CALL ReceivePoints() CALL CheckBuffer(1) !Just frees the buffer space for future subroutines. diff --git a/fem/src/MeshPartition.F90 b/fem/src/MeshPartition.F90 index 47e16f25b4..805735de94 100644 --- a/fem/src/MeshPartition.F90 +++ b/fem/src/MeshPartition.F90 @@ -471,9 +471,27 @@ SUBROUTINE Zoltan_Interface( Model, Mesh, SerialMode, NoPartitions, PartitionCan CALL Info(FuncName, Message) IF(ImbalanceTol < 1.0) CALL FATAL(FuncName, 'Unable to rebalance successfully') DEALLOCATE(ElemAdj, ElemAdjProc, ElemStart, PartSuccess) + + ! release zoltan input/output arrays + IF(numImport > 0) & + zierr = Zoltan_LB_Free_Part(importGlobalGids,importLocalGids,importProcs,importToPart) + IF(numExport > 0) & + zierr = Zoltan_LB_Free_Part(exportGlobalGids,exportLocalGids,exportProcs,exportToPart) + + ! release zoltan object and mpi communicators + CALL Zoltan_Destroy(zz_obj) GOTO 10 END IF END IF + + ! release zoltan input/output arrays + IF(numImport > 0) & + zierr = Zoltan_LB_Free_Part(importGlobalGids,importLocalGids,importProcs,importToPart) + IF(numExport > 0) & + zierr = Zoltan_LB_Free_Part(exportGlobalGids,exportLocalGids,exportProcs,exportToPart) + + ! release zoltan object and mpi communicators + CALL Zoltan_Destroy(zz_obj) CALL Info(FuncName,'Finished Zoltan partitioning',Level=10) @@ -620,11 +638,12 @@ SUBROUTINE GlobalElemAdjacency( Mesh, ElemAdj, ElemAdjProc, ElemStart, DIM ) NBulk = Mesh % NumberOfBulkElements !Find and globally number mesh faces - IF(Nbulk == 0 ) THEN - CONTINUE - - ELSE IF(DIM == 3) THEN - CALL FindMeshFaces3D(Mesh) + IF(DIM == 3) THEN + IF( NBulk /= 0 ) THEN + CALL FindMeshFaces3D(Mesh) + ELSE + CALL AllocateVector( Mesh % Faces, 6*Mesh % NumberOfBulkElements, 'FindMeshFaces3D' ) + END IF CALL FindMeshEdges3D(Mesh) CALL SParFaceNumbering(Mesh) MFacePtr => Mesh % Faces @@ -687,7 +706,7 @@ SUBROUTINE GlobalElemAdjacency( Mesh, ElemAdj, ElemAdjProc, ElemStart, DIM ) !don't know at this point the required size of work_int, if there are lots !of bulk elements, probably NBulk is sufficient, but for only a few !disconnected elements, maybe not. So set min size = 1000 - work_size = MAX(NBulk, 1000) + work_size = MAX(NBulk, 5000) ALLOCATE(SendFaces(ParEnv % PEs),RecvFaces(ParEnv % PEs),work_int(work_size)) RecvFaces % Count = 0 @@ -2350,7 +2369,7 @@ END FUNCTION ElementPartitions !> sending to another partition. !------------------------------------------------------------------------------ SUBROUTINE PackMeshPieces(Model, Mesh, NewPart, ParallelMesh, NoPartitions, & - SentPack, dim, NodalVals ) + SentPack, dim, NodalVals ) IMPLICIT NONE @@ -3041,6 +3060,7 @@ SUBROUTINE UnpackMeshPieces(Model, Mesh, NewMesh, NewPart, & REAL(KIND=dp), POINTER, OPTIONAL :: NodalVals(:,:) !---------------------------------- TYPE( MeshPack_t), ALLOCATABLE, TARGET :: RecPack(:) + !---------------------------------- TYPE(Element_t), POINTER :: Element, Element0 INTEGER :: i,j,k,n,t,nbulk,nbdry,allocstat,part,elemcode,elemindex,geom_id,sweep,partindex INTEGER :: gind,lind,rcount,icount,lcount,minelem,maxelem,newnbdry,newnodes,newnbulk @@ -3327,8 +3347,7 @@ SUBROUTINE UnpackMeshPieces(Model, Mesh, NewMesh, NewPart, & END IF IF( GlobalToLocal(k) <= 0 .OR. GlobalToLocal(k) > newnodes ) THEN errcount = errcount + 1 - PRINT *,'Local index out of bounds:',ParEnv % Mype,Element % PartIndex, k,minind,maxind,& - SIZE(GlobalToLocal), GlobalToLocal(k) + PRINT *,'Local index out of bounds:',ParEnv % Mype,Element % PartIndex, k,minind,maxind,GlobalToLocal(k) END IF Element % NodeIndexes(j) = GlobalToLocal(k) END DO @@ -3438,12 +3457,7 @@ SUBROUTINE UnpackMeshPieces(Model, Mesh, NewMesh, NewPart, & NewMesh % Nodes % y(k) = PPack % rdata(rcount+2) IF( dim == 3 ) NewMesh % Nodes % z(k) = PPack % rdata(rcount+3) rcount = rcount + dim - - IF(nVals > 0 ) THEN - NewVals(k,1:nVals) = PPack % rdata(rcount+1:rcount+nVals) - rcount = rcount + nVals - END IF - + NewMesh % ParallelInfo % GInterface(k) = PPack % ldata(lcount+1) lcount = lcount + 1 END DO @@ -3452,11 +3466,14 @@ SUBROUTINE UnpackMeshPieces(Model, Mesh, NewMesh, NewPart, & IF( errcount > 0 ) THEN CALL Fatal(Caller,'Encountered '//I2S(errcount)//' indexing issues in nodes') END IF - + + + n = COUNT( NewMesh % ParallelInfo % GInterface ) CALL Info(Caller,'Potential interface nodes '//I2S(n)//' out of '& //I2S(NewMesh % NumberOfNodes),Level=20) - + + CALL Info(Caller,'Creating local to global numbering for '& //I2S(newnodes)//' nodes',Level=20) ALLOCATE( NewMesh % ParallelInfo % GlobalDofs( newnodes ), STAT = allocstat) diff --git a/fem/src/MeshRemeshing.F90 b/fem/src/MeshRemeshing.F90 index 9d302a215b..809fe80330 100644 --- a/fem/src/MeshRemeshing.F90 +++ b/fem/src/MeshRemeshing.F90 @@ -560,9 +560,9 @@ SUBROUTINE Set_PMMG_Parameters(SolverParams, ReTrial ) ! If this is a ReTrial then we only change some real valued keywords! - IF( PRESENT( ReTrial ) ) THEN - IF( ReTrial ) RETURN - END IF + !IF( PRESENT( ReTrial ) ) THEN + ! IF( ReTrial ) RETURN + !END IF !!! PARAMS: generic options (debug, mem, verbosity) ! [val] Set the verbosity level to n diff --git a/fem/src/MeshUtils.F90 b/fem/src/MeshUtils.F90 index df8dc8ece3..25bbbe1bb3 100644 --- a/fem/src/MeshUtils.F90 +++ b/fem/src/MeshUtils.F90 @@ -1409,6 +1409,8 @@ SUBROUTINE EnlargeCoordinates(Mesh) Mesh % MaxBDOFs * Mesh % NumberOFBulkElements n0 = SIZE( Mesh % Nodes % x ) + IF(.NOT. ASSOCIATED(Mesh % Nodes % x)) n0 = 0 + pelementsPresent = .FALSE. DO i=1,Mesh % NumberOfBulkElements IF(isPelement(Mesh % Elements(i))) THEN @@ -24555,6 +24557,14 @@ SUBROUTINE ReleaseMesh( Mesh ) IF ( ASSOCIATED( Mesh % ParallelInfo % FaceInterface ) ) & DEALLOCATE( Mesh % ParallelInfo % FaceInterface ) + IF ( ASSOCIATED( Mesh % ParallelInfo % FaceNeighbourList ) ) THEN + DO i=1,Mesh % NumberOfNodes + IF(ASSOCIATED( Mesh % ParallelInfo % FaceNeighbourList(i) % Neighbours ) ) & + DEALLOCATE( Mesh % ParallelInfo % FaceNeighbourList(i) % Neighbours ) + END DO + DEALLOCATE( Mesh % ParallelInfo % FaceNeighbourList ) + END IF + IF ( ASSOCIATED( Mesh % ParallelInfo % EdgeNeighbourList ) ) THEN DO i=1,Mesh % NumberOfNodes IF(ASSOCIATED( Mesh % ParallelInfo % EdgeNeighbourList(i) % Neighbours ) ) &