From 10b6f5977eb8baad1f9a5a79f676cca765c92b0c Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Fri, 30 Aug 2019 16:04:07 -0600 Subject: [PATCH 1/9] Adding boundary conditions info whe printing the ascii mesh. Minor bug fix for floating condition --- .../Interface_velocity_solver.cpp | 39 +++++++++++++------ 1 file changed, 28 insertions(+), 11 deletions(-) diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.cpp b/src/core_landice/mode_forward/Interface_velocity_solver.cpp index 6612b7bba0..ca3e511c54 100644 --- a/src/core_landice/mode_forward/Interface_velocity_solver.cpp +++ b/src/core_landice/mode_forward/Interface_velocity_solver.cpp @@ -1545,7 +1545,7 @@ void import2DFields(std::map bdExtensionMap, double const* bedTopograp // Check that this margin location is not below sea level! if (bedTopographyData[iV] < 0.0) { // This check is probably redundant... thicknessData[iV] = eps*3.0; // insert special value here to make identifying these points easier in exo output - elevationData[iV] = (rho_ocean / rho_ice - 1.0) * thicknessData[iV]; // floating surface + elevationData[iV] = (1.0 - rho_ice / rho_ocean) * thicknessData[iV]; // floating surface } } } else { @@ -1561,7 +1561,7 @@ void import2DFields(std::map bdExtensionMap, double const* bedTopograp //for (int i = 0; i < dirichletNodesIDs.size(); i++) // std::cout << dirichletNodesIDs.at(i) << ' '; // print entire list of Diri nodes for debugging. thicknessData[iV] = eps*2.0; // insert special small value here to make identifying these points easier in exo output - elevationData[iV] = (rho_ocean / rho_ice - 1.0) * thicknessData[iV]; // floating surface + elevationData[iV] = (1.0 - rho_ice / rho_ocean) * thicknessData[iV]; // floating surface } } // if below sea level } // floating or not @@ -2157,21 +2157,32 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex) std::cout << "Writing mesh to albany.msh." << std::endl; // msh file std::ofstream outfile; + outfile.precision(15); outfile.open ("albany.msh", std::ios::out | std::ios::trunc); if (outfile.is_open()) { int nVerticesBoundaryEdge = 0; - for (int index = 0; index < nVertices; index++) { + for (int index = 0; index < nEdges; index++) { if (isBoundaryEdge[index]) nVerticesBoundaryEdge += 1; } - std::vector verticesOnBoundaryEdge; - verticesOnBoundaryEdge.resize(2 * nVerticesBoundaryEdge); + + //creating set from vector so that we can find elements in the set in log time + std::set floatingEdgesSet; + for (int i = 0; i < floatingEdgesIds.size(); i++) + floatingEdgesSet.insert(floatingEdgesIds[i]); + + + std::vector boundaryEdges; //list of edge vertices and edge label + boundaryEdges.resize(3 * nVerticesBoundaryEdge); int iVerticesBoundaryEdge = 0; - for (int index = 0; index < nVertices; index++) { + for (int index = 0; index < nEdges; index++) { if (isBoundaryEdge[index]) { - verticesOnBoundaryEdge[0 + 2 * iVerticesBoundaryEdge] = verticesOnEdge[0 + 2 * index]; - verticesOnBoundaryEdge[1 + 2 * iVerticesBoundaryEdge] = verticesOnEdge[1 + 2 * index]; - iVerticesBoundaryEdge += 1; + boundaryEdges[0 + 3 * iVerticesBoundaryEdge] = verticesOnEdge[0 + 2 * index]; + boundaryEdges[1 + 3 * iVerticesBoundaryEdge] = verticesOnEdge[1 + 2 * index]; + auto search = floatingEdgesSet.find(indexToEdgeID[index]); //slow but executed only when printing + //boundary edges labels: 2 if floating, 1 otherwise + boundaryEdges[2 + 3 * iVerticesBoundaryEdge] = (search != floatingEdgesSet.end()) ? 2 : 1; + iVerticesBoundaryEdge ++; } } //std::cout<<"final count: "< Date: Fri, 8 Nov 2019 14:01:53 -0700 Subject: [PATCH 2/9] Modified the logic to determined whether the finite element mesh needs to be recomputed or not The mesh needs to be recomputed non only of the vertices ice mask has changed, but also if the floating edges mask and the cell ice mask has changed. Note: in order to compare the clee mask with the previous time steps, we would need the cell mask to have 2 time levels. Howver, if we set the cell mask to have 2 time levels we get different results. This needs to be fixed. --- src/core_landice/Registry.xml | 11 +++- .../analysis_members/mpas_li_global_stats.F | 2 +- .../analysis_members/mpas_li_regional_stats.F | 2 +- .../mode_forward/mpas_li_advection.F | 4 +- .../mode_forward/mpas_li_calving.F | 20 +++---- .../mode_forward/mpas_li_diagnostic_vars.F | 3 +- .../mode_forward/mpas_li_iceshelf_melt.F | 4 +- src/core_landice/mode_forward/mpas_li_sia.F | 2 +- .../mode_forward/mpas_li_statistics.F | 2 +- .../mode_forward/mpas_li_subglacial_hydro.F | 16 +++--- .../mode_forward/mpas_li_thermal.F | 2 +- .../mode_forward/mpas_li_velocity.F | 57 ++++++++++++++++--- .../mode_forward/mpas_li_velocity_external.F | 14 +++-- 13 files changed, 95 insertions(+), 44 deletions(-) diff --git a/src/core_landice/Registry.xml b/src/core_landice/Registry.xml index 7062959278..c534626142 100644 --- a/src/core_landice/Registry.xml +++ b/src/core_landice/Registry.xml @@ -931,7 +931,7 @@ is the value of that variable from the *previous* time level! - + + diff --git a/src/core_landice/analysis_members/mpas_li_global_stats.F b/src/core_landice/analysis_members/mpas_li_global_stats.F index 41ddc9aebd..b367875979 100644 --- a/src/core_landice/analysis_members/mpas_li_global_stats.F +++ b/src/core_landice/analysis_members/mpas_li_global_stats.F @@ -283,7 +283,7 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel=1) call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied) call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied) call mpas_pool_get_array(geometryPool, 'basalMassBal', basalMassBal) diff --git a/src/core_landice/analysis_members/mpas_li_regional_stats.F b/src/core_landice/analysis_members/mpas_li_regional_stats.F index 0e0543df57..4c6df7def9 100644 --- a/src/core_landice/analysis_members/mpas_li_regional_stats.F +++ b/src/core_landice/analysis_members/mpas_li_regional_stats.F @@ -275,7 +275,7 @@ subroutine li_compute_regional_stats(domain, memberName, timeLevel, err) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel=1) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied) call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied) diff --git a/src/core_landice/mode_forward/mpas_li_advection.F b/src/core_landice/mode_forward/mpas_li_advection.F index 4b728d3d29..4c0576ab89 100644 --- a/src/core_landice/mode_forward/mpas_li_advection.F +++ b/src/core_landice/mode_forward/mpas_li_advection.F @@ -257,7 +257,7 @@ subroutine li_advection_thickness_tracers(& call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) call mpas_pool_get_array(geometryPool, 'layerThicknessEdge', layerThicknessEdge) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'dynamicThickening', dynamicThickening) call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness) @@ -630,7 +630,7 @@ subroutine li_update_geometry(geometryPool) call mpas_pool_get_dimension(geometryPool, 'nCells', nCells) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface) call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface) diff --git a/src/core_landice/mode_forward/mpas_li_calving.F b/src/core_landice/mode_forward/mpas_li_calving.F index 41e0523d8d..ac7d032ab0 100644 --- a/src/core_landice/mode_forward/mpas_li_calving.F +++ b/src/core_landice/mode_forward/mpas_li_calving.F @@ -402,7 +402,7 @@ subroutine li_restore_calving_front(domain, err) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'restoreThickness', restoreThickness) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) ! get required fields from the thermal pool call mpas_pool_get_array(thermalPool, 'temperature', temperature) @@ -870,7 +870,7 @@ subroutine floating_calving(domain, calvingFraction, err) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'thickness', thickness) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) calvingThickness = 0.0_RKIND @@ -929,7 +929,7 @@ subroutine remove_small_islands(meshPool, geometryPool) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'thickness', thickness) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) do iCell = 1, nCellsSolve @@ -1039,7 +1039,7 @@ subroutine topographic_calving(domain, calvingFraction, err) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) calvingThickness = 0.0_RKIND @@ -1152,7 +1152,7 @@ subroutine eigencalving(domain, err) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) call mpas_pool_get_array(geometryPool, 'eigencalvingParameter', eigencalvingParameter) @@ -1336,7 +1336,7 @@ subroutine distribute_calving_flux(meshPool, geometryPool, scratchPool, calvingF ! get fields call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_array(meshPool, 'deltat', deltat) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'requiredCalvingVolumeRate', requiredCalvingVolumeRate) @@ -1566,7 +1566,7 @@ subroutine calculate_calving_front_mask(meshPool, geometryPool, calvingFrontMask call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) calvingFrontMask = 0 !initialize @@ -1656,7 +1656,7 @@ subroutine remove_icebergs(domain) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'contiguousDynamicIce', contiguousDynamicIce) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(geometryPool, 'nCellsSolve', nCellsSolve) @@ -1710,7 +1710,7 @@ subroutine remove_icebergs(domain) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'contiguousDynamicIce', contiguousDynamicIce) call mpas_pool_get_array(geometryPool, 'contiguousDynamicIceOld', contiguousDynamicIceOld) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) @@ -1765,7 +1765,7 @@ subroutine remove_icebergs(domain) block => domain % blocklist do while (associated(block)) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'contiguousDynamicIce', contiguousDynamicIce) diff --git a/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F b/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F index ab47f0d226..a8f2664137 100644 --- a/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F +++ b/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F @@ -258,10 +258,9 @@ subroutine li_calculate_apparent_diffusivity(meshPool, velocityPool, scratchPool call mpas_pool_get_array(geometryPool, 'thickness', thickness, timeLevel=1) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness, timeLevel=1) call mpas_pool_get_array(geometryPool, 'apparentDiffusivity', apparentDiffusivity) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(velocityPool, 'uReconstructZonal', uReconstructAxis1) call mpas_pool_get_array(velocityPool, 'uReconstructMeridional', uReconstructAxis2) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_field(scratchPool, 'workCell', slopeReconstructXField) call mpas_allocate_scratch_field(slopeReconstructXField, .true.) call mpas_pool_get_field(scratchPool, 'workCell2', slopeReconstructYField) diff --git a/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F b/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F index 95f6716fa8..0b99e9d843 100644 --- a/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F +++ b/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F @@ -247,7 +247,7 @@ subroutine li_basal_melt_floating_ice(domain, err) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal) ! get fields from the scratch pool @@ -853,7 +853,7 @@ subroutine calc_iceshelf_draft_info(domain, GLdepth, CFdepth, err) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) call mpas_pool_get_array(meshPool, 'areaCell', areaCell) diff --git a/src/core_landice/mode_forward/mpas_li_sia.F b/src/core_landice/mode_forward/mpas_li_sia.F index b22001c4f2..d400c12e35 100644 --- a/src/core_landice/mode_forward/mpas_li_sia.F +++ b/src/core_landice/mode_forward/mpas_li_sia.F @@ -316,7 +316,7 @@ subroutine li_sia_solve(meshPool, geometryPool, thermalPool, velocityPool, err) call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'slopeEdge', slopeEdge) call mpas_pool_get_array(geometryPool, 'normalSlopeEdge', normalSlopeEdge) call mpas_pool_get_array(thermalPool, 'temperature', temperature) diff --git a/src/core_landice/mode_forward/mpas_li_statistics.F b/src/core_landice/mode_forward/mpas_li_statistics.F index 47a7be44e8..99de7c4939 100644 --- a/src/core_landice/mode_forward/mpas_li_statistics.F +++ b/src/core_landice/mode_forward/mpas_li_statistics.F @@ -235,7 +235,7 @@ subroutine li_compute_statistics(domain, itimestep) ! Geometry arrays call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'thickness', thickness, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness, timeLevel = 1) diff --git a/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F b/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F index 9e2e4f79ca..596119d01c 100644 --- a/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F +++ b/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F @@ -163,7 +163,7 @@ subroutine li_SGH_init(domain, err) waterPressure = max(0.0_RKIND, waterPressure) waterPressure = min(waterPressure, rhoi * gravity * thickness) ! set pressure correctly under floating ice and open ocean - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) where ( (li_mask_is_floating_ice(cellMask)) .or. & ( (.not. li_mask_is_ice(cellMask)) .and. (bedTopography < config_sea_level) ) ) @@ -379,7 +379,7 @@ subroutine li_SGH_solve(domain, err) call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH) call mpas_pool_get_array(hydroPool, 'basalMeltInput', basalMeltInput) call mpas_pool_get_array(hydroPool, 'externalWaterInput', externalWaterInput) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) WtillOld = Wtill Wtill = Wtill + deltatSGH * ( (basalMeltInput + externalWaterInput) / rho_water - Cd) @@ -420,7 +420,7 @@ subroutine li_SGH_solve(domain, err) call mpas_pool_get_array(hydroPool, 'waterVelocity', waterVelocity) call mpas_pool_get_array(hydroPool, 'waterVelocityCellX', waterVelocityCellX) call mpas_pool_get_array(hydroPool, 'waterVelocityCellY', waterVelocityCellY) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) allocate(cellJunk(nCells+1)) call mpas_reconstruct(meshPool, waterVelocity, waterVelocityCellX, waterVelocityCellY, & cellJunk, cellJunk, cellJunk) @@ -476,7 +476,7 @@ subroutine li_SGH_solve(domain, err) call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) divergence(:) = 0.0_RKIND ! zero div before starting ! loop over locally owned cells @@ -551,7 +551,7 @@ subroutine li_SGH_solve(domain, err) call mpas_pool_get_array(hydroPool, 'divergence', divergence) call mpas_pool_get_array(hydroPool, 'divergenceChannel', divergenceChannel) call mpas_pool_get_array(hydroPool, 'channelAreaChangeCell', channelAreaChangeCell) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) waterThicknessOld = waterThickness waterThickness = waterThicknessOld + deltatSGH * ( (basalMeltInput + externalWaterInput) / rho_water - divergence & @@ -755,7 +755,7 @@ subroutine calc_edge_quantities(block, err) call mpas_pool_get_array(hydroPool, 'hydropotentialBaseSlopeNormal', hydropotentialBaseSlopeNormal) call mpas_pool_get_array(hydroPool, 'waterPressureSlopeNormal', waterPressureSlopeNormal) call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell) call mpas_pool_get_array(hydroPool, 'hydropotentialBaseSlopeTangent', hydropotentialBaseSlopeTangent) @@ -1217,7 +1217,7 @@ subroutine calc_pressure(block, err) call mpas_pool_get_array(velocityPool, 'basalSpeed', basalSpeed) call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA) call mpas_pool_get_array(geometryPool, 'thickness', thickness) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) openingRate(:) = bedRough * basalSpeed(:) * (bedRoughMax - waterThickness(:)) @@ -1338,7 +1338,7 @@ subroutine calc_pressure_diag_vars(block, err) call mpas_pool_get_array(hydroPool, 'hydropotentialBase', hydropotentialBase) call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness) call mpas_pool_get_array(hydroPool, 'hydropotential', hydropotential) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) effectivePressure = rhoi * gravity * thickness - waterPressure ! < this should evalute to 0 for floating ice if Pw set correctly there. diff --git a/src/core_landice/mode_forward/mpas_li_thermal.F b/src/core_landice/mode_forward/mpas_li_thermal.F index 383ba22dae..6a4807a1dd 100644 --- a/src/core_landice/mode_forward/mpas_li_thermal.F +++ b/src/core_landice/mode_forward/mpas_li_thermal.F @@ -831,7 +831,7 @@ subroutine li_thermal_solver(domain, err) ! get fields from the geometry pool call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'groundedBasalMassBal', groundedBasalMassBal) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal) call mpas_pool_get_array(geometryPool, 'basalWaterThickness', basalWaterThickness) diff --git a/src/core_landice/mode_forward/mpas_li_velocity.F b/src/core_landice/mode_forward/mpas_li_velocity.F index 4f5f6a2820..59a7ed5e2e 100644 --- a/src/core_landice/mode_forward/mpas_li_velocity.F +++ b/src/core_landice/mode_forward/mpas_li_velocity.F @@ -268,8 +268,9 @@ subroutine li_velocity_solve(domain, solveVelo, err) integer, pointer :: nEdges integer, pointer :: nVertInterfaces integer, pointer :: nCells - integer, dimension(:), pointer :: edgeMask, cellMask, vertexMask, vertexMaskOld + integer, dimension(:), pointer :: edgeMask, cellMask, cellMaskOld, vertexMask, vertexMaskOld integer, dimension(:,:), pointer :: dirichletVelocityMaskOld, dirichletVelocityMaskNew + integer, dimension(:), pointer :: floatingEdgesOld, floatingEdgesNew real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, normalVelocityInitial real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ, & uReconstructZonal, uReconstructMeridional @@ -287,7 +288,9 @@ subroutine li_velocity_solve(domain, solveVelo, err) real (kind=RKIND), dimension(:), pointer :: upperSurface integer, dimension(:), pointer :: indexToEdgeID integer, pointer :: anyDynamicVertexMaskChanged + integer, pointer :: cellIceMaskChanged integer, pointer :: dirichletMaskChanged + integer, pointer :: floatingEdgesChanged ! truly local variables integer :: cell1, cell2 integer :: cell3, cell4, thisCell @@ -300,7 +303,9 @@ subroutine li_velocity_solve(domain, solveVelo, err) real (kind=RKIND) :: maxThicknessOnProc, maxThicknessAllProcs real (kind=RKIND) :: xVelEdge, yVelEdge integer :: blockDynamicVertexMaskChanged, procDynamicVertexMaskChanged + integer :: blockCellIceMaskChanged, procCellIceMaskChanged integer :: blockDirichletMaskChanged, procDirichletMaskChanged + integer :: blockFloatingEdgesChanged, procFloatingEdgesChanged call mpas_timer_start("velocity solve") @@ -385,7 +390,9 @@ subroutine li_velocity_solve(domain, solveVelo, err) ! Don't bother checking this with SIA because it requires an extra global reduce maxThicknessOnProc = 0.0_RKIND ! initialize to procDynamicVertexMaskChanged = 0 + procCellIceMaskChanged = 0 procDirichletMaskChanged = 0 + procFloatingEdgesChanged = 0 block => domain % blocklist do while (associated(block)) @@ -397,7 +404,7 @@ subroutine li_velocity_solve(domain, solveVelo, err) maxThicknessOnProc = max(maxThicknessOnProc, maxval(thickness)) ! The interface expects an array where 1's are floating edges and 0's are non-floating edges. - call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges) + call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) floatingEdges = li_mask_is_floating_ice_int(edgeMask) call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel=1) @@ -407,7 +414,7 @@ subroutine li_velocity_solve(domain, solveVelo, err) ! TODO: There may be some aspects of the mask that are ok change for external dycores, ! but for now just check the whole thing. call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMaskOld, timeLevel=2) - if ( sum(li_mask_is_dynamic_ice_int(vertexMask) - li_mask_is_dynamic_ice_int(vertexMaskOld)) /= 0 ) then + if ( sum(abs(li_mask_is_dynamic_ice_int(vertexMask) - li_mask_is_dynamic_ice_int(vertexMaskOld))) /= 0 ) then blockDynamicVertexMaskChanged = 1 else blockDynamicVertexMaskChanged = 0 @@ -416,11 +423,26 @@ subroutine li_velocity_solve(domain, solveVelo, err) ! Determine if any blocks on this processor had a change to the vertex mask procDynamicVertexMaskChanged = max(procDynamicVertexMaskChanged, blockDynamicVertexMaskChanged) !print *,'procVertexMaskChanged', procVertexMaskChanged + + ! Determine if the cell mask changed during this time step for this block (needed for external dycores) + ! TODO: There may be some aspects of the mask that are ok change for external dycores, + ! but for now just check the whole thing. + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel=1) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMaskOld, timeLevel=2) + if (( sum(abs(li_mask_is_dynamic_ice_int(cellMask) - li_mask_is_dynamic_ice_int(cellMaskOld))) /= 0 ) .or. & + ( sum(abs(li_mask_is_ice_int(cellMask) - li_mask_is_ice_int(cellMaskOld))) /= 0 )) then + blockCellIceMaskChanged = 1 + else + blockCellIceMaskChanged = 0 + endif + + ! Determine if any blocks on this processor had a change to the cell mask + procCellIceMaskChanged = max(procCellIceMaskChanged, blockCellIceMaskChanged) - ! Also check to see if the Dirichlet b.c. mask has changed + ! Check to see if the Dirichlet b.c. mask has changed call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMaskOld, timeLevel=2) call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMaskNew, timeLevel=1) - if ( sum(dirichletVelocityMaskNew - dirichletVelocityMaskOld) /= 0 ) then + if ( sum(abs(dirichletVelocityMaskNew - dirichletVelocityMaskOld)) /= 0 ) then blockDirichletMaskChanged = 1 else blockDirichletMaskChanged = 0 @@ -428,6 +450,17 @@ subroutine li_velocity_solve(domain, solveVelo, err) ! Determine if any blocks on this processor had a change to the vertex mask procDirichletMaskChanged = max(procDirichletMaskChanged, blockDirichletMaskChanged) + ! Check to see if the Floating Edges Ids have changed + call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdgesOld, timeLevel=2) + call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdgesNew, timeLevel=1) + if ( sum(abs(floatingEdgesNew - floatingEdgesOld)) /= 0 ) then + blockFloatingEdgesChanged = 1 + else + blockFloatingEdgesChanged = 0 + endif + ! Determine if any blocks on this processor had a change to the vertex mask + procFloatingEdgesChanged = max(procFloatingEdgesChanged, blockFloatingEdgesChanged) + ! Set beta for solver to use ! (this could potentially be applied to an SIA solver, so it is calculated in this module) call calculate_beta(block, err_tmp) @@ -451,10 +484,16 @@ subroutine li_velocity_solve(domain, solveVelo, err) call mpas_pool_get_array(velocityPool, 'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged) call mpas_dmpar_max_int(domain % dminfo, procDynamicVertexMaskChanged, anyDynamicVertexMaskChanged) !print *,'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged + ! Do the same for the cell mask + call mpas_pool_get_array(velocityPool, 'cellIceMaskChanged', cellIceMaskChanged) + call mpas_dmpar_max_int(domain % dminfo, procCellIceMaskChanged, cellIceMaskChanged) ! Do the same for the Dirichlet b.c. mask call mpas_pool_get_array(velocityPool, 'dirichletMaskChanged', dirichletMaskChanged) call mpas_dmpar_max_int(domain % dminfo, procDirichletMaskChanged, dirichletMaskChanged) - !print *,'dirichletMaskChanged', dirichletMaskChanged + ! Do the same for the Floating edges + call mpas_pool_get_array(velocityPool, 'floatingEdgesChanged', floatingEdgesChanged) + call mpas_dmpar_max_int(domain % dminfo, procFloatingEdgesChanged, floatingEdgesChanged) + !print *,'floatingEdgesChanged', floatingEdgesChanged endif @@ -478,7 +517,7 @@ subroutine li_velocity_solve(domain, solveVelo, err) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(meshPool, 'cellsOnVertex', cellsOnVertex) call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface) @@ -823,7 +862,7 @@ subroutine calculate_beta(block, err) call mpas_pool_get_array(velocityPool, 'betaSolve', betaSolve) call mpas_pool_get_array(velocityPool, 'beta', beta) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) betaSolve = beta @@ -999,7 +1038,7 @@ subroutine remove_fast_ice(meshPool, uReconstructX, uReconstructY, normalVelocit call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface) diff --git a/src/core_landice/mode_forward/mpas_li_velocity_external.F b/src/core_landice/mode_forward/mpas_li_velocity_external.F index a214560efe..0235126917 100644 --- a/src/core_landice/mode_forward/mpas_li_velocity_external.F +++ b/src/core_landice/mode_forward/mpas_li_velocity_external.F @@ -452,7 +452,9 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro logical, pointer :: config_output_external_velocity_solver_data real (kind=RKIND), pointer :: config_ice_density integer, pointer :: anyDynamicVertexMaskChanged + integer, pointer :: cellIceMaskChanged integer, pointer :: dirichletMaskChanged + integer, pointer :: floatingEdgesChanged integer, pointer :: nEdges integer, pointer :: nVertLevels integer, pointer :: timestepNumber @@ -486,7 +488,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface) call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface) call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel = 1) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal) @@ -501,9 +503,11 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call mpas_pool_get_array(velocityPool, 'uReconstructZ', uReconstructZ) call mpas_pool_get_array(velocityPool, 'betaSolve', betaSolve) call mpas_pool_get_array(velocityPool, 'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged) + call mpas_pool_get_array(velocityPool, 'cellIceMaskChanged', cellIceMaskChanged) call mpas_pool_get_array(velocityPool, 'dirichletMaskChanged', dirichletMaskChanged) call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1) - call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges) + call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges, timeLevel = 1) + call mpas_pool_get_array(velocityPool, 'floatingEdgesChanged', floatingEdgesChanged) call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor) ! Hydro variables @@ -521,7 +525,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro ! Note these functions will always be called on the first solve because we ! initialize vertexMask to garbage which sets anyDynamicVertexMaskChanged to 1. if ((anyDynamicVertexMaskChanged == 1) .or. (config_always_compute_fem_grid) .or. & - (dirichletMaskChanged == 1) ) then + (dirichletMaskChanged == 1) .or. (floatingEdgesChanged == 1) .or. (cellIceMaskChanged == 1)) then call mpas_log_write("Generating new external velocity solver FEM grid.", flushNow=.true.) call generate_fem_grid(config_velocity_solver, vertexMask, cellMask, dirichletVelocityMask, & floatingEdges, layerThicknessFractions, lowerSurface, thickness, err) @@ -809,13 +813,13 @@ subroutine li_velocity_external_write_albany_mesh(domain) call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal) call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel = 1) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) ! Velocity variables call mpas_pool_get_array(velocityPool, 'beta', beta) - call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges) + call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges, timeLevel = 1) call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1) call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor) From 297eb8825548081f874f09588f226ccab62f204a Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Fri, 8 Nov 2019 14:25:05 -0700 Subject: [PATCH 3/9] Remove old code for LIFEV from the Interface_velocty_solver.*pp. Cleaning. --- .../Interface_velocity_solver.cpp | 291 +----------------- .../Interface_velocity_solver.hpp | 53 ---- 2 files changed, 5 insertions(+), 339 deletions(-) diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.cpp b/src/core_landice/mode_forward/Interface_velocity_solver.cpp index ca3e511c54..58b13e1d9a 100644 --- a/src/core_landice/mode_forward/Interface_velocity_solver.cpp +++ b/src/core_landice/mode_forward/Interface_velocity_solver.cpp @@ -15,9 +15,6 @@ distributed with this code, or at http://mpas-dev.github.com/license.html #include #include #include "Interface_velocity_solver.hpp" -//#include -//#include -//#include // =================================================== //! Namespaces @@ -131,11 +128,6 @@ void velocity_solver_export_2d_data(double const* lowerSurface_F, double const* thickness_F, double const* beta_F) { if (isDomainEmpty) return; -#ifdef LIFEV - import2DFields(lowerSurface_F, thickness_F, beta_F, minThickneess); - velocity_solver_export_2d_data__(reducedComm, elevationData, thicknessData, - betaData, indexToVertexID); -#endif } void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F, @@ -210,23 +202,6 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F, } void velocity_solver_init_l1l2(double const* levelsRatio_F) { -#ifdef LIFEV - velocityOnVertices.resize(2 * nVertices * (nLayers + 1), 0.); - velocityOnCells.resize(2 * nCells_F * (nLayers + 1), 0.); - - if (isDomainEmpty) - return; - - layersRatio.resize(nLayers); - // !!Indexing of layers is reversed - for (int i = 0; i < nLayers; i++) - layersRatio[i] = levelsRatio_F[nLayers - 1 - i]; - //std::copy(levelsRatio_F, levelsRatio_F+nLayers, layersRatio.begin()); - mapCellsToVertices(velocityOnCells, velocityOnVertices, 2, nLayers, Ordering); - - velocity_solver_init_l1l2__(layersRatio, velocityOnVertices, initialize_velocity); - initialize_velocity = false; -#endif } @@ -236,82 +211,9 @@ void velocity_solver_solve_l1l2(double const* lowerSurface_F, double const* thickness_F, double const* beta_F, double const* temperature_F, double* const dirichletVelocityXValue, double* const dirichletVelocitYValue, double* u_normal_F, double* xVelocityOnCell, double* yVelocityOnCell) { - -#ifdef LIFEV - - std::fill(u_normal_F, u_normal_F + nEdges_F * (nLayers+1), 0.); - - double localSum(0), sum(0); - - for (int i = 0; i < nCellsSolve_F; i++) { - localSum = std::max(localSum, - std::fabs(thickness_F[i] - thicknessOnCells[i])); - } - - MPI_Allreduce(&localSum, &sum, 1, MPI_DOUBLE, MPI_MAX, comm); - - std::cout << "Thickness change: " << sum << std::endl; - std::copy(thickness_F, thickness_F + nCellsSolve_F, &thicknessOnCells[0]); - - if (!isDomainEmpty) { - std::vector temperatureData(nLayers * nVertices); - - import2DFields(lowerSurface_F, thickness_F, beta_F, minThickness); - - for (int index = 0; index < nVertices; index++) { - int iCell = vertexToFCell[index]; - for (int il = 0; il < nLayers; il++) { - temperatureData[index + il * nVertices] = temperature_F[iCell * nLayers - + (nLayers - il - 1)]; - } - } - - - velocity_solver_solve_l1l2__(elevationData, thicknessData, betaData, - temperatureData, indexToVertexID, velocityOnVertices); - } - - - mapVerticesToCells (velocityOnVertices, &velocityOnCells[0], 2, nLayers, Ordering); - - //computing x, yVelocityOnCell - int sizeVelOnCell = nCells_F * (nLayers + 1); - for(int iCell=0; iCell regulThk(thicknessData); - for (int index = 0; index < nVertices; index++) - regulThk[index] = std::max(1e-4, thicknessData[index]); - - std::vector mpasIndexToVertexID(nVertices); - for (int i = 0; i < nVertices; i++) { - mpasIndexToVertexID[i] = indexToCellID_F[vertexToFCell[i]]; - } -#ifdef LIFEV - velocity_solver_export_l1l2_velocity__(layersRatio, elevationData, regulThk, mpasIndexToVertexID, reducedComm); -#endif -} +void velocity_solver_export_l1l2_velocity() {} void velocity_solver_init_fo(double const *levelsRatio_F) { @@ -325,15 +227,9 @@ void velocity_solver_init_fo(double const *levelsRatio_F) { // !!Indexing of layers is reversed for (int i = 0; i < nLayers; i++) layersRatio[i] = levelsRatio_F[nLayers - 1 - i]; - //std::copy(levelsRatio_F, levelsRatio_F+nLayers, layersRatio.begin()); - mapCellsToVertices(velocityOnCells, velocityOnVertices, 2, nLayers, Ordering); -#ifdef LIFEV - velocity_solver_init_fo__(layersRatio, velocityOnVertices, indexToVertexID, initialize_velocity); -#endif - // iceProblemPtr->initializeSolverFO(layersRatio, velocityOnVertices, thicknessData, elevationData, indexToVertexID, initialize_velocity); initialize_velocity = false; } @@ -368,22 +264,6 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower if (!isDomainEmpty) { -#ifdef LIFEV - double localSum(0), sum(0); - - for (int i = 0; i < nCellsSolve_F; i++) { - localSum = std::max(localSum, - std::fabs(thickness_F[i] - thicknessOnCells[i])); - } - - MPI_Allreduce(&localSum, &sum, 1, MPI_DOUBLE, MPI_MAX, comm); - - std::cout << "Thickness change: " << sum << std::endl; - std::copy(thickness_F, thickness_F + nCellsSolve_F, &thicknessOnCells[0]); -#endif - - - std::map bdExtensionMap; import2DFields(bdExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, temperature_F, smb_F, minThickness); @@ -444,65 +324,6 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower allToAll(u_normal_F, sendEdgesList_F, recvEdgesList_F, nLayers+1); first_time_step = false; - -#ifdef LIFEV - - std::vector edgesProcId(nEdges_F), trianglesProcIds(nVertices_F); - getProcIds(edgesProcId, recvEdgesList_F); - getProcIds(trianglesProcIds, recvVerticesList_F); - - int localSumInt(0), sumInt(0); - - for (int i = 0; i < nEdges; i++) { - for (int il = 0; il < 1; il++) { - if (std::fabs( - velOnEdges[i * (nLayers+1) + il] - - u_normal_F[edgeToFEdge[i] * nLayers + il]) > 1e-9) - // if(edgeToFEdge[i]>nEdgesSolve_F) - { - localSumInt++; - int edge = edgeToFEdge[i]; - int gEdge = indexToEdgeID[i]; - ID fVertex0 = verticesOnEdge_F[2 * edge] - 1; - ID fVertex1 = verticesOnEdge_F[2 * edge + 1] - 1; - ID triaId0 = fVertexToTriangleID[fVertex0]; - ID triaId1 = fVertexToTriangleID[fVertex1]; - ID procTria0 = trianglesProcIds[fVertex0]; - ID procTria1 = trianglesProcIds[fVertex1]; - std::cout << "vs( " << velOnEdges[i * (nLayers+1) + il] << ", " - << u_normal_F[edgeToFEdge[i] * nLayers + il] << ") "; - std::cout << "edge: " << edge << ", gEdge: " << gEdge << ", on proc: " - << edgesProcId[edgeToFEdge[i]]; - if (triaId0 != NotAnId) { - std::cout << ". first tria0: " << triaId0 << " on proc: " - << procTria0; - } - if (triaId1 != NotAnId) { - std::cout << ".. second tria0:" << std::endl; - } - if ((triaId0 == NotAnId) || (triaId1 == NotAnId)) { - std::cout << ". and to Tria: " << triaId1 << " on proc: " << procTria1 - << std::endl; - } - - } - - //localSum = std::max(localSum, std::fabs(velOnEdges[i*nLayers+il] - u_normal_F[edgeToFEdge[i]*nLayers+il])); - } - } - - MPI_Allreduce(&localSumInt, &sumInt, 1, MPI_INT, MPI_SUM, comm); - - int localNum(sendEdgesListReversed.size()), num(0); - - MPI_Allreduce(&localNum, &num, 1, MPI_INT, MPI_SUM, comm); - - std::cout << "Edges change: " << sumInt << " " << num << std::endl; - -#endif - - - } @@ -526,7 +347,7 @@ void velocity_solver_finalize() { /*duality: * - * mpas(F) | lifev + * mpas(F) | Albany LandIce (C++) * ---------|--------- * cell | vertex * vertex | triangle @@ -722,8 +543,6 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel fCellsToReceive.clear(); - // if(! isDomainEmpty) - // { fCellsToReceive.reserve(nCells_F - nCellsSolve_F); for (int i = 0; i < nCells_F; i++) { bool isMine = i < nCellsSolve_F; @@ -749,7 +568,6 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel vertexToFCell.push_back(i); } } - // } //Compute the global number of vertices, and the localOffset on the local processor, such that a globalID = localOffset + index computeLocalOffset(fCellsToSend.size(), localOffset, nGlobalVertices); @@ -887,43 +705,7 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel if (isDomainEmpty) return; -#ifdef LIFEV - std::vector verticesCoords(3 * nVertices); - - for (int index = 0; index < nVertices; index++) { - int iCell = vertexToFCell[index]; - verticesCoords[index * 3] = xCell_F[iCell] / unit_length; - verticesCoords[index * 3 + 1] = yCell_F[iCell] / unit_length; - verticesCoords[index * 3 + 2] = zCell_F[iCell] / unit_length; - } - - velocity_solver_compute_2d_grid__(nGlobalTriangles, - nGlobalVertices, nGlobalEdges, indexToVertexID, verticesCoords, - isVertexBoundary, verticesOnTria,isBoundaryEdge, trianglesOnEdge, - trianglesPositionsOnEdge, verticesOnEdge, indexToEdgeID, - indexToTriangleID, procOnInterfaceEdge ); -#else - velocity_solver_compute_2d_grid__(reducedComm); -#endif - - /* - - //initialize the mesh - iceProblemPtr->mesh2DPtr.reset (new RegionMesh() ); - - //construct the mesh nodes - constructNodes ( * (iceProblemPtr->mesh2DPtr), indexToVertexID, verticesCoords, isVertexBoundary, nGlobalVertices, 3); - - //construct the mesh elements - constructElements ( * (iceProblemPtr->mesh2DPtr), indexToTriangleID, verticesOnTria, nGlobalTriangles); - - //construct the mesh facets - constructFacets ( * (iceProblemPtr->mesh2DPtr), isBoundaryEdge, trianglesOnEdge, trianglesPositionsOnEdge, verticesOnEdge, indexToEdgeID, procOnInterfaceEdge, nGlobalEdges, 3); - - Switch sw; - std::vector elSign; - checkVolumes ( * (iceProblemPtr->mesh2DPtr), elSign, sw ); - */ + velocity_solver_compute_2d_grid__(reducedComm); } void velocity_solver_extrude_3d_grid(double const* levelsRatio_F, @@ -1124,7 +906,7 @@ void get_prism_velocity_on_FEdges(double * uNormal, //edge is on boundary and wasn't found by previous two cases //For boundary edges one of the two triangles sharing the edge won't be part of the velocity solver's mesh and the dynamic_ice_bit_value will be 0. - //Compute iCells containinig the vertices of the triangle that is part of the mesh and bcoords the corresponding barycentric coordinates + //Compute iCells containing the vertices of the triangle that is part of the mesh and bcoords the corresponding barycentric coordinates if(verticesMask_F[fVertex0] & dynamic_ice_bit_value) { belongToTria(e_mid, t0, bcoords); for (int j = 0; j < 3; j++) iCells[j] = cellsOnVertex_F[3 * fVertex0 + j] - 1; @@ -1222,7 +1004,6 @@ void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, getProcIds(cellsProcId, recvCellsList_F); getProcIds(trianglesProcIds, recvVerticesList_F); - //std::cout << "SendList " ; for (int i = 0; i < nVertices; i++) { int iCell = vertexToFCell[i]; if (iCell < nCellsSolve_F) @@ -1239,11 +1020,8 @@ void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, if (!belongToTriaOnSameProc) { sendMap[cellsProcId[iCell]].insert( std::make_pair(fCellToVertexID[iCell], iCell)); - // std::cout<< "(" << cellsProcId[iCell] << "," << iCell << ") "; } - } - //std::cout < >::const_iterator it = sendMap.begin(); it != sendMap.end(); it++) { @@ -1256,7 +1034,6 @@ void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, exchange(it->first, &sendVec[0], &sendVec[0] + sendVec.size())); } - //std::cout << "ReceiveList " ; for (UInt i = 0; i < fCellsToReceive.size(); i++) { int iCell = fCellsToReceive[i]; int nEdg = nEdgesOnCells_F[iCell]; @@ -1266,11 +1043,9 @@ void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, if (triaId != NotAnId) { receiveMap[trianglesProcIds[fVertex]].insert( std::make_pair(fCellToVertexID[iCell], iCell)); - // std::cout<< "(" << trianglesProcIds[fVertex] << "," << iCell << ") "; } } } - //std::cout < >::const_iterator it = receiveMap.begin(); it != receiveMap.end(); it++) { @@ -1297,7 +1072,6 @@ void createReverseEdgesExchangeLists(exchangeList_Type& sendListReverse_F, getProcIds(edgesProcId, recvEdgesList_F); getProcIds(trianglesProcIds, recvVerticesList_F); - //std::cout << "EdgesSendList " ; for (int i = 0; i < nEdges; i++) { int iEdge = edgeToFEdge[i]; if (iEdge < nEdgesSolve_F) @@ -1313,10 +1087,8 @@ void createReverseEdgesExchangeLists(exchangeList_Type& sendListReverse_F, if (!belongToTriaOnSameProc) { sendMap[edgesProcId[iEdge]].insert( std::make_pair(fEdgeToEdgeID[iEdge], iEdge)); - //std::cout<< "(" << edgesProcId[iEdge] << "," << iEdge << ") "; } } - //std::cout < >::const_iterator it = sendMap.begin(); it != sendMap.end(); it++) { @@ -1329,7 +1101,6 @@ void createReverseEdgesExchangeLists(exchangeList_Type& sendListReverse_F, exchange(it->first, &sendVec[0], &sendVec[0] + sendVec.size())); } - //std::cout << "EdgesReceiveList " ; for (UInt i = 0; i < edgesToReceive.size(); i++) { int iEdge = edgesToReceive[i]; for (int j = 0; j < 2; j++) { @@ -1338,11 +1109,9 @@ void createReverseEdgesExchangeLists(exchangeList_Type& sendListReverse_F, if (triaId != NotAnId) { receiveMap[trianglesProcIds[fVertex]].insert( std::make_pair(fEdgeToEdgeID[iEdge], iEdge)); - // std::cout<< "(" << trianglesProcIds[fVertex] << "," << iEdge << ") "; } } } - // std::cout < >::const_iterator it = receiveMap.begin(); it != receiveMap.end(); it++) { @@ -1417,25 +1186,6 @@ double signedTriangleAreaOnSphere(const double* x, const double* y, > 0) ? area : -area; } -//TO BE FIXED, Access To verticesOnCell_F is not correct -void extendMaskByOneLayer(int const* verticesMask_F, - std::vector& extendedFVerticesMask) { - extendedFVerticesMask.resize(nVertices_F); - extendedFVerticesMask.assign(&verticesMask_F[0], - &verticesMask_F[0] + nVertices_F); - for (int i = 0; i < nCells_F; i++) { - bool belongsToMarkedTriangle = false; - int nEdg = nEdgesOnCells_F[i]; - for (UInt k = 0; k < nEdg && !belongsToMarkedTriangle; k++) - belongsToMarkedTriangle = belongsToMarkedTriangle - || verticesMask_F[verticesOnCell_F[maxNEdgesOnCell_F * i + k] - 1]; - if (belongsToMarkedTriangle) - for (UInt k = 0; k < nEdg; k++) { - ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F * i + k] - 1); - extendedFVerticesMask[fVertex] = !isGhostTriangle(fVertex); - } - } -} void import2DFields(std::map bdExtensionMap, double const* bedTopography_F, double const * lowerSurface_F, double const * thickness_F, double const * beta_F, double const* stiffnessFactor_F, double const* effecPress_F, @@ -1482,7 +1232,7 @@ void import2DFields(std::map bdExtensionMap, double const* bedTopograp for (int iVertex = 0; iVertex < 3; iVertex++) { int v = verticesOnTria[iVertex + 3 * index]; int iCell = vertexToFCell[v]; - //compute temperature by averaging tmeperature values of triangles vertices where ice is present + //compute temperature by averaging temperature values of triangles vertices where ice is present if (cellsMask_F[iCell] & ice_present_bit_value) { temperature += temperature_F[iCell * nLayers + ilReversed]; nPoints++; @@ -1527,10 +1277,8 @@ void import2DFields(std::map bdExtensionMap, double const* bedTopograp int c0 = cellsOnEdge_F[2 * fEdge] - 1; int c1 = cellsOnEdge_F[2 * fEdge + 1] - 1; c = (fCellToVertex[c0] == iV) ? c1 : c0; - //if(!(cellsMask_F[c] & ice_present_bit_value)) continue; if((cellsMask_F[c] & dynamic_ice_bit_value)) { double elev = thickness_F[c] + lowerSurface_F[c]; // - 1e-8*std::sqrt(pow(xCell_F[c0],2)+std::pow(yCell_F[c0],2)); - //std::cout << " elev="< bdExtensionMap, double const* bedTopograp // Otherwise, it will have a surface elevation below sea level, // which is unphysical and can cause large slopes and other issues. if (bedTopographyData[iV] < 0.0) { - //std::cout<<"non-floating boundary below sea level at iV="< mpasIndexToVertexID(nVertices); -// for (int i = 0; i < nVertices; i++) { -// mpasIndexToVertexID[i] = indexToCellID_F[vertexToFCell[i]]; -// } -// -// mapVerticesToCells(velocityOnVertices, &velocityOnCells[0], 2, nLayers, -// Ordering); - - // Write out ASCII format std::cout << "Writing mesh to albany.msh." << std::endl; @@ -2185,7 +1905,6 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex) iVerticesBoundaryEdge ++; } } - //std::cout<<"final count: "< #include #include #include @@ -45,11 +44,6 @@ distributed with this code, or at http://mpas-dev.github.com/license.html #define write_ascii_mesh write_ascii_mesh_ #endif -//#include -//#include - -//#include - struct exchange { const int procID; const std::vector vec; @@ -126,8 +120,6 @@ void velocity_solver_export_l1l2_velocity(); void velocity_solver_export_fo_velocity(); -//void velocity_solver_estimate_SS_SMB (const double* u_normal_F, double* sfcMassBal); - void interface_init_log(); void interface_redirect_stdout(int const* iTimestep); @@ -151,21 +143,6 @@ void write_ascii_mesh(int const* indexToCellID_F, extern int velocity_solver_init_mpi__(int* fComm); extern void velocity_solver_finalize__(); -#ifdef LIFEV -extern void velocity_solver_init_l1l2__(const std::vector& layersRatio, const std::vector& velocityOnVertices, bool initialize_velocity); - -extern void velocity_solver_solve_l1l2__(const std::vector& elevationData, - const std::vector& thicknessData, const std::vector& betaData, - const std::vector& temperatureData, const std::vector& indexToVertexID, - std::vector& velocityOnVertices); - -extern void velocity_solver_init_fo__(const std::vector& layersRatio, const std::vector& velocityOnVertices, const std::vector& indexToVertexID, bool initialize_velocity); - -extern void velocity_solver_export_l1l2_velocity__(const std::vector& layersRatio, const std::vector& elevationData, const std::vector& regulThk, - const std::vector& mpasIndexToVertexID, MPI_Comm reducedComm); - -#endif - extern void velocity_solver_set_physical_parameters__(double const& gravity, double const& ice_density, double const& ocean_density, double const& sea_level, double const& flowParamA, double const& flowLawExponent, double const& dynamic_thickness, bool const& useGLP, double const& clausiusClapeyronCoeff); @@ -188,25 +165,7 @@ extern void velocity_solver_solve_fo__(int nLayers, int nGlobalVertices, int& error, const double& deltat = 0.0); - -#ifdef LIFEV -extern void velocity_solver_compute_2d_grid__(int nGlobalTriangles, - int nGlobalVertices, int nGlobalEdges, - const std::vector& indexToVertexID, - const std::vector& verticesCoords, - const std::vector& isVertexBoundary, - const std::vector& verticesOnTria, - const std::vector& isBoundaryEdge, - const std::vector& trianglesOnEdge, - const std::vector& trianglesPositionsOnEdge, - const std::vector& verticesOnEdge, - const std::vector& indexToEdgeID, - const std::vector& indexToTriangleID, - const std::vector < std::pair >& procOnInterfaceEdge); - -#else extern void velocity_solver_compute_2d_grid__(MPI_Comm); -#endif extern void velocity_solver_export_2d_data__(MPI_Comm reducedComm, @@ -231,17 +190,8 @@ extern void velocity_solver_extrude_3d_grid__(int nLayers, int nGlobalTriangles, const std::vector& dirichletNodes, const std::vector&floatingEdges); -//extern void velocity_solver_export_l1l2_velocity__(); - extern void velocity_solver_export_fo_velocity__(MPI_Comm reducedComm); - -#ifdef LIFEV -extern int velocity_solver_initialize_iceProblem__(bool keep_proc, MPI_Comm reducedComm); -#endif - -//extern void velocity_solver_estimate_SS_SMB__ (const double* u_normal_F, double* sfcMassBal); - exchangeList_Type unpackMpiArray(int const* array); bool isGhostTriangle(int i, double relTol = 1e-1); @@ -271,9 +221,6 @@ void write_ascii_mesh_field_int(std::vector fieldData, std::string filename std::vector extendMaskByOneLayer(int const* verticesMask_F); -void extendMaskByOneLayer(int const* verticesMask_F, - std::vector& extendedFVerticesMask); - void importP0Temperature(); void exportDissipationHeat(double * dissipationHeat_F); From fb222d444ed6534be7c627168408d5ac08cbdacf Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Fri, 8 Nov 2019 15:17:12 -0700 Subject: [PATCH 4/9] Modify write_ascii_mesh so that it can mesh can be printed in parallel. --- .../mode_forward/Interface_velocity_solver.cpp | 15 ++++++++++----- .../mode_forward/mpas_li_velocity_external.F | 4 ---- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.cpp b/src/core_landice/mode_forward/Interface_velocity_solver.cpp index 58b13e1d9a..9638a383df 100644 --- a/src/core_landice/mode_forward/Interface_velocity_solver.cpp +++ b/src/core_landice/mode_forward/Interface_velocity_solver.cpp @@ -1878,7 +1878,11 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex) // msh file std::ofstream outfile; outfile.precision(15); - outfile.open ("albany.msh", std::ios::out | std::ios::trunc); + std::stringstream name; + int me; + MPI_Comm_rank(comm, &me); + name << "albany.msh." << me; + outfile.open (name.str(), std::ios::out | std::ios::trunc); if (outfile.is_open()) { int nVerticesBoundaryEdge = 0; @@ -1906,7 +1910,8 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex) } } - outfile << "Triangle " << 3 << "\n"; // first line saying it is a mesh of triangles + outfile << "Format: " << 1 << "\n"; // first line stating the format (we are going to provide global ids) + outfile << "Triangle " << 3 << "\n"; // second line saying it is a mesh of triangles outfile << nVertices << " " << nTriangles << " " << nVerticesBoundaryEdge << "\n"; // second line for (int index = 0; index < nVertices; index++) { //coordinates lines @@ -1916,14 +1921,14 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex) int vertexLabel = (!isVertexBoundary[index]) ? 0 : isDirichletVertex ? 3 : 1; - outfile << xCell_F[iCell] / unit_length << " " << yCell_F[iCell] / unit_length << " " << vertexLabel << "\n" ; + outfile << indexToVertexID[index] << " " << xCell_F[iCell] / unit_length << " " << yCell_F[iCell] / unit_length << " " << vertexLabel << "\n" ; } for (int index = 0; index < nTriangles; index++) //triangles lines - outfile << verticesOnTria[0 + 3 * index] + 1 << " " << verticesOnTria[1 + 3 * index] + 1 << " " << verticesOnTria[2 + 3 * index] + 1 << " " << 1 << "\n"; // last digit can be used to specify a 'material'. Not used by Albany LandIce, so giving dummy value + outfile << indexToTriangleID[index] << " " << verticesOnTria[0 + 3 * index] + 1 << " " << verticesOnTria[1 + 3 * index] + 1 << " " << verticesOnTria[2 + 3 * index] + 1 << " " << 1 << "\n"; // last digit can be used to specify a 'material'. Not used by Albany LandIce, so giving dummy value for (int index = 0; index < nVerticesBoundaryEdge; index++) // boundary edges lines - outfile << boundaryEdges[0 + 3 * index] + 1 << " " << boundaryEdges[1 + 3 * index] + 1 << " " << boundaryEdges[2 + 3 * index] << "\n"; //last digit can be used to tell whether it's floating or not.. but let's worry about this later. + outfile << indexToEdgeID[index] << " " << boundaryEdges[0 + 3 * index] + 1 << " " << boundaryEdges[1 + 3 * index] + 1 << " " << boundaryEdges[2 + 3 * index] << "\n"; //last digit can be used to tell whether it's floating or not.. but let's worry about this later. outfile.close(); } diff --git a/src/core_landice/mode_forward/mpas_li_velocity_external.F b/src/core_landice/mode_forward/mpas_li_velocity_external.F index 0235126917..8587bd5d65 100644 --- a/src/core_landice/mode_forward/mpas_li_velocity_external.F +++ b/src/core_landice/mode_forward/mpas_li_velocity_external.F @@ -783,10 +783,6 @@ subroutine li_velocity_external_write_albany_mesh(domain) call mpas_log_write("config_velocity solver needs to be set to 'FO' for config_write_albany_ascii_mesh to work.", & MPAS_LOG_CRIT) endif - ! check nProcs - if (domain % dminfo % nProcs /= 1) then - call mpas_log_write("config_write_albany_ascii_mesh currently only works on 1 processor.", MPAS_LOG_CRIT) - endif ! check nBlocks if (domain % dminfo % total_blocks /= 1) then call mpas_log_write("config_write_albany_ascii_mesh currently only works on 1 block per processor.", MPAS_LOG_CRIT) From 69f9946d47d5ec526ada6e065d462a5eaf84dbe3 Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Fri, 8 Nov 2019 15:24:43 -0700 Subject: [PATCH 5/9] Remove functions computing the tet/prism mesh connectivity. These functions will be copied in Albany, so that the MALI can be built with shared libraries. --- .../Interface_velocity_solver.cpp | 157 +----------------- .../Interface_velocity_solver.hpp | 7 +- 2 files changed, 6 insertions(+), 158 deletions(-) diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.cpp b/src/core_landice/mode_forward/Interface_velocity_solver.cpp index 9638a383df..7d32951efb 100644 --- a/src/core_landice/mode_forward/Interface_velocity_solver.cpp +++ b/src/core_landice/mode_forward/Interface_velocity_solver.cpp @@ -741,9 +741,13 @@ void velocity_solver_extrude_3d_grid(double const* levelsRatio_F, verticesCoords[index * 3 + 2] = zCell_F[iCell] / unit_length; } + std::vector> procsSharingVertices(nVertices); + for(int i=0; i< nVertices; i++) + procsSharingVertex(i, procsSharingVertices[i]); + velocity_solver_extrude_3d_grid__(nLayers, nGlobalTriangles, nGlobalVertices, nGlobalEdges, Ordering, reducedComm, indexToVertexID, mpasIndexToVertexID, - verticesCoords, isVertexBoundary, verticesOnTria, isBoundaryEdge, + verticesCoords, isVertexBoundary, verticesOnTria, procsSharingVertices, isBoundaryEdge, trianglesOnEdge, trianglesPositionsOnEdge, verticesOnEdge, indexToEdgeID, indexToTriangleID, dirichletNodesIDs, floatingEdgesIds); } @@ -1691,157 +1695,6 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep ( bcoords[2] > -eps ); } -int prismType(long long int const* prismVertexMpasIds, int& minIndex) -{ - int PrismVerticesMap[6][6] = {{0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4}, {3, 5, 4, 0, 2, 1}, {4, 3, 5, 1, 0, 2}, {5, 4, 3, 2, 1, 0}}; - minIndex = std::min_element (prismVertexMpasIds, prismVertexMpasIds + 3) - prismVertexMpasIds; - - int v1 (prismVertexMpasIds[PrismVerticesMap[minIndex][1]]); - int v2 (prismVertexMpasIds[PrismVerticesMap[minIndex][2]]); - - return v1 > v2; -} - - void tetrasFromPrismStructured (long long int const* prismVertexMpasIds, long long int const* prismVertexGIds, long long int tetrasIdsOnPrism[][4]) - { - int PrismVerticesMap[6][6] = {{0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4}, {3, 5, 4, 0, 2, 1}, {4, 3, 5, 1, 0, 2}, {5, 4, 3, 2, 1, 0}}; - - int tetraOfPrism[2][3][4] = {{{0, 1, 2, 5}, {0, 1, 5, 4}, {0, 4, 5, 3}}, {{0, 1, 2, 4}, {0, 4, 2, 5}, {0, 4, 5, 3}}}; - - int tetraAdjacentToPrismLateralFace[2][3][2] = {{{1, 2}, {0, 1}, {0, 2}}, {{0, 2}, {0, 1}, {1, 2}}}; - int tetraFaceIdOnPrismLateralFace[2][3][2] = {{{0, 0}, {1, 1}, {2, 2}}, {{0, 0}, {1, 1}, {2, 2}}}; - int tetraAdjacentToBottomFace = 0; //does not depend on type; - int tetraAdjacentToUpperFace = 2; //does not depend on type; - int tetraFaceIdOnBottomFace = 3; //does not depend on type; - int tetraFaceIdOnUpperFace = 0; //does not depend on type; - - int minIndex; - int prismT = prismType(prismVertexMpasIds, minIndex); - - long long int reorderedPrismLIds[6]; - - for (int ii = 0; ii < 6; ii++) - { - reorderedPrismLIds[ii] = prismVertexGIds[PrismVerticesMap[minIndex][ii]]; - } - - for (int iTetra = 0; iTetra < 3; iTetra++) - for (int iVertex = 0; iVertex < 4; iVertex++) - { - tetrasIdsOnPrism[iTetra][iVertex] = reorderedPrismLIds[tetraOfPrism[prismT][iTetra][iVertex]]; - } - } - - - void computeMap() - { - int PrismVerticesMap[6][6] = {{0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4}, {3, 5, 4, 0, 2, 1}, {4, 3, 5, 1, 0, 2}, {5, 4, 3, 2, 1, 0}}; - - int tetraOfPrism[2][3][4] = {{{0, 1, 2, 5}, {0, 1, 5, 4}, {0, 4, 5, 3}}, {{0, 1, 2, 4}, {0, 4, 2, 5}, {0, 4, 5, 3}}}; - - int TetraFaces[4][3] = {{0 , 1 , 3}, {1 , 2 , 3}, {0 , 3 , 2}, {0 , 2 , 1}}; - - int PrismFaces[5][4] = {{0 , 1 , 4 , 3}, {1 , 2 , 5 , 4}, {0 , 3 , 5 , 2}, {0 , 2 , 1 , -1}, {3 , 4 , 5, -1}}; - - - for(int pType=0; pType<2; ++pType){ - std::cout<< "pType: " << pType < v2; - - for (int iTetra = 0; iTetra < 3; iTetra++) - for (int iVertex = 0; iVertex < 4; iVertex++) - { - tetrasIdsOnPrism[iTetra][iVertex] = prismVertexGIds[tetraOfPrism[prismType][iTetra][iVertex]]; - } - - // return; - - int reorderedPrismLIds[6]; - - for (int ii = 0; ii < 6; ii++) - { - reorderedPrismLIds[ii] = prismVertexGIds[PrismVerticesMap[minIndex][ii]]; - } - - for (int iTetra = 0; iTetra < 3; iTetra++) - for (int iVertex = 0; iVertex < 4; iVertex++) - { - tetrasIdsOnPrism[iTetra][iVertex] = reorderedPrismLIds[tetraOfPrism[prismType][iTetra][iVertex]]; - } - } - - - - - void setBdFacesOnPrism (const std::vector > >& prismStruct, const std::vector& prismFaceIds, std::vector& tetraPos, std::vector& facePos) - { - int numTriaFaces = prismFaceIds.size() - 2; - tetraPos.assign(numTriaFaces,-1); - facePos.assign(numTriaFaces,-1); - - - for (int iTetra (0), k (0); (iTetra < 3 && k < numTriaFaces); iTetra++) - { - bool found; - for (int jFaceLocalId = 0; jFaceLocalId < 4; jFaceLocalId++ ) - { - found = true; - for (int ip (0); ip < 3 && found; ip++) - { - int localId = prismStruct[iTetra][jFaceLocalId][ip]; - int j = 0; - found = false; - while ( (j < prismFaceIds.size()) && !found ) - { - found = (localId == prismFaceIds[j]); - j++; - } - } - if (found) - { - tetraPos[k] = iTetra; - facePos[k] = jFaceLocalId; - k += found; - break; - } - } - } - } - void procsSharingVertex(const int vertex, std::vector& procIds) { int fCell = vertexToFCell[vertex]; procIds.clear(); diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.hpp b/src/core_landice/mode_forward/Interface_velocity_solver.hpp index fbb64b4555..e7b37b5015 100644 --- a/src/core_landice/mode_forward/Interface_velocity_solver.hpp +++ b/src/core_landice/mode_forward/Interface_velocity_solver.hpp @@ -181,6 +181,7 @@ extern void velocity_solver_extrude_3d_grid__(int nLayers, int nGlobalTriangles, const std::vector& verticesCoords, const std::vector& isVertexBoundary, const std::vector& verticesOnTria, + const std::vector> procsSharingVertices, const std::vector& isBoundaryEdge, const std::vector& trianglesOnEdge, const std::vector& trianglesPositionsOnEdge, @@ -264,12 +265,6 @@ void allToAll(std::vector& field, exchangeList_Type const* sendList, void allToAll(double* field, exchangeList_Type const* sendList, exchangeList_Type const* recvList, int fieldDim = 1); -int prismType(long long int const* prismVertexMpasIds, int& minIndex); -void tetrasFromPrismStructured (long long int const* prismVertexMpasIds, long long int const* prismVertexGIds, long long int tetrasIdsOnPrism[][4]); -void computeMap(); - -void setBdFacesOnPrism (const std::vector > >& prismStruct, const std::vector& prismFaceIds, std::vector& tetraPos, std::vector& facePos); -void tetrasFromPrismStructured (int const* prismVertexMpasIds, int const* prismVertexGIds, int tetrasIdsOnPrism[][4]); void procsSharingVertex(const int vertex, std::vector& procIds); bool belongToTria(double const* x, double const* t, double bcoords[3], double eps = 1e-3); From 7dbf97bc46d2932053f790cb487c15aa9a0b3ca9 Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Fri, 8 Nov 2019 16:08:15 -0700 Subject: [PATCH 6/9] Several changes to the Interface_velocity_solver: 1. sendig now MPAS global IDs to Albany, instead of creating contiguous id (this was created for overcoming some issues with solver. It is not needed any more now. Old functionality can be restored by defining the useContiguousIds preprocessor directive.) 2. Change ownership to Albany triangles (MPAS vertices). This allows to better distribute the triangles along the interface. A triangle is set to belong to a processor P if: 1. P owns at least 2 nodes of the triangle associated to that vertex, OR 2. all the nodes of the triangle belong to three different processors, and P owns the fortran vertex and a node OR 3. the three nodes of the triangle and the fortran vertex belong to four different processors, and P owns the triangle node with the minimum ID It is possible to avoid changing ownership of triangles by undefining the preprocessor directive changeTrianglesOwnership. 3. Each node is owned by a processor that owns a triangle that contain the node This avoids situations problematic to the solver The Albany mesh will be constructed keeping the ownership of nodes and triangles as defined in MPAS Interface. Ownership of edges is not dictated by the MPAS Interface and can be arbitrarily chosen by Albany Cleaning --- .../Interface_velocity_solver.cpp | 566 +++++++++--------- .../Interface_velocity_solver.hpp | 28 +- .../mode_forward/mpas_li_velocity_external.F | 4 +- 3 files changed, 302 insertions(+), 296 deletions(-) diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.cpp b/src/core_landice/mode_forward/Interface_velocity_solver.cpp index 7d32951efb..f586626585 100644 --- a/src/core_landice/mode_forward/Interface_velocity_solver.cpp +++ b/src/core_landice/mode_forward/Interface_velocity_solver.cpp @@ -14,28 +14,30 @@ distributed with this code, or at http://mpas-dev.github.com/license.html #include #include #include +#include +#include #include "Interface_velocity_solver.hpp" // =================================================== //! Namespaces // =================================================== -//typedef std::list exchangeList_Type; +//#define useContiguousIDs +#define changeTrianglesOwnership + // ice_problem pointer int Ordering = 0; //ordering ==0 means that the mesh is extruded layerwise, whereas ordering==1 means that the mesh is extruded columnwise. MPI_Comm comm, reducedComm; bool isDomainEmpty = true; -bool initialize_velocity = true; bool first_time_step = true; int nCells_F, nEdges_F, nVertices_F; int nCellsSolve_F, nEdgesSolve_F, nVerticesSolve_F; -int nVertices, nEdges, nTriangles, nGlobalVertices, nGlobalEdges, - nGlobalTriangles; +int nVertices, nEdges, nTriangles, globalVertexStride, globalEdgeStride,globalTriangleStride; int maxNEdgesOnCell_F; int const *cellsOnEdge_F, *cellsOnVertex_F, *verticesOnCell_F, - *verticesOnEdge_F, *edgesOnCell_F, *indexToCellID_F, *nEdgesOnCells_F, + *verticesOnEdge_F, *edgesOnCell_F, *indexToCellID_F, *indexToEdgeID_F, *indexToVertexID_F, *nEdgesOnCells_F, *verticesMask_F, *cellsMask_F, *dirichletCellsMask_F, *floatingEdgesMask_F; std::vector layersRatio, levelsNormalizedThickness; int nLayers; @@ -61,10 +63,10 @@ int Interface_stdout; // the location of stdout as we use it here //void *phgGrid = 0; std::vector edgesToReceive, fCellsToReceive, indexToTriangleID, - verticesOnTria, trianglesOnEdge, trianglesPositionsOnEdge, verticesOnEdge, - trianglesProcIds, reduced_ranks; -std::vector indexToVertexID, vertexToFCell, triangleToFVertex, indexToEdgeID, edgeToFEdge, - mask, fVertexToTriangleID, fCellToVertex, floatingEdgesIds, dirichletNodesIDs; + verticesOnTria, trianglesOnEdge, verticesOnEdge, + trianglesProcIds, reduced_ranks; +std::vector indexToVertexID, vertexToFCell, vertexProcIDs, triangleToFVertex, indexToEdgeID, edgeToFEdge, + fVertexToTriangleID, fCellToVertex, floatingEdgesIds, dirichletNodesIDs; std::vector temperatureOnTetra, dissipationHeatOnTetra, velocityOnVertices, velocityOnCells, elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, temperatureData, smbData, thicknessOnCells; std::vector isVertexBoundary, isBoundaryEdge; @@ -83,8 +85,8 @@ double radius; exchangeList_Type const *sendCellsList_F = 0, *recvCellsList_F = 0; exchangeList_Type const *sendEdgesList_F = 0, *recvEdgesList_F = 0; exchangeList_Type const *sendVerticesList_F = 0, *recvVerticesList_F = 0; -exchangeList_Type sendCellsListReversed, recvCellsListReversed, - sendEdgesListReversed, recvEdgesListReversed; +exchangeList_Type sendVerticesListReversed, recvVerticesListReversed, + sendCellsListReversed, recvCellsListReversed; exchange::exchange(int _procID, int const* vec_first, int const* vec_last, int fieldDim) : @@ -103,7 +105,7 @@ int velocity_solver_init_mpi(int* fComm) { // get MPI_Comm from Fortran comm = MPI_Comm_f2c(*fComm); reducedComm = MPI_COMM_NULL; // initialize to null so we can check if set - return velocity_solver_init_mpi__(fComm); + return velocity_solver_init_mpi__(comm); } @@ -138,6 +140,8 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F, int const* _verticesOnCell_F, int const* _verticesOnEdge_F, int const* _edgesOnCell_F, int const* _nEdgesOnCells_F, int const* _indexToCellID_F, + int const* _indexToEdgeID_F, + int const* _indexToVertexID_F, double const* _xCell_F, double const* _yCell_F, double const* _zCell_F, double const* _xVertex_F, double const* _yVertex_F, double const* _zVertex_F, double const* _areaTriangle_F, @@ -160,7 +164,9 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F, verticesOnEdge_F = _verticesOnEdge_F; edgesOnCell_F = _edgesOnCell_F; nEdgesOnCells_F = _nEdgesOnCells_F; + indexToEdgeID_F = _indexToEdgeID_F; indexToCellID_F = _indexToCellID_F; + indexToVertexID_F = _indexToVertexID_F; xCell_F = _xCell_F; yCell_F = _yCell_F; zCell_F = _zCell_F; @@ -168,7 +174,6 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F, yVertex_F = _yVertex_F; zVertex_F = _zVertex_F; areaTriangle_F = _areaTriangle_F; - mask.resize(nVertices_F); thicknessOnCells.resize(nCellsSolve_F); @@ -181,9 +186,6 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F, recvVerticesList_F = new exchangeList_Type( unpackMpiArray(recvVerticesArray_F)); - trianglesProcIds.resize(nVertices_F); - getProcIds(trianglesProcIds, recvVerticesList_F); - if (radius > 10) { xCellProjected.resize(nCells_F); yCellProjected.resize(nCells_F); @@ -229,8 +231,6 @@ void velocity_solver_init_fo(double const *levelsRatio_F) { layersRatio[i] = levelsRatio_F[nLayers - 1 - i]; mapCellsToVertices(velocityOnCells, velocityOnVertices, 2, nLayers, Ordering); - - initialize_velocity = false; } void velocity_solver_solve_fo(double const* bedTopography_F, double const* lowerSurface_F, @@ -279,7 +279,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower double dt = (*deltat)/secondsInAYear; int albany_error; - velocity_solver_solve_fo__(nLayers, nGlobalVertices, nGlobalTriangles, + velocity_solver_solve_fo__(nLayers, globalVertexStride, globalTriangleStride, Ordering, first_time_step, indexToVertexID, indexToTriangleID, minBeta, regulThk, levelsNormalizedThickness, elevationData, thicknessData, betaData, bedTopographyData, smbData, @@ -308,21 +308,14 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower xVelocityOnCell[indexReversed] = velocityOnCells[index]; yVelocityOnCell[indexReversed] = velocityOnCells[index+sizeVelOnCell]; } - - if (!isDomainEmpty) get_prism_velocity_on_FEdges(u_normal_F, velocityOnCells, edgeToFEdge); - + allToAll(u_normal_F, sendEdgesList_F, recvEdgesList_F, nLayers+1); std::vector velOnEdges(nEdges * (nLayers+1)); for (int i = 0; i < nEdges; i++) { for (int il = 0; il < nLayers+1; il++) { velOnEdges[i * (nLayers+1) + il] = u_normal_F[edgeToFEdge[i] * (nLayers+1) + il]; } } - - allToAll(u_normal_F, &sendEdgesListReversed, &recvEdgesListReversed, nLayers+1); - - allToAll(u_normal_F, sendEdgesList_F, recvEdgesList_F, nLayers+1); - first_time_step = false; } @@ -357,7 +350,6 @@ void velocity_solver_finalize() { void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cellsMask_F, int const* _dirichletCellsMask_F, int const* _floatingEdgesMask_F) { int numProcs, me; - verticesMask_F = _verticesMask_F; cellsMask_F = _cellsMask_F; verticesMask_F = _verticesMask_F; @@ -370,62 +362,138 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel numProcs + 1), globalOffsetVertices(numProcs + 1), globalOffsetEdge( numProcs + 1); + + + // First, we compute the FE triangles belonging to this processor. + // If changeTrianglesOwnership is not define, the triangles belonging to this + // processor will be the subset of the triangles (MPAS vertices) owned by this proc + // that contain dynamic ice. + // If changeTrianglesOwnership is define, we rearrange the ownership of the triangles + // to improve the quality of the FE mesh and avoid corner cases (see below). + triangleToFVertex.clear(); triangleToFVertex.reserve(nVertices_F); std::vector fVertexToTriangle(nVertices_F, NotAnId); - bool changed = false; - for (int i(0); i < nVerticesSolve_F; i++) { - if ((verticesMask_F[i] & dynamic_ice_bit_value) && !isGhostTriangle(i)) { - fVertexToTriangle[i] = triangleToFVertex.size(); - triangleToFVertex.push_back(i); - } - changed = changed || (verticesMask_F[i] != mask[i]); - } - for (int i(0); i < nVertices_F; i++) - mask[i] = verticesMask_F[i]; + //vector containing proc ranks for owned and shared FE triangles + trianglesProcIds.assign(nVertices_F,NotAnId); + + //vector containing proc ranks for owned and shared MPAS cells + std::vector fCellsProcIds(nCells_F); + getProcIds(fCellsProcIds, recvCellsList_F); +#ifdef changeTrianglesOwnership + std::vector fVerticesProcIds(nVertices_F); + getProcIds(fVerticesProcIds, recvVerticesList_F); + for (int i(0); i < nVertices_F; i++) { + int cellWithMinID=nCellsSolve_F; + if ((verticesMask_F[i] & dynamic_ice_bit_value)) { + int minCellId = std::numeric_limits::max(); + int minCellIdProc(0); + + int cellProc[3]; + bool invalidCell=false; + for (int j = 0; j < 3; j++) { + int iCell = cellsOnVertex_F[3 * i + j] - 1; + if(iCell >= nCells_F) { + invalidCell = true; + break; + } + int cellID = indexToCellID_F[iCell]; + cellProc[j] = fCellsProcIds[iCell]; + if(cellID < minCellId) { + minCellId = cellID; + cellWithMinID = iCell; + minCellIdProc = cellProc[j]; + } + } + + if(invalidCell) continue; + + // the proc that owns at least 2 nodes of the triangle i. If all nodes belong to different procs, procOwns2Nodes is set to -1 + int procOwns2Nodes = ((cellProc[0] == cellProc[1]) || (cellProc[0] == cellProc[2])) ? cellProc[0] : + (cellProc[1] == cellProc[2]) ? cellProc[1] : -1; + + int vertexProc = fVerticesProcIds[i]; + bool triangleOwnsANode = (cellProc[0] == vertexProc) || (cellProc[1] == vertexProc) || (cellProc[2] == vertexProc); + + //A triangle will be owned by a proc if: + // 1. the proc owns at least 2 nodes of the triangle associated to that vertex, OR + // 2. all the nodes of the triangle belong to three different procs, and the proc owns the fortran vertex and a node OR + // 3. the three nodes of the triangle and the fortran vertex belong to four different procs, and the proc owns the node with the minimum ID - if (changed) - std::cout << "mask changed!!" << std::endl; + trianglesProcIds[i] = (procOwns2Nodes != -1) ? procOwns2Nodes : + triangleOwnsANode ? vertexProc : + minCellIdProc; - if ((me == 0) && (triangleToFVertex.size() == 0)) - for (int i(0); i < nVerticesSolve_F; i++) { - if (!isGhostTriangle(i)) { + if (trianglesProcIds[i] == me) { fVertexToTriangle[i] = triangleToFVertex.size(); triangleToFVertex.push_back(i); - break; - } + } + } + } +#else + //in this case we just set the proc ranks for owned and shared FE triangles to the be the same as MPAS owned and shared vertices + getProcIds(trianglesProcIds, recvVerticesList_F); + for (int i(0); i < nVerticesSolve_F; i++) { + if (verticesMask_F[i] & dynamic_ice_bit_value) { + fVertexToTriangle[i] = triangleToFVertex.size(); + triangleToFVertex.push_back(i); } + } +#endif nTriangles = triangleToFVertex.size(); + //Initialize the ice sheet problem with the number of FE triangles on this prov initialize_iceProblem(nTriangles); - //Compute the global number of triangles, and the localOffset on the local processor, such that a globalID = localOffset + index + //If useContiguousIDs is defined, create a contiguous list of global triangle IDs, otherwise use MPAS vertices IDs + fVertexToTriangleID.assign(nVertices_F, NotAnId); +#ifdef useContiguousIDs + // Compute the localOffset on the local processor, such that a globalID = localOffset + index int localOffset(0); - nGlobalTriangles = 0; + int nGlobalTriangles = 0; computeLocalOffset(nTriangles, localOffset, nGlobalTriangles); - - //Communicate the globalIDs, computed locally, to the other processors. - indexToTriangleID.resize(nTriangles); - - //To make local, not used - fVertexToTriangleID.assign(nVertices_F, NotAnId); - // std::vector fVertexToTriangleID(nVertices_F, NotAnId); for (int index(0); index < nTriangles; index++) fVertexToTriangleID[triangleToFVertex[index]] = index + localOffset; - +#else + for (int index(0); index < nTriangles; index++) { + fVertexToTriangleID[triangleToFVertex[index]] = indexToVertexID_F[triangleToFVertex[index]]; + } +#endif + + //Communicate the local and global triangle IDs to the other processors. +#ifdef changeTrianglesOwnership + // because we change the ownership of some triangles, we need to first communicate back to the processors that used to own those triangles + // the data of the newly owned triangles. We do this by defining "reversed" send and receive lists, communicate back using those lists, and + // then communicate "forward" using the usual send and receive lists. + // We could join these two step in one communication, but for the moment we do that separately + createReverseVerticesExchangeLists(sendVerticesListReversed, recvVerticesListReversed, trianglesProcIds, indexToVertexID_F); + allToAll(fVertexToTriangleID, &sendVerticesListReversed, &recvVerticesListReversed); + allToAll(fVertexToTriangle, &sendVerticesListReversed, &recvVerticesListReversed); + allToAll(trianglesProcIds, sendVerticesList_F, recvVerticesList_F); +#endif allToAll(fVertexToTriangleID, sendVerticesList_F, recvVerticesList_F); + allToAll(fVertexToTriangle, sendVerticesList_F, recvVerticesList_F); - for (int index(0); index < nTriangles; index++) + + //we define the vector of global triangles Ids and compute the stride between the largest and the smallest Id globally + //This will be needed by the velocity solver to create the 3D FE mesh. + indexToTriangleID.resize(nTriangles); + int maxTriangleID=std::numeric_limits::min(), minTriangleID=std::numeric_limits::max(), maxGlobalTriangleID, minGlobalTriangleID; + for (int index(0); index < nTriangles; index++) { indexToTriangleID[index] = fVertexToTriangleID[triangleToFVertex[index]]; + maxTriangleID = (indexToTriangleID[index] > maxTriangleID) ? indexToTriangleID[index] : maxTriangleID; + minTriangleID = (indexToTriangleID[index] < minTriangleID) ? indexToTriangleID[index] : minTriangleID; + } - //Compute triangle edges - std::vector fEdgeToEdge(nEdges_F), edgesToSend, trianglesProcIds( - nVertices_F); - getProcIds(trianglesProcIds, recvVerticesList_F); + MPI_Allreduce(&maxTriangleID, &maxGlobalTriangleID, 1, MPI_INT, MPI_MAX, comm); + MPI_Allreduce(&minTriangleID, &minGlobalTriangleID, 1, MPI_INT, MPI_MIN, comm); + globalTriangleStride = maxGlobalTriangleID - minGlobalTriangleID +1; - int interfaceSize(0); + // Second, we compute the FE edges belonging to the FE triangles owned by this processor. + // We first compute boundary edges, and then all the other edges. + std::vector fEdgeToEdge(nEdges_F), edgesToSend; std::vector fEdgeToEdgeID(nEdges_F, NotAnId); edgesToReceive.clear(); @@ -439,19 +507,26 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel edgesToSend.reserve(nEdgesSolve_F); isBoundaryEdge.reserve(nEdges_F); - //first, we compute boundary edges (boundary edges must be the first edges) + //we compute boundary edges (boundary edges must be the first edges) for (int i = 0; i < nEdges_F; i++) { ID fVertex1(verticesOnEdge_F[2 * i] - 1), fVertex2( verticesOnEdge_F[2 * i + 1] - 1); + + // skip the (shared) edge when the associated MPAS vertices are not valid + if((fVertex1>=nVertices_F) || (fVertex2>=nVertices_F)) + continue; + ID triaId_1 = fVertexToTriangleID[fVertex1]; ID triaId_2 = fVertexToTriangleID[fVertex2]; bool isboundary = (triaId_1 == NotAnId) || (triaId_2 == NotAnId); ID iTria1 = fVertexToTriangle[fVertex1]; ID iTria2 = fVertexToTriangle[fVertex2]; - if (iTria1 == NotAnId) + if (trianglesProcIds[fVertex1] != me) { + std::swap(fVertex1, fVertex2); std::swap(iTria1, iTria2); - bool belongsToLocalTriangle = (iTria1 != NotAnId) || (iTria2 != NotAnId); + } + bool belongsToLocalTriangle = (trianglesProcIds[fVertex1] == me); if (belongsToLocalTriangle) { if (isboundary) { @@ -460,36 +535,36 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel trianglesOnEdge.push_back(iTria1); trianglesOnEdge.push_back(iTria2); isBoundaryEdge.push_back(true); - } else - interfaceSize += (iTria2 == NotAnId); + } } } numBoundaryEdges = edgeToFEdge.size(); - //procOnInterfaceEdge contains the pairs . - std::vector < std::pair > procOnInterfaceEdge; - procOnInterfaceEdge.reserve(interfaceSize); - //then, we compute the other edges for (int i = 0; i < nEdges_F; i++) { ID fVertex1(verticesOnEdge_F[2 * i] - 1), fVertex2( verticesOnEdge_F[2 * i + 1] - 1); + + // skip the (shared) edge when the associated MPAS vertices are not valid + if((fVertex1>=nVertices_F) || (fVertex2>=nVertices_F)) + continue; + ID iTria1 = fVertexToTriangle[fVertex1]; ID iTria2 = fVertexToTriangle[fVertex2]; ID triaId_1 = fVertexToTriangleID[fVertex1]; //global Triangle ID triaId_2 = fVertexToTriangleID[fVertex2]; //global Triangle - if (iTria1 == NotAnId) { - std::swap(iTria1, iTria2); + if (trianglesProcIds[fVertex1] != me) { std::swap(fVertex1, fVertex2); + std::swap(iTria1, iTria2); } bool belongsToAnyTriangle = (triaId_1 != NotAnId) || (triaId_2 != NotAnId); bool isboundary = (triaId_1 == NotAnId) || (triaId_2 == NotAnId); - bool belongsToLocalTriangle = (iTria1 != NotAnId); + bool belongsToLocalTriangle = (trianglesProcIds[fVertex1] == me); bool isMine = i < nEdgesSolve_F; if (belongsToLocalTriangle && !isboundary) { @@ -498,9 +573,6 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel trianglesOnEdge.push_back(iTria1); trianglesOnEdge.push_back(iTria2); isBoundaryEdge.push_back(false); - if (iTria2 == NotAnId) - procOnInterfaceEdge.push_back( - std::make_pair(fEdgeToEdge[i], trianglesProcIds[fVertex2])); } if (belongsToAnyTriangle && isMine) { @@ -511,27 +583,42 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel } - //Compute the global number of edges, and the localOffset on the local processor, such that a globalID = localOffset + index + //If useContiguousIDs is defined, create a contiguous list of global edge IDs, otherwise use MPAS edges IDs +#ifdef useContiguousIDs + //Compute the localOffset on the local processor, such that a globalID = localOffset + index + int nGlobalEdges; computeLocalOffset(edgesToSend.size(), localOffset, nGlobalEdges); - - //Communicate the globalIDs, computed locally, to the other processors. for (ID index = 0; index < edgesToSend.size(); index++) fEdgeToEdgeID[edgesToSend[index]] = index + localOffset; - + //Communicate the globalIDs, computed locally, to the other processors. allToAll(fEdgeToEdgeID, sendEdgesList_F, recvEdgesList_F); +#else + for (int fEdge = 0; fEdge < nEdges_F; fEdge++) + fEdgeToEdgeID[fEdge] = indexToEdgeID_F[fEdge]; +#endif nEdges = edgeToFEdge.size(); indexToEdgeID.resize(nEdges); floatingEdgesIds.clear(); floatingEdgesIds.reserve(nEdges); + int maxEdgeID=std::numeric_limits::min(), minEdgeID=std::numeric_limits::max(), maxGlobalEdgeID, minGlobalEdgeID; for (int index = 0; index < nEdges; index++) { int fEdge = edgeToFEdge[index]; indexToEdgeID[index] = fEdgeToEdgeID[fEdge]; + maxEdgeID = (indexToEdgeID[index] > maxEdgeID) ? indexToEdgeID[index] : maxEdgeID; + minEdgeID = (indexToEdgeID[index] < minEdgeID) ? indexToEdgeID[index] : minEdgeID; if((floatingEdgesMask_F[fEdge]!=0)&&(index fCellsToSend; fCellsToSend.reserve(nCellsSolve_F); @@ -542,21 +629,42 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel std::vector fCellToVertexID(nCells_F, NotAnId); fCellsToReceive.clear(); + vertexProcIDs.clear(); + + std::vector verticesProcIds(nCells_F, NotAnId); fCellsToReceive.reserve(nCells_F - nCellsSolve_F); + vertexProcIDs.reserve(nCells_F); for (int i = 0; i < nCells_F; i++) { bool isMine = i < nCellsSolve_F; bool belongsToLocalTriangle = false; bool belongsToAnyTriangle = false; int nEdg = nEdgesOnCells_F[i]; + bool invalidVertex = false; + int minTriangleProcId = std::numeric_limits::max(); + bool nodeOwnedByATriaProc = false; for (int j = 0; j < nEdg; j++) { ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F * i + j] - 1); - ID iTria = fVertexToTriangle[fVertex]; + if(fVertex >= nVertices_F) { + invalidVertex = true; + break; + } + ID triaId = fVertexToTriangleID[fVertex]; - belongsToLocalTriangle = belongsToLocalTriangle || (iTria != NotAnId); + belongsToLocalTriangle = belongsToLocalTriangle || (trianglesProcIds[fVertex] == me); belongsToAnyTriangle = belongsToAnyTriangle || (triaId != NotAnId); + + if(triaId != NotAnId) { + nodeOwnedByATriaProc = nodeOwnedByATriaProc || (trianglesProcIds[fVertex] == fCellsProcIds[i]); + minTriangleProcId = (trianglesProcIds[fVertex] < minTriangleProcId) ? trianglesProcIds[fVertex] : minTriangleProcId; + } } + if(invalidVertex) continue; + + if(belongsToAnyTriangle) + verticesProcIds[i] = nodeOwnedByATriaProc ? fCellsProcIds[i] : minTriangleProcId; + if (belongsToAnyTriangle && isMine) { fCellsToSend.push_back(i); if (!belongsToLocalTriangle) @@ -566,42 +674,73 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel if (belongsToLocalTriangle) { fCellToVertex[i] = vertexToFCell.size(); vertexToFCell.push_back(i); + vertexProcIDs.push_back(reduced_ranks[verticesProcIds[i]]); } } + nVertices = vertexToFCell.size(); - //Compute the global number of vertices, and the localOffset on the local processor, such that a globalID = localOffset + index + //If useContiguousIDs is defined, create a contiguous list of global FE vertices IDs, otherwise use MPAS cells IDs +#ifdef useContiguousIDs + //Compute the localOffset on the local processor, such that a globalID = localOffset + index + int nGlobalVertices; computeLocalOffset(fCellsToSend.size(), localOffset, nGlobalVertices); - - //Communicate the globalIDs, computed locally, to the other processors. for (int index = 0; index < int(fCellsToSend.size()); index++) fCellToVertexID[fCellsToSend[index]] = index + localOffset; - + //Communicate the globalIDs, computed locally, to the other processors. allToAll(fCellToVertexID, sendCellsList_F, recvCellsList_F); +#else + for (int fcell = 0; fcell < nCells_F; fcell++) + fCellToVertexID[fcell] = indexToCellID_F[fcell]; +#endif - nVertices = vertexToFCell.size(); - int vertexColumnShift = (Ordering == 1) ? 1 : nGlobalVertices; - int vertexLayerShift = (Ordering == 0) ? 1 : nLayers + 1; - std::cout << "\n nvertices: " << nVertices << " " << nGlobalVertices << "\n" - << std::endl; - + int maxVertexID=std::numeric_limits::min(), minVertexID=std::numeric_limits::max(), maxGlobalVertexID, minGlobalVertexID; indexToVertexID.resize(nVertices); + for (int index = 0; index < nVertices; index++) { + int fCell = vertexToFCell[index]; + indexToVertexID[index] = fCellToVertexID[fCell]; + maxVertexID = (indexToVertexID[index] > maxVertexID) ? indexToVertexID[index] : maxVertexID; + minVertexID = (indexToVertexID[index] < minVertexID) ? indexToVertexID[index] : minVertexID; + } + + MPI_Allreduce(&maxVertexID, &maxGlobalVertexID, 1, MPI_INT, MPI_MAX, comm); + MPI_Allreduce(&minVertexID, &minGlobalVertexID, 1, MPI_INT, MPI_MIN, comm); + globalVertexStride = maxGlobalVertexID - minGlobalVertexID + 1; + + + int vertexColumnShift = (Ordering == 1) ? 1 : globalVertexStride; + int vertexLayerShift = (Ordering == 0) ? 1 : nLayers + 1; dirichletNodesIDs.clear(); dirichletNodesIDs.reserve(nVertices); //need to improve storage efficiency + isVertexBoundary.assign(nVertices, false); for (int index = 0; index < nVertices; index++) { int fCell = vertexToFCell[index]; - indexToVertexID[index] = fCellToVertexID[fCell]; for(int il=0; il< nLayers+1; ++il) { int imask_F = il+(nLayers+1)*fCell; if(dirichletCellsMask_F[imask_F]!=0) dirichletNodesIDs.push_back((nLayers-il)*vertexColumnShift+indexToVertexID[index]*vertexLayerShift); } + + int nEdg = nEdgesOnCells_F[fCell]; + int j = 0; + bool isBoundary; + do { + int fVertex = verticesOnCell_F[maxNEdgesOnCell_F * fCell + j++] - 1; + isBoundary = !(verticesMask_F[fVertex] & dynamic_ice_bit_value); + } while ((j < nEdg) && (!isBoundary)); + isVertexBoundary[index] = isBoundary; } + // because we change the ownership of some vertices, we need to first communicate back to the processors that used to own those vertices + // the data of the newly owned vertices. We do this by defining "reversed" send and receive lists, communicate back using that list, and + // then communicate forward with the usual send and receive lists. + // We could join these two step in one communication, but for the moment we do that separately. + // We need to communicate info about the vertices when we get the ice velocity on vertices form the velocity solver/ + createReverseCellsExchangeLists(sendCellsListReversed, recvCellsListReversed, - fVertexToTriangleID, fCellToVertexID); + verticesProcIds, indexToCellID_F); //construct the local vector vertices on triangles making sure the area is positive verticesOnTria.resize(nTriangles * 3); @@ -620,96 +759,24 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel std::swap(verticesOnTria[3 * index + 1], verticesOnTria[3 * index + 2]); } - //construct the local vector vertices on edges - trianglesPositionsOnEdge.resize(2 * nEdges); - isVertexBoundary.assign(nVertices, false); - verticesOnEdge.resize(2 * nEdges); - - //contains the local id of a triangle and the global id of the edges of the triangle. - //dataForGhostTria[4*i] contains the triangle id - //dataForGhostTria[4*i+1+k] contains the global id of the edge (at position k = 0,1,2) of the triangle. - //Possible Optimization: for our purposes it would be enough to store two of the three edges of a triangle. - std::vector dataForGhostTria(nVertices_F * 4, NotAnId); - - //* - for (int iV = 0; iV < nVertices; iV++) { - int fCell = vertexToFCell[iV]; - int nEdg = nEdgesOnCells_F[fCell]; - int j = 0; - bool isBoundary; - do { - int fVertex = verticesOnCell_F[maxNEdgesOnCell_F * fCell + j++] - 1; - isBoundary = !(verticesMask_F[fVertex] & dynamic_ice_bit_value); - } while ((j < nEdg) && (!isBoundary)); - isVertexBoundary[iV] = isBoundary; - } - /*/ - for(int index=0; index= 0) { @@ -856,14 +918,9 @@ void get_prism_velocity_on_FEdges(double * uNormal, UInt nPoints3D = nCells_F * (nLayers + 1); - // Loop over all edges of the triangulation - MPAS will decide which edges it should use. - for (int i = 0; i < nEdges; i++) { - - //identifying vertices on the edge - ID lId0 = verticesOnEdge[2 * i]; - ID lId1 = verticesOnEdge[2 * i + 1]; - int iCell0 = vertexToFCell[lId0]; - int iCell1 = vertexToFCell[lId1]; + for (int iEdge = 0; iEdge < nEdgesSolve_F; iEdge++) { + int iCell0 = cellsOnEdge_F[2 * iEdge] - 1; + int iCell1 = cellsOnEdge_F[2 * iEdge + 1] - 1; //computing normal to the cell edge (dual of triangular edge) double nx = xCell_F[iCell1] - xCell_F[iCell0]; @@ -873,9 +930,14 @@ void get_prism_velocity_on_FEdges(double * uNormal, ny /= n; //identifying triangles that shares the edge - ID iEdge = edgeToFEdge[i]; ID fVertex0 = verticesOnEdge_F[2 * iEdge] - 1; ID fVertex1 = verticesOnEdge_F[2 * iEdge + 1] - 1; + + + int iTria0 = fVertexToTriangleID[fVertex0]; + int iTria1 = fVertexToTriangleID[fVertex1]; + if((iTria0 == NotAnId) && (iTria1 == NotAnId)) continue; + double t0[2*3], t1[2*3]; //t0[0] contains the x-coords of vertices of triangle 0 and t0[1] its y-coords. for (int j = 0; j < 3; j++) { int iCell = cellsOnVertex_F[3 * fVertex0 + j] - 1; @@ -906,7 +968,7 @@ void get_prism_velocity_on_FEdges(double * uNormal, for (int j = 0; j < 3; j++) iCells[j] = cellsOnVertex_F[3 * fVertex1 + j] - 1; } - else if(i& velocityOnVertices, } } -void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, +void createReverseVerticesExchangeLists(exchangeList_Type& sendListReverse_F, exchangeList_Type& receiveListReverse_F, - const std::vector& fVertexToTriangleID, - const std::vector& fCellToVertexID) { + const std::vector& trianglesProcIds, const int* indexToVertexID_F) { sendListReverse_F.clear(); receiveListReverse_F.clear(); - //std::map > sendMap; std::map > sendMap, receiveMap; - std::vector cellsProcId(nCells_F), trianglesProcIds(nVertices_F); - getProcIds(cellsProcId, recvCellsList_F); - getProcIds(trianglesProcIds, recvVerticesList_F); - - for (int i = 0; i < nVertices; i++) { - int iCell = vertexToFCell[i]; - if (iCell < nCellsSolve_F) - continue; - bool belongToTriaOnSameProc = false; - int j(0); - int nEdg = nEdgesOnCells_F[iCell]; - do { - ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F * iCell + j] - 1); - ID triaId = fVertexToTriangleID[fVertex]; - belongToTriaOnSameProc = (triaId != NotAnId) - && (trianglesProcIds[fVertex] == cellsProcId[iCell]); - } while ((belongToTriaOnSameProc == false) && (++j < nEdg)); - if (!belongToTriaOnSameProc) { - sendMap[cellsProcId[iCell]].insert( - std::make_pair(fCellToVertexID[iCell], iCell)); + std::vector verticesProcIds(nVertices_F); + getProcIds(verticesProcIds, recvVerticesList_F); + int me; + MPI_Comm_rank(comm, &me); + for (int i = 0; i < nTriangles; i++) { + int iVertex = triangleToFVertex[i]; + if ((iVertex >= nVerticesSolve_F) && (trianglesProcIds[iVertex] == me)) { + sendMap[verticesProcIds[iVertex]].insert( + std::make_pair(indexToVertexID_F[iVertex]-1, iVertex)); } } @@ -1038,16 +1086,10 @@ void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, exchange(it->first, &sendVec[0], &sendVec[0] + sendVec.size())); } - for (UInt i = 0; i < fCellsToReceive.size(); i++) { - int iCell = fCellsToReceive[i]; - int nEdg = nEdgesOnCells_F[iCell]; - for (int j = 0; j < nEdg; j++) { - ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F * iCell + j] - 1); - ID triaId = fVertexToTriangleID[fVertex]; - if (triaId != NotAnId) { - receiveMap[trianglesProcIds[fVertex]].insert( - std::make_pair(fCellToVertexID[iCell], iCell)); - } + for (int iVertex = 0; iVertex < nVerticesSolve_F; iVertex++) { + if((trianglesProcIds[iVertex] != NotAnId) && (trianglesProcIds[iVertex] != me)) { + receiveMap[trianglesProcIds[iVertex]].insert( + std::make_pair(indexToVertexID_F[iVertex]-1, iVertex)); } } @@ -1064,33 +1106,21 @@ void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, } } -void createReverseEdgesExchangeLists(exchangeList_Type& sendListReverse_F, +void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, exchangeList_Type& receiveListReverse_F, - const std::vector& fVertexToTriangleID, - const std::vector& fEdgeToEdgeID) { + const std::vector& verticesProcIds, const int* indexToCellID_F) { sendListReverse_F.clear(); receiveListReverse_F.clear(); - //std::map > sendMap; std::map > sendMap, receiveMap; - std::vector edgesProcId(nEdges_F), trianglesProcIds(nVertices_F); - getProcIds(edgesProcId, recvEdgesList_F); - getProcIds(trianglesProcIds, recvVerticesList_F); - - for (int i = 0; i < nEdges; i++) { - int iEdge = edgeToFEdge[i]; - if (iEdge < nEdgesSolve_F) - continue; - bool belongToTriaOnSameProc = false; - int j(0); - do { - ID fVertex(verticesOnEdge_F[2 * iEdge + j] - 1); - ID triaId = fVertexToTriangleID[fVertex]; - belongToTriaOnSameProc = (triaId != NotAnId) - && (trianglesProcIds[fVertex] == edgesProcId[iEdge]); - } while ((belongToTriaOnSameProc == false) && (++j < 2)); - if (!belongToTriaOnSameProc) { - sendMap[edgesProcId[iEdge]].insert( - std::make_pair(fEdgeToEdgeID[iEdge], iEdge)); + std::vector CellsProcIds(nCells_F); + getProcIds(CellsProcIds, recvCellsList_F); + int me; + MPI_Comm_rank(comm, &me); + for (int i = 0; i < nVertices; i++) { + int iCell = vertexToFCell[i]; + if ((iCell >= nCellsSolve_F) && (verticesProcIds[iCell] == me)) { + sendMap[CellsProcIds[iCell]].insert( + std::make_pair(indexToCellID_F[iCell]-1, iCell)); } } @@ -1105,15 +1135,10 @@ void createReverseEdgesExchangeLists(exchangeList_Type& sendListReverse_F, exchange(it->first, &sendVec[0], &sendVec[0] + sendVec.size())); } - for (UInt i = 0; i < edgesToReceive.size(); i++) { - int iEdge = edgesToReceive[i]; - for (int j = 0; j < 2; j++) { - ID fVertex(verticesOnEdge_F[2 * iEdge + j] - 1); - ID triaId = fVertexToTriangleID[fVertex]; - if (triaId != NotAnId) { - receiveMap[trianglesProcIds[fVertex]].insert( - std::make_pair(fEdgeToEdgeID[iEdge], iEdge)); - } + for (int iCell = 0; iCell < nCellsSolve_F; iCell++) { + if((verticesProcIds[iCell] != NotAnId) && (verticesProcIds[iCell] != me)) { + receiveMap[verticesProcIds[iCell]].insert( + std::make_pair(indexToCellID_F[iCell]-1, iCell)); } } @@ -1154,19 +1179,6 @@ void mapCellsToVertices(const std::vector& velocityOnCells, } } -bool isGhostTriangle(int i, double relTol) { - double x[3], y[3], area; - - for (int j = 0; j < 3; j++) { - int iCell = cellsOnVertex_F[3 * i + j] - 1; - x[j] = xCell_F[iCell]; - y[j] = yCell_F[iCell]; - } - - area = std::fabs(signedTriangleArea(x, y)); - return false; //(std::fabs(areaTriangle_F[i]-area)/areaTriangle_F[i] > relTol); -} - double signedTriangleArea(const double* x, const double* y) { double u[2] = { x[1] - x[0], y[1] - y[0] }; double v[2] = { x[2] - x[0], y[2] - y[0] }; @@ -1273,11 +1285,6 @@ void import2DFields(std::map bdExtensionMap, double const* bedTopograp bool foundNeighbor = false; for (int j = 0; j < nEdg; j++) { int fEdge = edgesOnCell_F[maxNEdgesOnCell_F * fCell + j] - 1; - //bool keep = (mask[verticesOnEdge_F[2 * fEdge] - 1] & dynamic_ice_bit_value) - // && (mask[verticesOnEdge_F[2 * fEdge + 1] - 1] & dynamic_ice_bit_value); - //if (!keep) - // continue; - int c0 = cellsOnEdge_F[2 * fEdge] - 1; int c1 = cellsOnEdge_F[2 * fEdge + 1] - 1; c = (fCellToVertex[c0] == iV) ? c1 : c0; @@ -1426,7 +1433,9 @@ void exportDissipationHeat(double * dissipationHeat_F) { } } - //allToAll (dissipationHeat_F, &sendVerticesListReversed, &recvEdgesListReversed, nLayers+1); +#ifdef changeTrianglesOwnership + allToAll (dissipationHeat_F, &sendVerticesListReversed, &recvVerticesListReversed, nLayers); +#endif allToAll (dissipationHeat_F, sendVerticesList_F, recvVerticesList_F, nLayers); } @@ -1441,7 +1450,7 @@ void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id) { int nonEmpty = int(nLocalEntities > 0); MPI_Allgather(&nonEmpty, 1, MPI_INT, &haveElements[0], 1, MPI_INT, comm); std::vector ranks; - reduced_ranks.resize(numProcs,0); + reduced_ranks.assign(numProcs,-1); for (int i = 0; i < numProcs; i++) { if (haveElements[i]) { reduced_ranks[i] = ranks.size(); @@ -1701,7 +1710,6 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep int nEdg = nEdgesOnCells_F[fCell]; int me; MPI_Comm_rank(comm, &me); - procIds.reserve(nEdg); for(int i=0; i& betaData, const std::vector& indexToVertexID); -extern void velocity_solver_extrude_3d_grid__(int nLayers, int nGlobalTriangles, - int nGlobalVertices, int nGlobalEdges, int Ordering, MPI_Comm reducedComm, +extern void velocity_solver_extrude_3d_grid__( + int nLayers, int globalTriangleStride, int globalVertexStride, int globalEdgeStride, + int Ordering, MPI_Comm reducedComm, const std::vector& indexToVertexID, const std::vector& mpasIndexToVertexID, + const std::vector& vertexProcIDs, const std::vector& verticesCoords, - const std::vector& isVertexBoundary, const std::vector& verticesOnTria, const std::vector> procsSharingVertices, const std::vector& isBoundaryEdge, const std::vector& trianglesOnEdge, - const std::vector& trianglesPositionsOnEdge, const std::vector& verticesOnEdge, const std::vector& indexToEdgeID, const std::vector& indexToTriangleID, @@ -195,8 +196,6 @@ extern void velocity_solver_export_fo_velocity__(MPI_Comm reducedComm); exchangeList_Type unpackMpiArray(int const* array); -bool isGhostTriangle(int i, double relTol = 1e-1); - double signedTriangleArea(const double* x, const double* y); double signedTriangleArea(const double* x, const double* y, const double* z); @@ -232,15 +231,14 @@ void get_prism_velocity_on_FEdges(double* uNormal, int initialize_iceProblem(int nTriangles); -void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, +void createReverseVerticesExchangeLists(exchangeList_Type& sendListReverse_F, exchangeList_Type& receiveListReverse_F, - const std::vector& fVertexToTriangleID, - const std::vector& fCellToVertexID); + const std::vector& fVertexToTriangleID, const int* indexToVertexID_F); -void createReverseEdgesExchangeLists(exchangeList_Type& sendListReverse_F, +void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, exchangeList_Type& receiveListReverse_F, - const std::vector& fVertexToTriangleID, - const std::vector& fEdgeToEdgeID); + const std::vector& fCellToVertexID, + const int* indexToCellID_F); void mapCellsToVertices(const std::vector& velocityOnCells, std::vector& velocityOnVertices, int fieldDim, int numLayers, diff --git a/src/core_landice/mode_forward/mpas_li_velocity_external.F b/src/core_landice/mode_forward/mpas_li_velocity_external.F index 8587bd5d65..f1fa8bd529 100644 --- a/src/core_landice/mode_forward/mpas_li_velocity_external.F +++ b/src/core_landice/mode_forward/mpas_li_velocity_external.F @@ -327,7 +327,7 @@ subroutine li_velocity_external_block_init(block, err) call velocity_solver_set_grid_data(nCells, nEdges, nVertices, nVertInterfaces, & nCellsSolve, nEdgesSolve, nVerticesSolve, maxNEdgesOnCell, radius, & cellsOnEdge, cellsOnVertex, verticesOnCell, verticesOnEdge, edgesOnCell, & - nEdgesOnCell, indexToCellID, & + nEdgesOnCell, indexToCellID, indexToEdgeID, indexToVertexID, & xCell, yCell, zCell, xVertex, yVertex, zVertex, areaTriangle, & sendCellsArray, recvCellsArray, & sendEdgesArray, recvEdgesArray, & @@ -1022,7 +1022,7 @@ subroutine generate_fem_grid(config_velocity_solver, vertexMask, cellMask, diric case ('FO') ! =============================================== #ifdef USE_EXTERNAL_FIRSTORDER call mpas_timer_start("velocity_solver_extrude_3d_grid") - call velocity_solver_extrude_3d_grid(layerThicknessFractions, lowerSurface, thickness) + call velocity_solver_extrude_3d_grid(layerThicknessFractions) call mpas_timer_stop("velocity_solver_extrude_3d_grid") call mpas_timer_start("velocity_solver_init_FO") call velocity_solver_init_FO(layerThicknessFractions) From 66e0b31b1223fd586ec2f5c17c951b9d7f1d7e52 Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Mon, 11 Nov 2019 13:53:49 -0700 Subject: [PATCH 7/9] Remove option to send to Albany list of contiguous Ids --- .../Interface_velocity_solver.cpp | 110 +++--------------- .../Interface_velocity_solver.hpp | 4 - 2 files changed, 13 insertions(+), 101 deletions(-) diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.cpp b/src/core_landice/mode_forward/Interface_velocity_solver.cpp index f586626585..cbbdfcfcac 100644 --- a/src/core_landice/mode_forward/Interface_velocity_solver.cpp +++ b/src/core_landice/mode_forward/Interface_velocity_solver.cpp @@ -22,7 +22,6 @@ distributed with this code, or at http://mpas-dev.github.com/license.html //! Namespaces // =================================================== -//#define useContiguousIDs #define changeTrianglesOwnership @@ -62,7 +61,7 @@ int original_stderr; // the location of stderr before we captured it int Interface_stdout; // the location of stdout as we use it here //void *phgGrid = 0; -std::vector edgesToReceive, fCellsToReceive, indexToTriangleID, +std::vector indexToTriangleID, verticesOnTria, trianglesOnEdge, verticesOnEdge, trianglesProcIds, reduced_ranks; std::vector indexToVertexID, vertexToFCell, vertexProcIDs, triangleToFVertex, indexToEdgeID, edgeToFEdge, @@ -288,13 +287,8 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower albany_error, dt); *error=albany_error; } - - exportDissipationHeat(dissipation_heat_F); - std::vector mpasIndexToVertexID(nVertices); - for (int i = 0; i < nVertices; i++) { - mpasIndexToVertexID[i] = indexToCellID_F[vertexToFCell[i]]; - } + exportDissipationHeat(dissipation_heat_F); mapVerticesToCells(velocityOnVertices, &velocityOnCells[0], 2, nLayers, Ordering); @@ -363,7 +357,6 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel numProcs + 1); - // First, we compute the FE triangles belonging to this processor. // If changeTrianglesOwnership is not define, the triangles belonging to this // processor will be the subset of the triangles (MPAS vertices) owned by this proc @@ -447,22 +440,12 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel //Initialize the ice sheet problem with the number of FE triangles on this prov initialize_iceProblem(nTriangles); - //If useContiguousIDs is defined, create a contiguous list of global triangle IDs, otherwise use MPAS vertices IDs + //Create a list of global IDs for FE triangles, using MPAS vertices IDs fVertexToTriangleID.assign(nVertices_F, NotAnId); -#ifdef useContiguousIDs - // Compute the localOffset on the local processor, such that a globalID = localOffset + index - int localOffset(0); - int nGlobalTriangles = 0; - computeLocalOffset(nTriangles, localOffset, nGlobalTriangles); - for (int index(0); index < nTriangles; index++) - fVertexToTriangleID[triangleToFVertex[index]] = index + localOffset; -#else for (int index(0); index < nTriangles; index++) { fVertexToTriangleID[triangleToFVertex[index]] = indexToVertexID_F[triangleToFVertex[index]]; } -#endif - //Communicate the local and global triangle IDs to the other processors. #ifdef changeTrianglesOwnership // because we change the ownership of some triangles, we need to first communicate back to the processors that used to own those triangles // the data of the newly owned triangles. We do this by defining "reversed" send and receive lists, communicate back using those lists, and @@ -476,7 +459,6 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel allToAll(fVertexToTriangleID, sendVerticesList_F, recvVerticesList_F); allToAll(fVertexToTriangle, sendVerticesList_F, recvVerticesList_F); - //we define the vector of global triangles Ids and compute the stride between the largest and the smallest Id globally //This will be needed by the velocity solver to create the 3D FE mesh. indexToTriangleID.resize(nTriangles); @@ -493,18 +475,15 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel // Second, we compute the FE edges belonging to the FE triangles owned by this processor. // We first compute boundary edges, and then all the other edges. - std::vector fEdgeToEdge(nEdges_F), edgesToSend; + std::vector fEdgeToEdge(nEdges_F); std::vector fEdgeToEdgeID(nEdges_F, NotAnId); - edgesToReceive.clear(); edgeToFEdge.clear(); isBoundaryEdge.clear(); trianglesOnEdge.clear(); - edgesToReceive.reserve(nEdges_F - nEdgesSolve_F); edgeToFEdge.reserve(nEdges_F); trianglesOnEdge.reserve(nEdges_F * 2); - edgesToSend.reserve(nEdgesSolve_F); isBoundaryEdge.reserve(nEdges_F); //we compute boundary edges (boundary edges must be the first edges) @@ -574,28 +553,10 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel trianglesOnEdge.push_back(iTria2); isBoundaryEdge.push_back(false); } - - if (belongsToAnyTriangle && isMine) { - edgesToSend.push_back(i); - if (!belongsToLocalTriangle) - edgesToReceive.push_back(i); - } - } - //If useContiguousIDs is defined, create a contiguous list of global edge IDs, otherwise use MPAS edges IDs -#ifdef useContiguousIDs - //Compute the localOffset on the local processor, such that a globalID = localOffset + index - int nGlobalEdges; - computeLocalOffset(edgesToSend.size(), localOffset, nGlobalEdges); - for (ID index = 0; index < edgesToSend.size(); index++) - fEdgeToEdgeID[edgesToSend[index]] = index + localOffset; - //Communicate the globalIDs, computed locally, to the other processors. - allToAll(fEdgeToEdgeID, sendEdgesList_F, recvEdgesList_F); -#else for (int fEdge = 0; fEdge < nEdges_F; fEdge++) fEdgeToEdgeID[fEdge] = indexToEdgeID_F[fEdge]; -#endif nEdges = edgeToFEdge.size(); indexToEdgeID.resize(nEdges); @@ -619,8 +580,6 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel // Third, we compute the FE vertices belonging to the FE triangles owned by this processor. // We need to make sure that an FE vertex is owned by a proc that owns a FE triangle that contain that vertex // Otherwise we might end up with weird situation where a vertex could belong to a process with no associated triangle. - std::vector fCellsToSend; - fCellsToSend.reserve(nCellsSolve_F); vertexToFCell.clear(); vertexToFCell.reserve(nCells_F); @@ -628,12 +587,10 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel fCellToVertex.assign(nCells_F, NotAnId); std::vector fCellToVertexID(nCells_F, NotAnId); - fCellsToReceive.clear(); vertexProcIDs.clear(); std::vector verticesProcIds(nCells_F, NotAnId); - fCellsToReceive.reserve(nCells_F - nCellsSolve_F); vertexProcIDs.reserve(nCells_F); for (int i = 0; i < nCells_F; i++) { bool isMine = i < nCellsSolve_F; @@ -665,12 +622,6 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel if(belongsToAnyTriangle) verticesProcIds[i] = nodeOwnedByATriaProc ? fCellsProcIds[i] : minTriangleProcId; - if (belongsToAnyTriangle && isMine) { - fCellsToSend.push_back(i); - if (!belongsToLocalTriangle) - fCellsToReceive.push_back(i); - } - if (belongsToLocalTriangle) { fCellToVertex[i] = vertexToFCell.size(); vertexToFCell.push_back(i); @@ -679,21 +630,8 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel } nVertices = vertexToFCell.size(); - //If useContiguousIDs is defined, create a contiguous list of global FE vertices IDs, otherwise use MPAS cells IDs -#ifdef useContiguousIDs - //Compute the localOffset on the local processor, such that a globalID = localOffset + index - int nGlobalVertices; - computeLocalOffset(fCellsToSend.size(), localOffset, nGlobalVertices); - for (int index = 0; index < int(fCellsToSend.size()); index++) - fCellToVertexID[fCellsToSend[index]] = index + localOffset; - //Communicate the globalIDs, computed locally, to the other processors. - allToAll(fCellToVertexID, sendCellsList_F, recvCellsList_F); -#else for (int fcell = 0; fcell < nCells_F; fcell++) - fCellToVertexID[fcell] = indexToCellID_F[fcell]; -#endif - - + fCellToVertexID[fcell] = indexToCellID_F[fcell]; int maxVertexID=std::numeric_limits::min(), minVertexID=std::numeric_limits::max(), maxGlobalVertexID, minGlobalVertexID; indexToVertexID.resize(nVertices); @@ -793,10 +731,6 @@ void velocity_solver_extrude_3d_grid(double const* levelsRatio_F) { levelsNormalizedThickness[i + 1] = levelsNormalizedThickness[i] + layersRatio[i]; - std::vector mpasIndexToVertexID(nVertices); - for (int i = 0; i < nVertices; i++) - mpasIndexToVertexID[i] = indexToCellID_F[vertexToFCell[i]]; - //construct the local vector of coordinates std::vector verticesCoords(3 * nVertices); @@ -812,7 +746,7 @@ void velocity_solver_extrude_3d_grid(double const* levelsRatio_F) { procsSharingVertex(i, procsSharingVertices[i]); velocity_solver_extrude_3d_grid__(nLayers, globalTriangleStride, globalVertexStride, - globalEdgeStride, Ordering, reducedComm, indexToVertexID, mpasIndexToVertexID, vertexProcIDs, + globalEdgeStride, Ordering, reducedComm, indexToVertexID, vertexProcIDs, verticesCoords, verticesOnTria, procsSharingVertices, isBoundaryEdge, trianglesOnEdge, verticesOnEdge, indexToEdgeID, indexToTriangleID, dirichletNodesIDs, floatingEdgesIds); @@ -992,7 +926,7 @@ void get_prism_velocity_on_FEdges(double * uNormal, //Scale the coordinates so that they sum to 1. for(int j=0; j<3; j++) bcoords[j] /= sum; - } + } else { //error, edge midpont does not belong to either triangle std::cout << "Error, edge midpont does not belong to either triangle" << std::endl; @@ -1249,7 +1183,7 @@ void import2DFields(std::map bdExtensionMap, double const* bedTopograp int v = verticesOnTria[iVertex + 3 * index]; int iCell = vertexToFCell[v]; //compute temperature by averaging temperature values of triangles vertices where ice is present - if (cellsMask_F[iCell] & ice_present_bit_value) { + if (cellsMask_F[iCell] & ice_present_bit_value) { temperature += temperature_F[iCell * nLayers + ilReversed]; nPoints++; } @@ -1463,24 +1397,6 @@ void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id) { MPI_Comm_create(comm, reduced_group_id, &reduced_comm_id); } -void computeLocalOffset(int nLocalEntities, int& localOffset, - int& nGlobalEntities) { - int numProcs, me; - MPI_Comm_size(comm, &numProcs); - MPI_Comm_rank(comm, &me); - std::vector offsetVec(numProcs); - - MPI_Allgather(&nLocalEntities, 1, MPI_INT, &offsetVec[0], 1, MPI_INT, comm); - - localOffset = 0; - for (int i = 0; i < me; i++) - localOffset += offsetVec[i]; - - nGlobalEntities = localOffset; - for (int i = me; i < numProcs; i++) - nGlobalEntities += offsetVec[i]; -} - void getProcIds(std::vector& field, int const * recvArray) { int me; MPI_Comm_rank(comm, &me); @@ -1799,9 +1715,9 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep // individual field values // Call needed functions to process MPAS fields to Albany units/format - + std::map bdExtensionMap; // local map to be created by import2DFields - import2DFields(bdExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, + import2DFields(bdExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, temperature_F, smb_F, minThickness); import2DFieldsObservations(bdExtensionMap, @@ -1837,13 +1753,13 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep outfile.open ("temperature.ascii", std::ios::out | std::ios::trunc); if (outfile.is_open()) { outfile << nTriangles << " " << nLayers << "\n"; //number of triangles and number of layers on first line - + double midLayer = 0; for (int il = 0; il < nLayers; il++) { //sigma coordinates for temperature midLayer += layersRatio[il]/2.0; outfile << midLayer << "\n"; midLayer += layersRatio[il]/2.0; - } + } for(int il = 0; il& indexToVertexID, - const std::vector& mpasIndexToVertexID, const std::vector& vertexProcIDs, const std::vector& verticesCoords, const std::vector& verticesOnTria, @@ -247,9 +246,6 @@ void mapCellsToVertices(const std::vector& velocityOnCells, void mapVerticesToCells(const std::vector& velocityOnVertices, double* velocityOnCells, int fieldDim, int numLayers, int ordering); -void computeLocalOffset(int nLocalEntities, int& localOffset, - int& nGlobalEntities); - void getProcIds(std::vector& field, int const* recvArray); void getProcIds(std::vector& field, exchangeList_Type const* recvList); From 43996c34a9202301fc0a8f55dbb0b350e595588b Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Mon, 11 Nov 2019 13:56:59 -0700 Subject: [PATCH 8/9] Interface_velocity_solver: consolidate functions to create reversed comm lists for vertices and cells into one function --- .../Interface_velocity_solver.cpp | 160 +++++------------- .../Interface_velocity_solver.hpp | 17 +- 2 files changed, 52 insertions(+), 125 deletions(-) diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.cpp b/src/core_landice/mode_forward/Interface_velocity_solver.cpp index cbbdfcfcac..671ed32bcd 100644 --- a/src/core_landice/mode_forward/Interface_velocity_solver.cpp +++ b/src/core_landice/mode_forward/Interface_velocity_solver.cpp @@ -66,8 +66,8 @@ std::vector indexToTriangleID, trianglesProcIds, reduced_ranks; std::vector indexToVertexID, vertexToFCell, vertexProcIDs, triangleToFVertex, indexToEdgeID, edgeToFEdge, fVertexToTriangleID, fCellToVertex, floatingEdgesIds, dirichletNodesIDs; -std::vector temperatureOnTetra, dissipationHeatOnTetra, velocityOnVertices, velocityOnCells, - elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, temperatureData, smbData, thicknessOnCells; +std::vector dissipationHeatOnPrisms, velocityOnVertices, velocityOnCells, + elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, temperatureDataOnPrisms, smbData, thicknessOnCells; std::vector isVertexBoundary, isBoundaryEdge; // only needed for creating ASCII mesh @@ -264,15 +264,13 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower if (!isDomainEmpty) { std::map bdExtensionMap; - import2DFields(bdExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, temperature_F, smb_F, minThickness); + importFields(bdExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, temperature_F, smb_F, minThickness); std::vector regulThk(thicknessData); for (int index = 0; index < nVertices; index++) regulThk[index] = std::max(1e-4, thicknessData[index]); - importP0Temperature(); - - dissipationHeatOnTetra.resize(3 * nLayers * indexToTriangleID.size()); + dissipationHeatOnPrisms.resize(nLayers * indexToTriangleID.size()); std::cout << "\n\nTimeStep: "<< *deltat << "\n\n"<< std::endl; @@ -283,7 +281,7 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower regulThk, levelsNormalizedThickness, elevationData, thicknessData, betaData, bedTopographyData, smbData, stiffnessFactorData, effecPressData, - temperatureOnTetra, dissipationHeatOnTetra, velocityOnVertices, + temperatureDataOnPrisms, dissipationHeatOnPrisms, velocityOnVertices, albany_error, dt); *error=albany_error; } @@ -451,7 +449,7 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel // the data of the newly owned triangles. We do this by defining "reversed" send and receive lists, communicate back using those lists, and // then communicate "forward" using the usual send and receive lists. // We could join these two step in one communication, but for the moment we do that separately - createReverseVerticesExchangeLists(sendVerticesListReversed, recvVerticesListReversed, trianglesProcIds, indexToVertexID_F); + createReverseExchangeLists(sendVerticesListReversed, recvVerticesListReversed, trianglesProcIds, indexToVertexID_F, recvVerticesList_F); allToAll(fVertexToTriangleID, &sendVerticesListReversed, &recvVerticesListReversed); allToAll(fVertexToTriangle, &sendVerticesListReversed, &recvVerticesListReversed); allToAll(trianglesProcIds, sendVerticesList_F, recvVerticesList_F); @@ -576,7 +574,6 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel MPI_Allreduce(&minEdgeID, &minGlobalEdgeID, 1, MPI_INT, MPI_MIN, comm); globalEdgeStride = maxGlobalEdgeID - minGlobalEdgeID + 1; - // Third, we compute the FE vertices belonging to the FE triangles owned by this processor. // We need to make sure that an FE vertex is owned by a proc that owns a FE triangle that contain that vertex // Otherwise we might end up with weird situation where a vertex could belong to a process with no associated triangle. @@ -677,8 +674,7 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel // We could join these two step in one communication, but for the moment we do that separately. // We need to communicate info about the vertices when we get the ice velocity on vertices form the velocity solver/ - createReverseCellsExchangeLists(sendCellsListReversed, recvCellsListReversed, - verticesProcIds, indexToCellID_F); + createReverseExchangeLists(sendCellsListReversed, recvCellsListReversed, verticesProcIds, indexToCellID_F, recvCellsList_F); //construct the local vector vertices on triangles making sure the area is positive verticesOnTria.resize(nTriangles * 3); @@ -991,70 +987,22 @@ void mapVerticesToCells(const std::vector& velocityOnVertices, } } -void createReverseVerticesExchangeLists(exchangeList_Type& sendListReverse_F, - exchangeList_Type& receiveListReverse_F, - const std::vector& trianglesProcIds, const int* indexToVertexID_F) { - sendListReverse_F.clear(); - receiveListReverse_F.clear(); - std::map > sendMap, receiveMap; - std::vector verticesProcIds(nVertices_F); - getProcIds(verticesProcIds, recvVerticesList_F); - int me; - MPI_Comm_rank(comm, &me); - for (int i = 0; i < nTriangles; i++) { - int iVertex = triangleToFVertex[i]; - if ((iVertex >= nVerticesSolve_F) && (trianglesProcIds[iVertex] == me)) { - sendMap[verticesProcIds[iVertex]].insert( - std::make_pair(indexToVertexID_F[iVertex]-1, iVertex)); - } - } - - for (std::map >::const_iterator it = sendMap.begin(); - it != sendMap.end(); it++) { - std::vector sendVec(it->second.size()); - int i = 0; - for (std::map::const_iterator iter = it->second.begin(); - iter != it->second.end(); iter++) - sendVec[i++] = iter->second; - sendListReverse_F.push_back( - exchange(it->first, &sendVec[0], &sendVec[0] + sendVec.size())); - } - - for (int iVertex = 0; iVertex < nVerticesSolve_F; iVertex++) { - if((trianglesProcIds[iVertex] != NotAnId) && (trianglesProcIds[iVertex] != me)) { - receiveMap[trianglesProcIds[iVertex]].insert( - std::make_pair(indexToVertexID_F[iVertex]-1, iVertex)); - } - } - - for (std::map >::const_iterator it = - receiveMap.begin(); it != receiveMap.end(); it++) { - std::vector receiveVec(it->second.size()); - int i = 0; - for (std::map::const_iterator iter = it->second.begin(); - iter != it->second.end(); iter++) - receiveVec[i++] = iter->second; - receiveListReverse_F.push_back( - exchange(it->first, &receiveVec[0], - &receiveVec[0] + receiveVec.size())); - } -} -void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, +void createReverseExchangeLists(exchangeList_Type& sendListReverse_F, exchangeList_Type& receiveListReverse_F, - const std::vector& verticesProcIds, const int* indexToCellID_F) { + const std::vector& newProcIds, const int* indexToID_F, exchangeList_Type const * recvList_F) { sendListReverse_F.clear(); receiveListReverse_F.clear(); std::map > sendMap, receiveMap; - std::vector CellsProcIds(nCells_F); - getProcIds(CellsProcIds, recvCellsList_F); + int nFEntities = newProcIds.size(); + std::vector procIds(nFEntities); + getProcIds(procIds, recvList_F); int me; MPI_Comm_rank(comm, &me); - for (int i = 0; i < nVertices; i++) { - int iCell = vertexToFCell[i]; - if ((iCell >= nCellsSolve_F) && (verticesProcIds[iCell] == me)) { - sendMap[CellsProcIds[iCell]].insert( - std::make_pair(indexToCellID_F[iCell]-1, iCell)); + for (int fEntity = 0; fEntity < nFEntities; fEntity++) { + if ((procIds[fEntity] != me) && (newProcIds[fEntity] == me)) { + sendMap[procIds[fEntity]].insert( + std::make_pair(indexToID_F[fEntity], fEntity)); } } @@ -1069,10 +1017,10 @@ void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, exchange(it->first, &sendVec[0], &sendVec[0] + sendVec.size())); } - for (int iCell = 0; iCell < nCellsSolve_F; iCell++) { - if((verticesProcIds[iCell] != NotAnId) && (verticesProcIds[iCell] != me)) { - receiveMap[verticesProcIds[iCell]].insert( - std::make_pair(indexToCellID_F[iCell]-1, iCell)); + for (int fEntity = 0; fEntity < nFEntities; fEntity++) { + if((procIds[fEntity] == me) && (newProcIds[fEntity] != NotAnId) && (newProcIds[fEntity] != me)) { + receiveMap[newProcIds[fEntity]].insert( + std::make_pair(indexToID_F[fEntity], fEntity)); } } @@ -1137,7 +1085,7 @@ double signedTriangleAreaOnSphere(const double* x, const double* y, } -void import2DFields(std::map bdExtensionMap, double const* bedTopography_F, double const * lowerSurface_F, double const * thickness_F, +void importFields(std::map bdExtensionMap, double const* bedTopography_F, double const * lowerSurface_F, double const * thickness_F, double const * beta_F, double const* stiffnessFactor_F, double const* effecPress_F, double const * temperature_F, double const * smb_F, double eps) { @@ -1152,7 +1100,7 @@ void import2DFields(std::map bdExtensionMap, double const* bedTopograp if (effecPress_F != 0) effecPressData.assign(nVertices, 1e10); if(temperature_F != 0) - temperatureData.assign(nLayers * nTriangles, 1e10); + temperatureDataOnPrisms.assign(nLayers * nTriangles, 1e10); if (smb_F != 0) smbData.assign(nVertices, 1e10); @@ -1173,6 +1121,8 @@ void import2DFields(std::map bdExtensionMap, double const* bedTopograp effecPressData[index] = effecPress_F[iCell] / unit_length; } + int lElemColumnShift = (Ordering == 1) ? 1 : nTriangles; + int elemLayerShift = (Ordering == 0) ? 1 : nLayers; if(temperature_F != 0) { for (int index = 0; index < nTriangles; index++) { for (int il = 0; il < nLayers; il++) { @@ -1189,14 +1139,13 @@ void import2DFields(std::map bdExtensionMap, double const* bedTopograp } } if (nPoints == 0) //if triangle is in an ice-free area, set the temperature to T0 - temperatureData[index+il*nTriangles] = T0; + temperatureDataOnPrisms[index*elemLayerShift + il*lElemColumnShift] = T0; else - temperatureData[index+il*nTriangles] = temperature / nPoints; + temperatureDataOnPrisms[index*elemLayerShift + il*lElemColumnShift] = temperature / nPoints; } } } - //extend thickness elevation and basal friction data to the border for floating vertices std::set::const_iterator iter; @@ -1317,7 +1266,7 @@ void import2DFieldsObservations(std::map bdExtensionMap, indexToCellIDData[index] = indexToCellID_F[iCell]; } - //extend to the border for floating vertices (using map created by import2DFields above) + //extend to the border for floating vertices (using map created by importFields above) for (std::map::iterator it = bdExtensionMap.begin(); it != bdExtensionMap.end(); ++it) { int iv = it->first; @@ -1338,33 +1287,15 @@ void import2DFieldsObservations(std::map bdExtensionMap, } } - -void importP0Temperature() { - int lElemColumnShift = (Ordering == 1) ? 3 : 3 * indexToTriangleID.size(); - int elemLayerShift = (Ordering == 0) ? 3 : 3 * nLayers; - temperatureOnTetra.resize(3 * nLayers * indexToTriangleID.size()); - for (int index = 0; index < nTriangles; index++) { - for (int il = 0; il < nLayers; il++) { - for (int k = 0; k < 3; k++) - temperatureOnTetra[index * elemLayerShift + il * lElemColumnShift + k] = - temperatureData[index+il*nTriangles]; - } - } -} - void exportDissipationHeat(double * dissipationHeat_F) { std::fill(dissipationHeat_F, dissipationHeat_F + nVertices_F * (nLayers), 0.); - int lElemColumnShift = (Ordering == 1) ? 3 : 3 * indexToTriangleID.size(); - int elemLayerShift = (Ordering == 0) ? 3 : 3 * nLayers; + int lElemColumnShift = (Ordering == 1) ? 1 : indexToTriangleID.size(); + int elemLayerShift = (Ordering == 0) ? 1 : nLayers; for (int index = 0; index < nTriangles; index++) { for (int il = 0; il < nLayers; il++) { int ilReversed = nLayers - il - 1; int fVertex = triangleToFVertex[index]; - double dissipationHeat = 0; - for (int k = 0; k < 3; k++) - dissipationHeat += dissipationHeatOnTetra[index * elemLayerShift + il * lElemColumnShift + k]/3; //TODO: should be weighted average based on tets volumes - dissipationHeat_F[fVertex * nLayers + ilReversed] = dissipationHeat; - + dissipationHeat_F[fVertex * nLayers + ilReversed] = dissipationHeatOnPrisms[index * elemLayerShift + il * lElemColumnShift]; } } #ifdef changeTrianglesOwnership @@ -1716,8 +1647,8 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep // individual field values // Call needed functions to process MPAS fields to Albany units/format - std::map bdExtensionMap; // local map to be created by import2DFields - import2DFields(bdExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, + std::map bdExtensionMap; // local map to be created by importFields + importFields(bdExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, temperature_F, smb_F, minThickness); import2DFieldsObservations(bdExtensionMap, @@ -1752,20 +1683,23 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep std::cout << "Writing temperature.ascii." << std::endl; outfile.open ("temperature.ascii", std::ios::out | std::ios::trunc); if (outfile.is_open()) { - outfile << nTriangles << " " << nLayers << "\n"; //number of triangles and number of layers on first line + outfile << nTriangles << " " << nLayers << "\n"; //number of triangles and number of layers on first line - double midLayer = 0; - for (int il = 0; il < nLayers; il++) { //sigma coordinates for temperature - midLayer += layersRatio[il]/2.0; - outfile << midLayer << "\n"; - midLayer += layersRatio[il]/2.0; - } + double midLayer = 0; + for (int il = 0; il < nLayers; il++) { //sigma coordinates for temperature + midLayer += layersRatio[il]/2.0; + outfile << midLayer << "\n"; + midLayer += layersRatio[il]/2.0; + } - for(int il = 0; il& smbData, const std::vector& stiffnessFactorData, const std::vector& effecPressData, - const std::vector& temperatureOnTetra, - std::vector& dissipationHeatOnTetra, + const std::vector& temperatureDataOnPrisms, + std::vector& dissipationHeatOnPrisms, std::vector& velocityOnVertices, int& error, const double& deltat = 0.0); @@ -201,7 +201,7 @@ double signedTriangleArea(const double* x, const double* y, const double* z); void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id); -void import2DFields(std::map bdExtensionMap, double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F, +void importFields(std::map bdExtensionMap, double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F, double const* beta_F = 0, double const* stiffnessFactor_F = 0, double const* effecPress_F = 0, double const* temperature_F = 0, double const* smb_F = 0, double eps = 0); void import2DFieldsObservations(std::map bdExtensionMap, @@ -220,8 +220,6 @@ void write_ascii_mesh_field_int(std::vector fieldData, std::string filename std::vector extendMaskByOneLayer(int const* verticesMask_F); -void importP0Temperature(); - void exportDissipationHeat(double * dissipationHeat_F); void get_prism_velocity_on_FEdges(double* uNormal, @@ -230,14 +228,9 @@ void get_prism_velocity_on_FEdges(double* uNormal, int initialize_iceProblem(int nTriangles); -void createReverseVerticesExchangeLists(exchangeList_Type& sendListReverse_F, - exchangeList_Type& receiveListReverse_F, - const std::vector& fVertexToTriangleID, const int* indexToVertexID_F); - -void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, +void createReverseExchangeLists(exchangeList_Type& sendListReverse_F, exchangeList_Type& receiveListReverse_F, - const std::vector& fCellToVertexID, - const int* indexToCellID_F); + const std::vector& newProcIds, const int* indexToID_F, exchangeList_Type const * recvList_F); void mapCellsToVertices(const std::vector& velocityOnCells, std::vector& velocityOnVertices, int fieldDim, int numLayers, From e5f9ac25de87e8f722350f15419f7a3a0355fc5b Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Thu, 14 Nov 2019 10:24:22 -0700 Subject: [PATCH 9/9] Modify albany input files of some tests to use Prisms and remove old options --- testing_and_setup/compass/landice/MISMIP+/albany_input.yaml | 2 +- .../compass/landice/Thwaites_variability/albany_input.yaml | 4 ++-- .../compass/landice/circular-shelf/albany_input.yaml | 2 -- .../compass/landice/confined-shelf/albany_input.yaml | 2 -- testing_and_setup/compass/landice/dome/albany_input.yaml | 1 + testing_and_setup/compass/landice/greenland/albany_input.yaml | 3 +-- .../compass/landice/greenland/albany_schoof_input.yaml | 2 -- 7 files changed, 5 insertions(+), 11 deletions(-) diff --git a/testing_and_setup/compass/landice/MISMIP+/albany_input.yaml b/testing_and_setup/compass/landice/MISMIP+/albany_input.yaml index c89dcb14ce..9676e6bbb9 100644 --- a/testing_and_setup/compass/landice/MISMIP+/albany_input.yaml +++ b/testing_and_setup/compass/landice/MISMIP+/albany_input.yaml @@ -35,7 +35,7 @@ ANONYMOUS: Body Force: Type: FO INTERP SURF GRAD Discretization: - Cubature Degree: 1 + Element Shape: Prism Piro: LOCA: Bifurcation: { } diff --git a/testing_and_setup/compass/landice/Thwaites_variability/albany_input.yaml b/testing_and_setup/compass/landice/Thwaites_variability/albany_input.yaml index af2a6dd48c..a41f919c8d 100644 --- a/testing_and_setup/compass/landice/Thwaites_variability/albany_input.yaml +++ b/testing_and_setup/compass/landice/Thwaites_variability/albany_input.yaml @@ -31,8 +31,8 @@ ANONYMOUS: Type: Lateral Immersed Ratio: 8.85000000000000009e-01 Discretization: - Cubature Degree: 11 - Interleaved Ordering: true + Element Shape: Tetrahedron + #Exodus Output File Name: albany_output.exo Piro: LOCA: Bifurcation: { } diff --git a/testing_and_setup/compass/landice/circular-shelf/albany_input.yaml b/testing_and_setup/compass/landice/circular-shelf/albany_input.yaml index 029369ad0c..3ac436fb05 100644 --- a/testing_and_setup/compass/landice/circular-shelf/albany_input.yaml +++ b/testing_and_setup/compass/landice/circular-shelf/albany_input.yaml @@ -17,8 +17,6 @@ ANONYMOUS: Type: FO INTERP SURF GRAD Discretization: Exodus Output File Name: albany_output.exo - Cubature Degree: 1 - Interleaved Ordering: false Piro: LOCA: Bifurcation: { } diff --git a/testing_and_setup/compass/landice/confined-shelf/albany_input.yaml b/testing_and_setup/compass/landice/confined-shelf/albany_input.yaml index 52809be2bf..1c75c642b0 100644 --- a/testing_and_setup/compass/landice/confined-shelf/albany_input.yaml +++ b/testing_and_setup/compass/landice/confined-shelf/albany_input.yaml @@ -17,8 +17,6 @@ ANONYMOUS: Type: FO INTERP SURF GRAD Discretization: Exodus Output File Name: albany_output.exo - Cubature Degree: 1 - Interleaved Ordering: false Piro: LOCA: Bifurcation: { } diff --git a/testing_and_setup/compass/landice/dome/albany_input.yaml b/testing_and_setup/compass/landice/dome/albany_input.yaml index 01ec65683a..3357ddccc3 100644 --- a/testing_and_setup/compass/landice/dome/albany_input.yaml +++ b/testing_and_setup/compass/landice/dome/albany_input.yaml @@ -11,6 +11,7 @@ ANONYMOUS: 'Glen''s Law Homotopy Parameter': 1.0 Discretization: Exodus Output File Name: albany_output.exo + Element Shape: Prism Piro: LOCA: Bifurcation: { } diff --git a/testing_and_setup/compass/landice/greenland/albany_input.yaml b/testing_and_setup/compass/landice/greenland/albany_input.yaml index d2fb7eab2a..67b76efabe 100644 --- a/testing_and_setup/compass/landice/greenland/albany_input.yaml +++ b/testing_and_setup/compass/landice/greenland/albany_input.yaml @@ -17,8 +17,7 @@ ANONYMOUS: Type: FO INTERP SURF GRAD Discretization: Exodus Output File Name: albany_output.exo - Cubature Degree: 1 - Interleaved Ordering: false + Element Shape: Tetrahedron Piro: LOCA: Bifurcation: { } diff --git a/testing_and_setup/compass/landice/greenland/albany_schoof_input.yaml b/testing_and_setup/compass/landice/greenland/albany_schoof_input.yaml index 3e1fa88dc5..045feb0c49 100644 --- a/testing_and_setup/compass/landice/greenland/albany_schoof_input.yaml +++ b/testing_and_setup/compass/landice/greenland/albany_schoof_input.yaml @@ -39,8 +39,6 @@ ANONYMOUS: Immersed Ratio: 8.8521e-01 Discretization: Exodus Output File Name: albany_output.exo - Cubature Degree: 1 - Interleaved Ordering: false Piro: LOCA: Bifurcation: { }