#include "../MEF90/mef90.inc" module MEF90_APPEND(m_MEF90_HeatXferAssembly,MEF90_DIM)D #include "petsc/finclude/petsc.h" use petscsnes use petsctao use m_MEF90_HeatXfer_class implicit none(type, external) private public MEF90HeatXferOperator, & MEF90HeatXferBilinearForm, & MEF90HeatXFerEnergy, & MEF90HeatXFerIFunction, & MEF90HeatXferIJacobian contains #undef __FUNCT__ #define __FUNCT__ "MEF90HeatXferOperator" !!! author: Blaise Bourdin (2012-14, bourdin@lsu.edu) !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90HeatXferOperator: Build the operator. When called in SNES, the solution time should always match the target time, !!! so there is no need for interpolation of the fluxes, external, and boundary values !!! subroutine MEF90HeatXferOperator(snesTemp, x, residual, MEF90HeatXferCtx, ierr) type(tSNES), intent(IN) :: snesTemp type(tVec), intent(IN) :: x type(tVec), intent(INOUT) :: residual type(MEF90HeatXfer_Type), intent(IN) :: MEF90HeatXferCtx PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmTemperature, dmFlux, dmBoundaryFlux, dmExternalTemperature type(tPetscSection) :: sectionFlux, sectionBoundaryFlux, sectionExternalTemperature type(tVec) :: locResidual type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder, vecOffset PetscInt :: cell, iDof, jDof, iGauss type(MEF90HeatXferCellSetOptions_Type) :: cellSetOptions type(MEF90_MATS) :: thermalConductivity type(MEF90HeatXferFaceSetOptions_Type) :: faceSetOptions character(len=MEF90MXSTRLEN) :: setPrefix type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elem type(MEF90ElementType) :: elementType DMPolytopeType :: cellType type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90HeatXferGlobalOptions_Type) :: MEF90HeatXferGlobalOptions PetscReal, dimension(:), pointer :: temperatureDof, fluxArray, boundaryFluxArray, externalTemperatureArray, residualDof type(MEF90_VECT) :: advectionVec PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90HeatXferCtx%MEF90Ctx%comm, trim(MEF90HeatXferCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) MEF90HeatXferGlobalOptions = MEF90HeatXferCtx%globalOptions PetscCall(VecGetDM(MEF90HeatXferCtx%temperatureLocal, dmTemperature, ierr)) PetscCall(VecGetDM(MEF90HeatXferCtx%fluxLocal, dmFlux, ierr)) PetscCall(DMGetLocalSection(dmFlux, sectionFlux, ierr)) PetscCall(VecGetDM(MEF90HeatXferCtx%boundaryFluxLocal, dmBoundaryFlux, ierr)) PetscCall(DMGetLocalSection(dmBoundaryFlux, sectionBoundaryFlux, ierr)) PetscCall(VecGetDM(MEF90HeatXferCtx%externalTemperatureLocal, dmExternalTemperature, ierr)) PetscCall(DMGetLocalSection(dmExternalTemperature, sectionExternalTemperature, ierr)) PetscCall(DMGetLocalVector(dmTemperature, locResidual, ierr)) PetscCall(DMGlobalToLocal(dmTemperature, x, INSERT_VALUES, MEF90HeatXferCtx%TemperatureLocal, ierr)) PetscCall(VecSet(residual, 0.0_kr, ierr)) PetscCall(VecSet(locResidual, 0.0_kr, ierr)) !! cell-based contributions PetscCall(DMGetLabelIdIS(dmTemperature, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90HeatXferCtx%MEF90Ctx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90HeatXferCtx%fluxLocal, fluxArray, ierr)) do set = 1, size(setID) PetscCall(DMGetStratumIS(dmTemperature, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (setPrefix, '(A,"cs",I4.4,"_")') trim(MEF90HeatXferCtx%prefix), setID(set) PetscCall(MEF90HeatXferCellSetOptionsSetFromOptions(MEF90HeatXferCtx%comm, setPrefix, MEF90_DIM, cellSetOptions, ierr)) select type (thermalConductivityMat => cellSetOptions%thermalConductivity) type is (MEF90_MATS) thermalConductivity = thermalConductivityMat end select PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmTemperature, setPointID(1), cellType, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellType, elementType, ierr)) QuadratureOrder = elementType%order * 2 PetscCall(MEF90ElementCreate(dmTemperature, setPointIS, elem, QuadratureOrder, elementType, ierr)) allocate (residualDof(size(elem(1)%BF(:, 1)))) do cell = 1, size(setPointID) residualDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, MEF90HeatXferCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 1, size(elem(cell)%BF(:, 1)) do iDof = 1, size(elem(cell)%BF(:, 1)) residualDof(jDof) = residualDof(jDof) + (thermalConductivity * temperatureDof(iDof) * elem(cell)%Grad_BF(iDof, iGauss) .dotP.elem(cell)%Grad_BF(jDof, iGauss)) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! jDof end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, MEF90HeatXferCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(DMPlexVecSetClosure(dmTemperature, PETSC_NULL_SECTION, locResidual, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell if (norm2(cellSetOptions%advectionVector) /= 0.0_kr) then advectionVec = cellSetOptions%advectionVector do cell = 1, size(setPointID) residualDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, MEF90HeatXferCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 1, size(elem(cell)%BF(:, 1)) do iDof = 1, size(elem(cell)%BF(:, 1)) residualDof(jDof) = residualDof(jDof) - (cellSetOptions%density * cellSetOptions%specificHeat * advectionVec.dotP.temperatureDof(iDof) * elem(cell)%Grad_BF(iDof, iGauss)) * elem(cell)%BF(jDof, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! jDof end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, MEF90HeatXferCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(DMPlexVecSetClosure(dmTemperature, PETSC_NULL_SECTION, locResidual, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell end if ! cellSetOptions%advectionVector if (cellSetOptions%flux /= 0.0_kr) then do cell = 1, size(setPointID) residualDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(PetscSectionGetOffset(sectionFlux, setPointID(cell), vecOffset, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) do iDof = 1, size(elem(cell)%BF(:, 1)) residualDof(iDof) = residualDof(iDof) - fluxArray(vecOffset + 1) * elem(cell)%BF(iDof, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! iGauss PetscCall(DMPlexVecSetClosure(dmTemperature, PETSC_NULL_SECTION, locResidual, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell end if ! cellSetOptions%flux deallocate (residualDof) PetscCall(MEF90ElementDestroy(elem, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! pointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90HeatXferCtx%fluxLocal, fluxArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS !! face-based contributions PetscCall(DMGetLabelIdIS(dmTemperature, MEF90FaceSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90HeatXferCtx%MEF90Ctx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90HeatXferCtx%boundaryFluxLocal, boundaryFluxArray, ierr)) PetscCall(VecGetArray(MEF90HeatXferCtx%externalTemperatureLocal, externalTemperatureArray, ierr)) do set = 1, size(setID) PetscCall(DMGetStratumIS(dmTemperature, MEF90FaceSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (setPrefix, '(A,"fs",I4.4,"_")') trim(MEF90HeatXferCtx%prefix), setID(set) PetscCall(MEF90HeatXferFaceSetOptionsSetFromOptions(MEF90HeatXferCtx%comm, setPrefix, faceSetOptions, ierr)) if (faceSetOptions%boundaryFlux /= 0.0_kr) then PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmTemperature, setPointID(1), cellType, ierr)) PetscCall(MEF90ElementGetTypeBoundary(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellType, elementType, ierr)) QuadratureOrder = elementType%order * 2 PetscCall(MEF90ElementCreate(dmTemperature, setPointIS, elem, QuadratureOrder, elementType, ierr)) allocate (residualDof(size(elem(1)%BF(:, 1)))) do cell = 1, size(setPointID) residualDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(PetscSectionGetOffset(sectionBoundaryFlux, setPointID(cell), vecOffset, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) do iDof = 1, size(elem(cell)%BF(:, 1)) residualDof(iDof) = residualDof(iDof) - boundaryFluxArray(vecOffset + 1) * elem(cell)%BF(iDof, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! iGauss PetscCall(DMPlexVecSetClosure(dmTemperature, PETSC_NULL_SECTION, locResidual, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell deallocate (residualDof) PetscCall(MEF90ElementDestroy(elem, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) end if ! faceSetOptions%boundaryFlux if (faceSetOptions%surfaceThermalConductivity /= 0.0_kr) then PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmTemperature, setPointID(1), cellType, ierr)) PetscCall(MEF90ElementGetTypeBoundary(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellType, elementType, ierr)) QuadratureOrder = elementType%order * 2 PetscCall(MEF90ElementCreate(dmTemperature, setPointIS, elem, QuadratureOrder, elementType, ierr)) allocate (residualDof(size(elem(1)%BF(:, 1)))) do cell = 1, size(setPointID) residualDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, MEF90HeatXferCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(PetscSectionGetOffset(sectionExternalTemperature, setPointID(cell), vecOffset, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 1, size(elem(cell)%BF(:, 1)) do iDof = 1, size(elem(cell)%BF(:, 1)) residualDof(jDof) = residualDof(jDof) + faceSetOptions%surfaceThermalConductivity * temperatureDof(iDof) * elem(cell)%BF(iDof, iGauss) * elem(cell)%BF(jDof, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof residualDof(jDof) = residualDof(jDof) - faceSetOptions%surfaceThermalConductivity * externalTemperatureArray(vecOffset + 1) * elem(cell)%BF(jDof, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! jDof end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, MEF90HeatXferCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(DMPlexVecSetClosure(dmTemperature, PETSC_NULL_SECTION, locResidual, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell deallocate (residualDof) PetscCall(MEF90ElementDestroy(elem, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) end if ! faceSetOptions%surfaceThermalConductivity PetscCall(ISDestroy(setPointIS, ierr)) end if ! pointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90HeatXferCtx%boundaryFluxLocal, boundaryFluxArray, ierr)) PetscCall(VecRestoreArray(MEF90HeatXferCtx%externalTemperatureLocal, externalTemperatureArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS PetscCall(DMLocalToGlobalBegin(dmTemperature, locResidual, ADD_VALUES, residual, ierr)) PetscCall(DMLocalToGlobalEnd(dmTemperature, locResidual, ADD_VALUES, residual, ierr)) !PetscCall(DMRestoreLocalVector(dmTemperature,locTemperature,ierr)) PetscCall(DMRestoreLocalVector(dmTemperature, locResidual, ierr)) end subroutine MEF90HeatXferOperator #undef __FUNCT__ #define __FUNCT__ "MEF90HeatXferBilinearForm" !!! author: Blaise Bourdin (2012-14, bourdin@lsu.edu) !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90HeatXferBilinearForm: !!! subroutine MEF90HeatXferBilinearForm(snesTemp, x, A, M, MEF90HeatXferCtx, ierr) type(tSNES), intent(IN) :: snesTemp type(tVec), intent(IN) :: x type(tMat), intent(INOUT) :: A, M type(MEF90HeatXfer_Type), intent(IN) :: MEF90HeatXferCtx PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmTemperature type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder PetscInt :: cell, iDof, jDof, iGauss, nbDof type(MEF90HeatXferCellSetOptions_Type) :: cellSetOptions type(MEF90_MATS) :: thermalConductivity type(MEF90HeatXferFaceSetOptions_Type) :: faceSetOptions character(len=MEF90MXSTRLEN) :: setPrefix type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elem type(MEF90ElementType) :: elementType DMPolytopeType :: cellType type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90HeatXferGlobalOptions_Type) :: MEF90HeatXferGlobalOptions PetscReal, dimension(:), pointer :: matDof type(MEF90_VECT) :: advectionVec PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90HeatXferCtx%MEF90Ctx%comm, trim(MEF90HeatXferCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) MEF90HeatXferGlobalOptions = MEF90HeatXferCtx%globalOptions PetscCall(VecGetDM(MEF90HeatXferCtx%temperatureLocal, dmTemperature, ierr)) PetscCall(MatZeroEntries(A, ierr)) !! cell-based gradient contributions PetscCall(DMGetLabelIdIS(dmTemperature, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90HeatXferCtx%MEF90Ctx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) PetscCall(DMGetStratumIS(dmTemperature, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (setPrefix, '(A,"cs",I4.4,"_")') trim(MEF90HeatXferCtx%prefix), setID(set) PetscCall(MEF90HeatXferCellSetOptionsSetFromOptions(MEF90HeatXferCtx%comm, setPrefix, MEF90_DIM, cellSetOptions, ierr)) select type (thermalConductivityMat => cellSetOptions%thermalConductivity) type is (MEF90_MATS) thermalConductivity = thermalConductivityMat end select PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmTemperature, setPointID(1), cellType, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellType, elementType, ierr)) QuadratureOrder = elementType%order * 2 PetscCall(MEF90ElementCreate(dmTemperature, setPointIS, elem, QuadratureOrder, elementType, ierr)) nbDof = size(elem(1)%BF(:, 1)) allocate (matDof(nbDof * nbDof)) do cell = 1, size(setPointID) matDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 0, nbDof - 1 do iDof = 1, nbDof matDof(jDof * nbDof + iDof) = matDof(jDof * nbDof + iDof) + (thermalConductivity * elem(cell)%Grad_BF(iDof, iGauss) .dotP.elem(cell)%Grad_BF(jDof + 1, iGauss)) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! jDof end do ! iGauss PetscCall(DMPlexMatSetClosure(dmTemperature, PETSC_NULL_SECTION, PETSC_NULL_SECTION, A, setPointID(cell), matDof, ADD_VALUES, ierr)) end do ! cell if (norm2(cellSetOptions%advectionVector) /= 0.0_kr) then advectionVec = cellSetOptions%advectionVector do cell = 1, size(setPointID) matDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 0, nbDof - 1 do iDof = 1, nbDof matDof(jDof * nbDof + iDof) = matDof(jDof * nbDof + iDof) - (cellSetOptions%density * cellSetOptions%specificHeat * advectionVec.dotP.elem(cell)%Grad_BF(iDof, iGauss)) * elem(cell)%BF(jDof + 1, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! jDof end do ! iGauss PetscCall(DMPlexMatSetClosure(dmTemperature, PETSC_NULL_SECTION, PETSC_NULL_SECTION, A, setPointID(cell), matDof, ADD_VALUES, ierr)) end do ! cell end if ! cellSetOptions%advectionVector deallocate (matDof) PetscCall(MEF90ElementDestroy(elem, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! pointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS !! face-based energies PetscCall(DMGetLabelIdIS(dmTemperature, MEF90FaceSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90HeatXferCtx%MEF90Ctx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) PetscCall(DMGetStratumIS(dmTemperature, MEF90FaceSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (setPrefix, '(A,"fs",I4.4,"_")') trim(MEF90HeatXferCtx%prefix), setID(set) PetscCall(MEF90HeatXferFaceSetOptionsSetFromOptions(MEF90HeatXferCtx%comm, setPrefix, faceSetOptions, ierr)) if (faceSetOptions%surfaceThermalConductivity /= 0.0_kr) then PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmTemperature, setPointID(1), cellType, ierr)) PetscCall(MEF90ElementGetTypeBoundary(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellType, elementType, ierr)) QuadratureOrder = elementType%order * 2 PetscCall(MEF90ElementCreate(dmTemperature, setPointIS, elem, QuadratureOrder, elementType, ierr)) nbDof = size(elem(1)%BF(:, 1)) allocate (matDof(nbDof * nbDof)) do cell = 1, size(setPointID) matDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 0, nbDof - 1 do iDof = 1, nbDof matDof(jDof * nbDof + iDof) = matDof(jDof * nbDof + iDof) + faceSetOptions%surfaceThermalConductivity * elem(cell)%BF(iDof, iGauss) * elem(cell)%BF(jDof + 1, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! jDof end do ! iGauss PetscCall(DMPlexMatSetClosure(dmTemperature, PETSC_NULL_SECTION, PETSC_NULL_SECTION, A, setPointID(cell), matDof, ADD_VALUES, ierr)) end do ! cell deallocate (matDof) PetscCall(MEF90ElementDestroy(elem, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) end if ! faceSetOptions%surfaceThermalConductivity PetscCall(ISDestroy(setPointIS, ierr)) end if ! pointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)) PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)) PetscCall(MatCopy(A, M, SAME_NONZERO_PATTERN, ierr)) end subroutine MEF90HeatXferBilinearForm #undef __FUNCT__ #define __FUNCT__ "MEF90HeatXFerEnergy" !!! author: Blaise Bourdin (2012-14, bourdin@lsu.edu) !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90HeatXFerEnergy: !!! subroutine MEF90HeatXFerEnergy(MEF90HeatXferCtx, energy, bodyWork, surfaceWork, ierr) type(MEF90HeatXfer_Type), intent(IN) :: MEF90HeatXferCtx PetscReal, dimension(:), pointer :: energy, bodyWork, surfaceWork PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmTemperature, dmFlux, dmBoundaryFlux type(tPetscSection) :: sectionFlux, sectionBoundaryFlux type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder, vecOffset PetscInt :: cell, iDof, iGauss type(MEF90HeatXferCellSetOptions_Type) :: cellSetOptions type(MEF90_MATS) :: thermalConductivity type(MEF90HeatXferFaceSetOptions_Type) :: faceSetOptions character(len=MEF90MXSTRLEN) :: setPrefix type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elem type(MEF90ElementType) :: elementType DMPolytopeType :: cellType type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90HeatXferGlobalOptions_Type) :: MEF90HeatXferGlobalOptions PetscReal, dimension(:), pointer :: TemperatureDof, fluxArray, boundaryFluxArray type(MEF90_VECT) :: gradTemperatureCell PetscReal :: bodyWorkCell, surfaceWorkCell PetscReal :: myEnergy, myBodyWork, mySurfaceWork energy = 0.0_kr bodyWork = 0.0_kr surfaceWork = 0.0_kr PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90HeatXferCtx%MEF90Ctx%comm, trim(MEF90HeatXferCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) MEF90HeatXferGlobalOptions = MEF90HeatXferCtx%globalOptions PetscCall(VecGetDM(MEF90HeatXferCtx%temperatureLocal, dmTemperature, ierr)) PetscCall(VecGetDM(MEF90HeatXferCtx%fluxLocal, dmFlux, ierr)) PetscCall(DMGetLocalSection(dmFlux, sectionFlux, ierr)) PetscCall(VecGetDM(MEF90HeatXferCtx%boundaryFluxLocal, dmBoundaryFlux, ierr)) PetscCall(DMGetLocalSection(dmBoundaryFlux, sectionBoundaryFlux, ierr)) !! cell-based energies PetscCall(DMGetLabelIdIS(dmTemperature, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90HeatXferCtx%MEF90Ctx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90HeatXferCtx%fluxLocal, fluxArray, ierr)) do set = 1, size(setID) !! myEnergy and myBodyWork are reduced over the whole communicator below, including over the ranks !! owning no point of this set, so they have to be initialized outside of the setPointIS test myEnergy = 0.0_kr myBodyWork = 0.0_kr PetscCall(DMGetStratumIS(dmTemperature, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) myenergy = 0.0_kr myBodyWork = 0.0_kr if (.not. PetscObjectIsNull(setPointIS)) then write (setPrefix, '(A,"cs",I4.4,"_")') trim(MEF90HeatXferCtx%prefix), setID(set) PetscCall(MEF90HeatXferCellSetOptionsSetFromOptions(MEF90HeatXferCtx%comm, setPrefix, MEF90_DIM, cellSetOptions, ierr)) select type (thermalConductivityMat => cellSetOptions%thermalConductivity) type is (MEF90_MATS) thermalConductivity = thermalConductivityMat end select PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmTemperature, setPointID(1), cellType, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellType, elementType, ierr)) QuadratureOrder = elementType%order * 2 PetscCall(MEF90ElementCreate(dmTemperature, setPointIS, elem, QuadratureOrder, elementType, ierr)) do cell = 1, size(setPointID) !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, MEF90HeatXferCtx%temperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) gradTemperatureCell = 0.0_kr do iDof = 1, size(elem(cell)%BF(:, 1)) gradTemperatureCell = gradTemperatureCell + temperatureDof(iDof) * elem(cell)%Grad_BF(iDof, iGauss) end do ! iDof myEnergy = myEnergy + ((thermalConductivity * gradTemperatureCell) .dotP.gradTemperatureCell) * elem(cell)%Gauss_C(iGauss) end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, MEF90HeatXferCtx%temperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) end do ! cell myEnergy = myEnergy * 0.5_kr if (cellSetOptions%flux /= 0.0_kr) then do cell = 1, size(setPointID) !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, MEF90HeatXferCtx%temperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(PetscSectionGetOffset(sectionFlux, setPointID(cell), vecOffset, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) bodyWorkCell = 0.0_kr do iDof = 1, size(elem(cell)%BF(:, 1)) bodyWorkCell = bodyWorkCell + fluxArray(vecOffset + 1) * temperatureDof(iDof) * elem(cell)%BF(iDof, iGauss) end do ! iDof myBodyWork = myBodyWork + bodyWorkCell * elem(cell)%Gauss_C(iGauss) end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, MEF90HeatXferCtx%temperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) end do ! cell end if ! cellSetOptions%flux PetscCall(MEF90ElementDestroy(elem, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! pointIS PetscCallMPI(MPI_AllReduce(myEnergy, energy(set), 1, MPIU_SCALAR, MPI_SUM, MEF90HeatXferCtx%MEF90Ctx%comm, ierr)) PetscCallMPI(MPI_AllReduce(myBodyWork, bodyWork(set), 1, MPIU_SCALAR, MPI_SUM, MEF90HeatXferCtx%MEF90Ctx%comm, ierr)) end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90HeatXferCtx%fluxLocal, fluxArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS !! face-based energies PetscCall(DMGetLabelIdIS(dmTemperature, MEF90FaceSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90HeatXferCtx%MEF90Ctx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90HeatXferCtx%boundaryFluxLocal, boundaryFluxArray, ierr)) do set = 1, size(setID) !! mySurfaceWork is reduced over the whole communicator below, including over the ranks !! owning no point of this set, so it has to be initialized outside of the setPointIS test mySurfaceWork = 0.0_kr PetscCall(DMGetStratumIS(dmTemperature, MEF90FaceSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (setPrefix, '(A,"fs",I4.4,"_")') trim(MEF90HeatXferCtx%prefix), setID(set) PetscCall(MEF90HeatXferFaceSetOptionsSetFromOptions(MEF90HeatXferCtx%comm, setPrefix, faceSetOptions, ierr)) if (faceSetOptions%boundaryFlux /= 0.0_kr) then PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmTemperature, setPointID(1), cellType, ierr)) PetscCall(MEF90ElementGetTypeBoundary(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellType, elementType, ierr)) QuadratureOrder = elementType%order * 2 PetscCall(MEF90ElementCreate(dmTemperature, setPointIS, elem, QuadratureOrder, elementType, ierr)) do cell = 1, size(setPointID) !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, MEF90HeatXferCtx%temperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(PetscSectionGetOffset(sectionBoundaryFlux, setPointID(cell), vecOffset, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) surfaceWorkCell = 0.0_kr do iDof = 1, size(elem(cell)%BF(:, 1)) surfaceWorkCell = surfaceWorkCell + boundaryFluxArray(vecOffset + 1) * temperatureDof(iDof) * elem(cell)%BF(iDof, iGauss) end do ! iDof mySurfaceWork = mySurfaceWork + surfaceWorkCell * elem(cell)%Gauss_C(iGauss) end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, MEF90HeatXferCtx%temperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) end do ! cell PetscCall(MEF90ElementDestroy(elem, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) end if ! faceSetOptions%boundaryFlux PetscCall(ISDestroy(setPointIS, ierr)) end if ! pointIS PetscCallMPI(MPI_AllReduce(mySurfaceWork, surfaceWork(set), 1, MPIU_SCALAR, MPI_SUM, MEF90HeatXferCtx%MEF90Ctx%comm, ierr)) end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90HeatXferCtx%boundaryFluxLocal, boundaryFluxArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS end subroutine MEF90HeatXFerEnergy #undef __FUNCT__ #define __FUNCT__ "MEF90HeatXFerIFunction" !!! author: Blaise Bourdin (2014, bourdin@lsu.edu) !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90HeatXFerIFunction: !!! subroutine MEF90HeatXFerIFunction(tempTS, time, x, xdot, F, MEF90HeatXferCtx, ierr) type(tTS), intent(IN) :: tempTS PetscReal, intent(IN) :: time type(tVec), intent(IN) :: x, xdot type(tVec), intent(INOUT) :: F type(MEF90HeatXfer_Type), intent(IN) :: MEF90HeatXferCtx PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmTemperature, dmFlux, dmBoundaryFlux, dmExternalTemperature type(tPetscSection) :: sectionFlux, sectionBoundaryFlux, sectionExternalTemperature type(tVec) :: locTemperature, locTemperatureDot, locF type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder, vecOffset PetscInt :: cell, iDof, jDof, iGauss type(MEF90HeatXferCellSetOptions_Type) :: cellSetOptions type(MEF90_MATS) :: thermalConductivity type(MEF90HeatXferFaceSetOptions_Type) :: faceSetOptions character(len=MEF90MXSTRLEN) :: setPrefix type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elem type(MEF90ElementType) :: elementType DMPolytopeType :: cellType type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90HeatXferGlobalOptions_Type) :: MEF90HeatXferGlobalOptions PetscReal, dimension(:), pointer :: temperatureDof, temperatureDotDof, fluxArray, boundaryFluxArray, externalTemperatureArray, residualDof type(MEF90_VECT) :: advectionVec PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90HeatXferCtx%MEF90Ctx%comm, trim(MEF90HeatXferCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) MEF90HeatXferGlobalOptions = MEF90HeatXferCtx%globalOptions PetscCall(VecGetDM(MEF90HeatXferCtx%temperatureLocal, dmTemperature, ierr)) PetscCall(VecGetDM(MEF90HeatXferCtx%fluxLocal, dmFlux, ierr)) PetscCall(DMGetLocalSection(dmFlux, sectionFlux, ierr)) PetscCall(VecGetDM(MEF90HeatXferCtx%boundaryFluxLocal, dmBoundaryFlux, ierr)) PetscCall(DMGetLocalSection(dmBoundaryFlux, sectionBoundaryFlux, ierr)) PetscCall(VecGetDM(MEF90HeatXferCtx%externalTemperatureLocal, dmExternalTemperature, ierr)) PetscCall(DMGetLocalSection(dmExternalTemperature, sectionExternalTemperature, ierr)) PetscCall(DMGetLocalVector(dmTemperature, locTemperature, ierr)) PetscCall(DMGetLocalVector(dmTemperature, locTemperatureDot, ierr)) PetscCall(DMGetLocalVector(dmTemperature, locF, ierr)) PetscCall(MEF90VecGlobalToLocalConstraint(x, MEF90HeatXferCtx%temperatureLocal, locTemperature, ierr)) PetscCall(VecSet(locTemperatureDot, 0.0_kr, ierr)) PetscCall(DMGlobalToLocalBegin(dmTemperature, xdot, INSERT_VALUES, locTemperatureDot, ierr)) PetscCall(DMGlobalToLocalEnd(dmTemperature, xdot, INSERT_VALUES, locTemperatureDot, ierr)) PetscCall(VecSet(F, 0.0_kr, ierr)) PetscCall(VecSet(locF, 0.0_kr, ierr)) !! cell-based contributions PetscCall(DMGetLabelIdIS(dmTemperature, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90HeatXferCtx%MEF90Ctx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90HeatXferCtx%fluxLocal, fluxArray, ierr)) do set = 1, size(setID) PetscCall(DMGetStratumIS(dmTemperature, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (setPrefix, '(A,"cs",I4.4,"_")') trim(MEF90HeatXferCtx%prefix), setID(set) PetscCall(MEF90HeatXferCellSetOptionsSetFromOptions(MEF90HeatXferCtx%comm, setPrefix, MEF90_DIM, cellSetOptions, ierr)) select type (thermalConductivityMat => cellSetOptions%thermalConductivity) type is (MEF90_MATS) thermalConductivity = thermalConductivityMat end select PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmTemperature, setPointID(1), cellType, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellType, elementType, ierr)) QuadratureOrder = elementType%order * 2 PetscCall(MEF90ElementCreate(dmTemperature, setPointIS, elem, QuadratureOrder, elementType, ierr)) allocate (residualDof(size(elem(1)%BF(:, 1)))) do cell = 1, size(setPointID) residualDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, locTemperatureDot, setPointID(cell), PETSC_NULL_INTEGER, temperatureDotDof, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 1, size(elem(cell)%BF(:, 1)) do iDof = 1, size(elem(cell)%BF(:, 1)) residualDof(jDof) = residualDof(jDof) + cellSetOptions%density * cellSetOptions%specificHeat * temperatureDotDof(iDof) * elem(cell)%BF(iDof, iGauss) * elem(cell)%BF(jDof, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! jDof end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, locTemperatureDot, setPointID(cell), PETSC_NULL_INTEGER, temperatureDotDof, ierr)) PetscCall(DMPlexVecSetClosure(dmTemperature, PETSC_NULL_SECTION, locF, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell do cell = 1, size(setPointID) residualDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, locTemperature, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 1, size(elem(cell)%BF(:, 1)) do iDof = 1, size(elem(cell)%BF(:, 1)) residualDof(jDof) = residualDof(jDof) + (thermalConductivity * temperatureDof(iDof) * elem(cell)%Grad_BF(iDof, iGauss) .dotP.elem(cell)%Grad_BF(jDof, iGauss)) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! jDof end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, locTemperature, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(DMPlexVecSetClosure(dmTemperature, PETSC_NULL_SECTION, locF, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell if (norm2(cellSetOptions%advectionVector) /= 0.0_kr) then advectionVec = cellSetOptions%advectionVector do cell = 1, size(setPointID) residualDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, locTemperature, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 1, size(elem(cell)%BF(:, 1)) do iDof = 1, size(elem(cell)%BF(:, 1)) residualDof(jDof) = residualDof(jDof) - (cellSetOptions%density * cellSetOptions%specificHeat * advectionVec.dotP.temperatureDof(iDof) * elem(cell)%Grad_BF(iDof, iGauss)) * elem(cell)%BF(jDof, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! jDof end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, locTemperature, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(DMPlexVecSetClosure(dmTemperature, PETSC_NULL_SECTION, locF, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell end if ! cellSetOptions%advectionVector if (cellSetOptions%flux /= 0.0_kr) then do cell = 1, size(setPointID) residualDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(PetscSectionGetOffset(sectionFlux, setPointID(cell), vecOffset, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) do iDof = 1, size(elem(cell)%BF(:, 1)) residualDof(iDof) = residualDof(iDof) - fluxArray(vecOffset + 1) * elem(cell)%BF(iDof, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! iGauss PetscCall(DMPlexVecSetClosure(dmTemperature, PETSC_NULL_SECTION, locF, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell end if ! cellSetOptions%flux deallocate (residualDof) PetscCall(MEF90ElementDestroy(elem, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! pointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90HeatXferCtx%fluxLocal, fluxArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS !! face-based contributions PetscCall(DMGetLabelIdIS(dmTemperature, MEF90FaceSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90HeatXferCtx%MEF90Ctx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90HeatXferCtx%boundaryFluxLocal, boundaryFluxArray, ierr)) PetscCall(VecGetArray(MEF90HeatXferCtx%externalTemperatureLocal, externalTemperatureArray, ierr)) do set = 1, size(setID) PetscCall(DMGetStratumIS(dmTemperature, MEF90FaceSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (setPrefix, '(A,"fs",I4.4,"_")') trim(MEF90HeatXferCtx%prefix), setID(set) PetscCall(MEF90HeatXferFaceSetOptionsSetFromOptions(MEF90HeatXferCtx%comm, setPrefix, faceSetOptions, ierr)) if (faceSetOptions%boundaryFlux /= 0.0_kr) then PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmTemperature, setPointID(1), cellType, ierr)) PetscCall(MEF90ElementGetTypeBoundary(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellType, elementType, ierr)) QuadratureOrder = elementType%order * 2 PetscCall(MEF90ElementCreate(dmTemperature, setPointIS, elem, QuadratureOrder, elementType, ierr)) allocate (residualDof(size(elem(1)%BF(:, 1)))) do cell = 1, size(setPointID) residualDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(PetscSectionGetOffset(sectionBoundaryFlux, setPointID(cell), vecOffset, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) do iDof = 1, size(elem(cell)%BF(:, 1)) residualDof(iDof) = residualDof(iDof) - boundaryFluxArray(vecOffset + 1) * elem(cell)%BF(iDof, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! iGauss PetscCall(DMPlexVecSetClosure(dmTemperature, PETSC_NULL_SECTION, locF, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell deallocate (residualDof) PetscCall(MEF90ElementDestroy(elem, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) end if ! faceSetOptions%boundaryFlux if (faceSetOptions%surfaceThermalConductivity /= 0.0_kr) then PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmTemperature, setPointID(1), cellType, ierr)) PetscCall(MEF90ElementGetTypeBoundary(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellType, elementType, ierr)) QuadratureOrder = elementType%order * 2 PetscCall(MEF90ElementCreate(dmTemperature, setPointIS, elem, QuadratureOrder, elementType, ierr)) allocate (residualDof(size(elem(1)%BF(:, 1)))) do cell = 1, size(setPointID) residualDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, locTemperature, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(PetscSectionGetOffset(sectionExternalTemperature, setPointID(cell), vecOffset, ierr)) do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 1, size(elem(cell)%BF(:, 1)) do iDof = 1, size(elem(cell)%BF(:, 1)) residualDof(jDof) = residualDof(jDof) + faceSetOptions%surfaceThermalConductivity * temperatureDof(iDof) * elem(cell)%BF(iDof, iGauss) * elem(cell)%BF(jDof, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof residualDof(jDof) = residualDof(jDof) - faceSetOptions%surfaceThermalConductivity * externalTemperatureArray(vecOffset + 1) * elem(cell)%BF(jDof, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! jDof end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, locTemperature, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(DMPlexVecSetClosure(dmTemperature, PETSC_NULL_SECTION, locF, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell deallocate (residualDof) PetscCall(MEF90ElementDestroy(elem, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) end if ! faceSetOptions%surfaceThermalConductivity PetscCall(ISDestroy(setPointIS, ierr)) end if ! pointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90HeatXferCtx%boundaryFluxLocal, boundaryFluxArray, ierr)) PetscCall(VecRestoreArray(MEF90HeatXferCtx%externalTemperatureLocal, externalTemperatureArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS PetscCall(DMLocalToGlobalBegin(dmTemperature, locF, ADD_VALUES, F, ierr)) PetscCall(DMLocalToGlobalEnd(dmTemperature, locF, ADD_VALUES, F, ierr)) PetscCall(DMRestoreLocalVector(dmTemperature, locTemperature, ierr)) PetscCall(DMRestoreLocalVector(dmTemperature, locTemperatureDot, ierr)) PetscCall(DMRestoreLocalVector(dmTemperature, locF, ierr)) end subroutine MEF90HeatXFerIFunction #undef __FUNCT__ #define __FUNCT__ "MEF90HeatXferIJacobian" !!! author: Blaise Bourdin (2014, bourdin@lsu.edu) !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90HeatXferIJacobian: !!! subroutine MEF90HeatXferIJacobian(tempTS, t, x, xdot, shift, A, M, MEF90HeatXferCtx, ierr) type(tTS), intent(IN) :: tempTS PetscReal, intent(IN) :: t type(tVec), intent(IN) :: x, xdot PetscReal, intent(IN) :: shift type(tMat), intent(INOUT) :: A, M type(MEF90HeatXfer_Type), intent(IN) :: MEF90HeatXferCtx PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmTemperature type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder PetscInt :: cell, iDof, jDof, iGauss, nbDof type(MEF90HeatXferCellSetOptions_Type) :: cellSetOptions type(MEF90_MATS) :: thermalConductivity type(MEF90HeatXferFaceSetOptions_Type) :: faceSetOptions character(len=MEF90MXSTRLEN) :: setPrefix type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elem type(MEF90ElementType) :: elementType DMPolytopeType :: cellType type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90HeatXferGlobalOptions_Type) :: MEF90HeatXferGlobalOptions PetscReal, dimension(:), pointer :: matDof type(MEF90_VECT) :: advectionVec PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90HeatXferCtx%MEF90Ctx%comm, trim(MEF90HeatXferCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) MEF90HeatXferGlobalOptions = MEF90HeatXferCtx%globalOptions PetscCall(VecGetDM(MEF90HeatXferCtx%temperatureLocal, dmTemperature, ierr)) PetscCall(MatZeroEntries(A, ierr)) !! cell-based gradient contributions PetscCall(DMGetLabelIdIS(dmTemperature, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90HeatXferCtx%MEF90Ctx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) PetscCall(DMGetStratumIS(dmTemperature, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (setPrefix, '(A,"cs",I4.4,"_")') trim(MEF90HeatXferCtx%prefix), setID(set) PetscCall(MEF90HeatXferCellSetOptionsSetFromOptions(MEF90HeatXferCtx%comm, setPrefix, MEF90_DIM, cellSetOptions, ierr)) select type (thermalConductivityMat => cellSetOptions%thermalConductivity) type is (MEF90_MATS) thermalConductivity = thermalConductivityMat end select PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmTemperature, setPointID(1), cellType, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellType, elementType, ierr)) QuadratureOrder = elementType%order * 2 PetscCall(MEF90ElementCreate(dmTemperature, setPointIS, elem, QuadratureOrder, elementType, ierr)) nbDof = size(elem(1)%BF(:, 1)) allocate (matDof(nbDof * nbDof)) do cell = 1, size(setPointID) matDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 0, nbDof - 1 do iDof = 1, nbDof matDof(jDof * nbDof + iDof) = matDof(jDof * nbDof + iDof) + shift * cellSetOptions%density * cellSetOptions%specificHeat * elem(cell)%BF(iDof, iGauss) * elem(cell)%BF(jDof + 1, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! jDof end do ! iGauss PetscCall(DMPlexMatSetClosure(dmTemperature, PETSC_NULL_SECTION, PETSC_NULL_SECTION, A, setPointID(cell), matDof, ADD_VALUES, ierr)) end do ! cell do cell = 1, size(setPointID) matDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 0, nbDof - 1 do iDof = 1, nbDof matDof(jDof * nbDof + iDof) = matDof(jDof * nbDof + iDof) + (thermalConductivity * elem(cell)%Grad_BF(iDof, iGauss) .dotP.elem(cell)%Grad_BF(jDof + 1, iGauss)) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! jDof end do ! iGauss PetscCall(DMPlexMatSetClosure(dmTemperature, PETSC_NULL_SECTION, PETSC_NULL_SECTION, A, setPointID(cell), matDof, ADD_VALUES, ierr)) end do ! cell if (norm2(cellSetOptions%advectionVector) /= 0.0_kr) then advectionVec = cellSetOptions%advectionVector do cell = 1, size(setPointID) matDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 0, nbDof - 1 do iDof = 1, nbDof matDof(jDof * nbDof + iDof) = matDof(jDof * nbDof + iDof) - (cellSetOptions%density * cellSetOptions%specificHeat * advectionVec.dotP.elem(cell)%Grad_BF(iDof, iGauss)) * elem(cell)%BF(jDof + 1, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! jDof end do ! iGauss PetscCall(DMPlexMatSetClosure(dmTemperature, PETSC_NULL_SECTION, PETSC_NULL_SECTION, A, setPointID(cell), matDof, ADD_VALUES, ierr)) end do ! cell end if ! cellSetOptions%advectionVector deallocate (matDof) PetscCall(MEF90ElementDestroy(elem, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! pointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS !! face-based energies PetscCall(DMGetLabelIdIS(dmTemperature, MEF90FaceSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90HeatXferCtx%MEF90Ctx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) PetscCall(DMGetStratumIS(dmTemperature, MEF90FaceSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (setPrefix, '(A,"fs",I4.4,"_")') trim(MEF90HeatXferCtx%prefix), setID(set) PetscCall(MEF90HeatXferFaceSetOptionsSetFromOptions(MEF90HeatXferCtx%comm, setPrefix, faceSetOptions, ierr)) if (faceSetOptions%surfaceThermalConductivity /= 0.0_kr) then PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmTemperature, setPointID(1), cellType, ierr)) PetscCall(MEF90ElementGetTypeBoundary(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellType, elementType, ierr)) QuadratureOrder = elementType%order * 2 PetscCall(MEF90ElementCreate(dmTemperature, setPointIS, elem, QuadratureOrder, elementType, ierr)) nbDof = size(elem(1)%BF(:, 1)) allocate (matDof(nbDof * nbDof)) do cell = 1, size(setPointID) matDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop do iGauss = 1, size(elem(cell)%Gauss_C) do jDof = 0, nbDof - 1 do iDof = 1, nbDof matDof(jDof * nbDof + iDof) = matDof(jDof * nbDof + iDof) + faceSetOptions%surfaceThermalConductivity * elem(cell)%BF(iDof, iGauss) * elem(cell)%BF(jDof + 1, iGauss) * elem(cell)%Gauss_C(iGauss) end do ! iDof end do ! jDof end do ! iGauss PetscCall(DMPlexMatSetClosure(dmTemperature, PETSC_NULL_SECTION, PETSC_NULL_SECTION, A, setPointID(cell), matDof, ADD_VALUES, ierr)) end do ! cell deallocate (matDof) PetscCall(MEF90ElementDestroy(elem, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) end if ! faceSetOptions%surfaceThermalConductivity PetscCall(ISDestroy(setPointIS, ierr)) end if ! pointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)) PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)) PetscCall(MatCopy(A, M, SAME_NONZERO_PATTERN, ierr)) end subroutine MEF90HeatXferIJacobian end module MEF90_APPEND(m_MEF90_HeatXferAssembly,MEF90_DIM)D