diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index 9f11d44ce..274ae93d0 100644 --- a/adflow/pyADflow.py +++ b/adflow/pyADflow.py @@ -18,21 +18,23 @@ v. 1.0 - Original pyAero Framework Implementation (RP,SM 2008) """ +import copy +import hashlib # ============================================================================= # Imports # ============================================================================= import os +import sys import time -import copy import types +from collections import OrderedDict + import numpy -import sys +from baseclasses import AeroProblem, AeroSolver, getPy3SafeString +from baseclasses.utils import CaseInsensitiveDict, Error from mpi4py import MPI -from baseclasses import AeroSolver, AeroProblem, getPy3SafeString -from baseclasses.utils import Error, CaseInsensitiveDict + from . import MExt -import hashlib -from collections import OrderedDict class ADFLOWWarning(object): @@ -3063,6 +3065,35 @@ def writeSurfaceSensitivity(self, fileName, func, groupName=None): f.close() # end if (root proc ) + def writeFamilySolution(self, familyName, outputDir=None, baseName=None, number=None): + if outputDir is None: + outputDir = self.getOption("outputDirectory") + + if baseName is None: + baseName = self.curAP.name + + if not familyName.lower() in self.families: + raise Error(f"Family {familyName} not found in the solver.") + + famList = self._getFamilyList(familyName) + numDigits = self.getOption("writeSolutionDigits") + famID = self.families[familyName.lower()][0] + isBC = True if famID in self._getFamilyList("allSurfaces") else False + + if number is not None: + baseName = f"{baseName}_{familyName.lower()}_{number:0{numDigits}d}" + else: + if self.getOption("numberSolutions"): + baseName = f"{baseName}_{familyName.lower()}_{self.curAP.adflowData.callCounter:0{numDigits}d}" + + fileName = os.path.join(outputDir, baseName) + fileName += ".dat" + + if isBC: + self.adflow.tecplotio.writebcsurfacesascii(fileName, famList) + else: + self.adflow.tecplotio.writeuserintsurf(familyName, fileName, famID) + def resetAdjoint(self, obj): """ Reset an adjoint 'obj' that is stored in the current diff --git a/src/f2py/adflow.pyf b/src/f2py/adflow.pyf index 4eec9145d..02775fe66 100644 --- a/src/f2py/adflow.pyf +++ b/src/f2py/adflow.pyf @@ -654,8 +654,41 @@ python module libadflow integer(kind=inttype), optional,intent(in),check(len(famlist)>=nfamlist),depend(famlist) :: nfamlist=len(famlist) end subroutine writetecplot + subroutine writebcsurfacesascii(filename,famlist) ! in output/tecplotIO.F90:tecplotio + use constants + use communication, only: myid,adflow_comm_world,nproc + use inputtimespectral, only: ntimeintervalsspectral + use inputphysics, only: equationmode + use inputiteration, only: printiterations + use inputio, only: precisionsurfgrid,precisionsurfsol + use outputmod, only: surfsolnames,numberofsurfsolvariables + use surfacefamilies, only: bcfamexchange,famnames,familyexchange + use utils, only: echk,setpointers,setbcpointers + use bcpointers, only: xx + use sorting, only: faminlist + use extraoutput, only: surfwriteblank + use oversetdata, only: zippermesh,zippermeshes + use surfaceutils + use petsc + character*(*) intent(in) :: filename + integer(kind=inttype) dimension(:),intent(in) :: famlist + end subroutine writebcsurfacesascii + subroutine initializeliftdistributiondata end subroutine initializeliftdistributiondata + + subroutine writeuserintsurf(slicename,filename,famid) ! in :test:tecplotIO.F90:tecplotio + use constants + use inputio + use commonformats, only: int5,sci6 + use usersurfaceintegrationdata, only: userintsurfs,nuserintsurfs,userintsurf + use usersurfaceintegrations, only: getusersurfacedata + use sorting, only: faminlist + character*(*) intent(in) :: slicename + character*(*) intent(in) :: filename + integer(kind=inttype) intent(in) :: famid + end subroutine writeuserintsurf + end subroutine tecplotio module surfaceintegrations @@ -681,7 +714,6 @@ python module libadflow subroutine interpolateintegrationsurfaces end subroutine interpolateintegrationsurfaces - end module usersurfaceintegrations module actuatorregion ! in :test:actuatorDiskRegion.F90 diff --git a/src/modules/commonFormats.F90 b/src/modules/commonFormats.F90 index ab6fd8bc0..b5f3eee42 100644 --- a/src/modules/commonFormats.F90 +++ b/src/modules/commonFormats.F90 @@ -35,4 +35,7 @@ module commonFormats ! Integers written with 5 characters character(len=maxStringLen) :: int5 = '(*(I5))' + ! Integers written with 7 characters + character(len=maxStringLen) :: int7 = '(*(I7))' + end module commonFormats diff --git a/src/output/tecplotIO.F90 b/src/output/tecplotIO.F90 index 1b5454e91..d433e4cd4 100644 --- a/src/output/tecplotIO.F90 +++ b/src/output/tecplotIO.F90 @@ -62,6 +62,7 @@ module tecplotIO integer(kind=intType), parameter :: nSliceMax = 1000 integer(kind=intType) :: nParaSlices = 0 integer(kind=intType) :: nAbsSlices = 0 + integer(kind=intType) :: nVolSlice = 0 type(slice), dimension(:, :), allocatable :: paraSlices, absSlices ! Data for the user supplied lift distributions @@ -2311,4 +2312,457 @@ subroutine writeSlice(slc, fileID, nFields) end if end subroutine writeSlice + subroutine writeUserIntSurf(sliceName, fileName, famID) + use constants + use inputIO + use commonFormats, only: int7 + use userSurfaceIntegrationData, only: userIntSurfs, nUserIntSurfs, userIntSurf + use userSurfaceIntegrations, only: getUserSurfaceData + use sorting, only: famInList + use communication, only: myID + use inputTimeSpectral, only: nTimeIntervalsSpectral + use inputIteration, only: printIterations + implicit none + + ! Input parameters + character(len=*), intent(in) :: sliceName + character(len=*), intent(in) :: fileName + integer(kind=intType), intent(in) :: famID + + ! Working variables + integer(kind=intType) :: i, j, iSurf, sps, file + real(kind=realType) :: tmp + real(kind=realType), dimension(:, :), allocatable :: vars + integer(kind=intType), dimension(:, :), allocatable :: conn + type(userIntSurf) :: surf + + if (myID == 0 .and. printIterations) then + print "(a)", "#" + print "(a)", "Writing volume slice" + end if + + timeLoop: do sps = 1, nTimeIntervalsSpectral + + ! This needs to happen on all processors because + ! getUserSurfaceData has to communicate the data for this surface. + ! Only the vars and conn on the root proc will have data. + do iSurf = 1, nUserIntSurfs + surf = userIntSurfs(iSurf) + if (surf%famID == famID) then + call getUserSurfaceData(iSurf, sps, vars, conn) + exit + end if + end do + + file = 11 + if (myID == 0) then + print "(a,4x,a)", "#", trim(fileName) + open (unit=file, file=trim(fileName)) + + ! Write the header information + write (file, *) "Title = ""ADflow Volume Slice""" + write (file, "(a)", advance='no') "Variables = " + write (file, "(a)", advance="no") " ""CoordinateX"" " + write (file, "(a)", advance="no") " ""CoordinateY"" " + write (file, "(a)", advance="no") " ""CoordinateZ"" " + write (file, "(a)", advance="no") " ""VelocityX"" " + write (file, "(a)", advance="no") " ""VelocityY"" " + write (file, "(a)", advance="no") " ""VelocityZ"" " + write (file, "(a)", advance="no") " ""Density"" " + write (file, "(a)", advance="no") " ""StaticPressure"" " + write (file, "(a)", advance="no") " ""TotalPressure"" " + write (file, "(a)", advance="no") " ""TotalTemperature"" " + write (file, "(a)", advance="no") " ""Mach"" " + + write (file, "(1x)") + + write (file, "(a,a,a)") "Zone T= """, trim(sliceName), """" + + if (size(vars, 2) > 0) then + write (file, *) "Nodes = ", size(vars, 2), " Elements= ", size(conn, 2), " ZONETYPE=FETRIANGLE" + write (file, *) "DATAPACKING=POINT" + + do i = 1, size(vars, 2) + ! Write the variables + do j = 1, size(vars, 1) + write (file, sci6, advance='no') vars(j, i) + end do + + write (file, "(1x)") + end do + + do i = 1, size(conn, 2) + write (file, int7) conn(1, i), conn(2, i), conn(3, i) + end do + else + write (file, *) "Nodes = ", 3, " Elements= ", 1, " ZONETYPE=FETRIANGLE" + write (file, *) "DATAPACKING=POINT" + do i = 1, 3 + do j = 1, 10 + write (file, sci6, advance='no') zero + end do + + write (file, "(1x)") + end do + write (file, int7) 1, 2, 3 + end if + end if + end do timeLoop + + end subroutine writeUserIntSurf + + subroutine writeBCSurfacesASCII(fileName, famList) + use constants + use communication, only: myid, adflow_comm_world, nProc + use commonFormats, only: int7 + use inputTimeSpectral, only: nTimeIntervalsSpectral + use inputPhysics, only: equationMode + use inputIteration, only: printIterations + use inputIO, only: precisionsurfgrid, precisionsurfsol + use outputMod, only: surfSolNames, numberOfSurfSolVariables + use surfaceFamilies, onlY: BCFamExchange, famNames, familyExchange + use utils, only: EChk, setPointers, setBCPointers + use BCPointers, only: xx + use sorting, only: famInList + use extraOutput, only: surfWriteBlank + use oversetData, only: zipperMesh, zipperMeshes + use surfaceUtils +#include + use petsc + implicit none + + ! Input Params + character(len=*), intent(in) :: fileName + integer(kind=intType), intent(in), dimension(:) :: famList + + ! Working parameters + integer(kind=intType) :: i, j, k, nn, mm, fileID, iVar, ii, ierr, iSize + integer(Kind=intType) :: nSolVar, iBeg, iEnd, jBeg, jEnd, sps, sizeNode, sizeCell + integer(kind=intType) :: iBCGroup, iFam, iProc, nCells, nNodes, nCellsToWrite, iZone, lastZoneSharing + character(len=maxStringLen) :: fname + character(len=7) :: intString + integer(kind=intType), dimension(:), allocatable :: nodeSizes, nodeDisps + integer(kind=intType), dimension(:), allocatable :: cellSizes, cellDisps + character(len=maxCGNSNameLen), dimension(:), allocatable :: solNames + real(kind=realType), dimension(:, :), allocatable :: nodalValues + integer(kind=intType), dimension(:, :), allocatable :: conn, localConn + real(kind=realType), dimension(:, :), allocatable :: vars + integer(kind=intType), dimension(:), allocatable :: mask, elemFam, localElemFam, cgnsBlockID + logical :: blankSave, BCGroupNeeded, dataWritten + type(zipperMesh), pointer :: zipper + type(familyExchange), pointer :: exch + if (myID == 0 .and. printIterations) then + print "(a)", "#" + print "(a)", "# Writing tecplot surface file(s):" + end if + + ! Number of surface variables. Note that we *explictly* + ! remove the potential for writing the surface blanks as + ! these are not necessary for the tecplot IO as we write the + ! zipper mesh directly. We must save and restore the + ! variable in case the CGNS otuput still wants to write it. + blankSave = surfWriteBlank + surfWriteBlank = .False. + call numberOfSurfSolVariables(nSolVar) + allocate (solNames(nSolVar)) + call surfSolNames(solNames) + + spectralLoop: do sps = 1, nTimeIntervalsSpectral + + ! If it is time spectral we need to agument the filename + if (equationMode == timeSpectral) then + write (intString, "(i7)") sps + intString = adjustl(intString) + fname = trim(fileName)//"Spectral"//trim(intString) + else + fname = fileName + end if + + fileID = 11 + ! Open file on root proc: + if (myid == 0) then + ! Print the filename to stdout + print "(a,4x,a)", "#", trim(fname) + + open (unit=fileID, file=trim(fname)) + + ! Write the title of the file + write (fileID, *) "Title = ""ADflow Surface Solution Data""" + + ! Write the variables header + write (fileID, "(a)", advance="no") ("Variables = ") + + ! Write the variable names + write (fileID, "(a)", advance="no") " ""CoordinateX"" " + write (fileID, "(a)", advance="no") " ""CoordinateY"" " + write (fileID, "(a)", advance="no") " ""CoordinateZ"" " + + ! Write the rest of the variables + do i = 1, nSolVar + write (fileID, "(a)", advance="no") " """//trim(solNames(i))//""" " + end do + write (fileID, "(1x)") + + deallocate (solNames) + + end if + + ! First pass through to generate and write header information + masterBCLoop1: do iBCGroup = 1, nFamExchange + + ! Pointers for easier reading + exch => BCFamExchange(iBCGroup, sps) + zipper => zipperMeshes(iBCGroup) + + ! First thing we do is figure out if we actually need to do + ! anything with this BCgroup at all. If none the requested + ! families are in this BCExcahnge we don't have to do + ! anything. + + BCGroupNeeded = .False. + do i = 1, size(famList) + if (famInLIst(famList(i), exch%famList)) then + BCGroupNeeded = .True. + end if + end do + + ! Keep going if we don't need this. + if (.not. BCGroupNeeded) then + cycle + end if + + ! Get the sizes of this BCGroup + call getSurfaceSize(sizeNode, sizeCell, exch%famList, size(exch%famList), .True.) + call mpi_reduce(sizeNode, nNodes, 1, adflow_integer, MPI_SUM, 0, adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! Now gather up the cell family info. + allocate (cellDisps(0:nProc), cellSizes(nProc)) + + call mpi_gather(sizeCell, 1, adflow_integer, & + cellSizes, 1, adflow_integer, 0, adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + allocate (localElemFam(sizeCell)) + call getSurfaceFamily(localElemFam, sizeCell, exch%famList, size(exch%famList), .True.) + + if (myid == 0) then + cellDisps(0) = 0 + do iProc = 1, nProc + cellDisps(iProc) = cellDisps(iProc - 1) + cellSizes(iProc) + end do + nCells = sum(cellSizes) + allocate (elemFam(nCells)) + end if + + call mpi_gatherv(localElemFam, & + size(localElemFam), adflow_integer, elemFam, & + cellSizes, cellDisps, adflow_integer, 0, adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! Deallocate temp memory + deallocate (localElemFam, cellSizes, cellDisps) + + rootProc: if (myid == 0 .and. nCells > 0) then + do iFam = 1, size(exch%famList) + + ! Check if we have to write this one: + famInclude: if (famInList(exch%famList(iFam), famList)) then + nCellsToWrite = 0 + do i = 1, nCells + ! Check if this elem is to be included + if (elemFam(i) == exch%famList(iFam)) then + nCellsToWrite = nCellsToWrite + 1 + end if + end do + + if (nCellsToWrite > 0) then + write (fileID, "(a,a,a)") "Zone T = """, trim(famNames(exch%famList(iFam))), """" + write (fileID, *) "Nodes = ", nNodes, " Elements = ", & + nCellsToWrite, " ZONETYPE=FEQUADRILATERAL" + write (fileID, *) "DATAPACKING=POINT" + + end if + end if famInclude + end do + end if rootProc + if (myid == 0) then + deallocate (elemFam) + end if + end do masterBCLoop1 + + ! Now do everything again but for real. + masterBCLoop: do iBCGroup = 1, nFamExchange + + ! Pointers for easier reading + exch => BCFamExchange(iBCGroup, sps) + zipper => zipperMeshes(iBCGroup) + + ! First thing we do is figure out if we actually need to do + ! anything with this BCgroup at all. If none the requested + ! families are in this BCExcahnge we don't have to do + ! anything. + + BCGroupNeeded = .False. + do i = 1, size(famList) + if (famInList(famList(i), exch%famList)) then + BCGroupNeeded = .True. + end if + end do + + ! Keep going if we don't need this. + if (.not. BCGroupNeeded) then + cycle + end if + + ! Get the sizes of this BCGroup + call getSurfaceSize(sizeNode, sizeCell, exch%famList, size(exch%famList), .True.) + allocate (nodalValues(max(sizeNode, 1), nSolVar + 3 + 6)) + ! Compute the nodal data + + call computeSurfaceOutputNodalData(BCFamExchange(iBCGroup, sps), & + zipperMeshes(iBCGroup), .False., nodalValues(:, :)) + + ! Gather up the number of nodes to be set to the root proc: + allocate (nodeSizes(nProc), nodeDisps(0:nProc)) + nodeSizes = 0 + nodeDisps = 0 + + call mpi_allgather(sizeNode, 1, adflow_integer, nodeSizes, 1, adflow_integer, & + adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + nodeDisps(0) = 0 + do iProc = 1, nProc + nodeDisps(iProc) = nodeDisps(iProc - 1) + nodeSizes(iProc) + end do + + iSize = 3 + 6 + nSolVar + if (myid == 0) then + nNodes = sum(nodeSizes) + else + nNodes = 1 + end if + + ! Only root proc actually has any space allocated + allocate (vars(nNodes, iSize)) + + ! Gather values to the root proc. + do i = 1, iSize + call mpi_gatherv(nodalValues(:, i), sizeNode, & + adflow_real, vars(:, i), nodeSizes, nodeDisps, & + adflow_real, 0, adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + end do + deallocate (nodalValues) + + ! Now gather up the connectivity + allocate (cellDisps(0:nProc), cellSizes(nProc)) + + call mpi_gather(sizeCell, 1, adflow_integer, & + cellSizes, 1, adflow_integer, 0, adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + if (allocated(cgnsBlockID)) then + deallocate (cgnsBlockID) + end if + + allocate (localConn(4, sizeCell), localElemFam(sizeCell), cgnsBlockID(sizeCell)) + call getSurfaceConnectivity(localConn, cgnsBlockID, sizeCell, exch%famList, size(exch%famList), .True.) + call getSurfaceFamily(localElemFam, sizeCell, exch%famList, size(exch%famList), .True.) + + if (myid == 0) then + cellDisps(0) = 0 + do iProc = 1, nProc + cellDisps(iProc) = cellDisps(iProc - 1) + cellSizes(iProc) + end do + nCells = sum(cellSizes) + allocate (conn(4, nCells)) + allocate (elemFam(nCells)) + end if + + ! We offset the conn array by nodeDisps(iProc) which + ! automagically adjusts the connectivity to account for the + ! number of nodes from different processors + + call mpi_gatherv(localConn + nodeDisps(myid), & + 4 * size(localConn, 2), adflow_integer, conn, & + cellSizes * 4, cellDisps * 4, adflow_integer, 0, adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + call mpi_gatherv(localElemFam, & + size(localElemFam), adflow_integer, elemFam, & + cellSizes, cellDisps, adflow_integer, 0, adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! Local values are finished + deallocate (localConn, localElemFam) + rootProc2: if (myid == 0 .and. nCells > 0) then + + allocate (mask(nCells)) + dataWritten = .False. + do iFam = 1, size(exch%famList) + + ! Check if we have to write this one: + famInclude2: if (famInList(exch%famList(iFam), famList)) then + + ! Create a temporary mask + mask = 0 + nCellsToWrite = 0 + do i = 1, nCells + ! Check if this elem is to be included + if (elemFam(i) == exch%famList(iFam)) then + mask(i) = 1 + nCellsToWrite = nCellsToWrite + 1 + end if + end do + + actualWrite2: if (nCellsToWrite > 0) then + + ! Dump the coordinates + do k = 1, nNodes + do j = 1, 3 + write (fileID, sci6, advance='no') vars(k, j) + end do + do j = 1, nSolVar + write (fileID, sci6, advance='no') vars(k, j + 9) + end do + write (fileID, "(1x)") + end do + + ! Dump the connectivity + do i = 1, nCells + ! Check if this elem is to be included + if (mask(i) == 1) then + do k = 1, 4 + write (fileID, int7, advance="no") conn(k, i) + end do + write (fileID, "(1x)") + end if + end do + + end if actualWrite2 + end if famInclude2 + end do + deallocate (mask) + end if rootProc2 + if (myid == 0) then + deallocate (conn, elemFam) + end if + deallocate (cellSizes, cellDisps, nodeSizes, nodeDisps, vars) + end do masterBCLoop + + if (myid == 0) then + close (fileID) + end if + + end do spectralLoop + + ! Restore the modified option + surfWriteBlank = blankSave + if (myID == 0 .and. printIterations) then + print "(a)", "# Tecplot surface file(s) written" + print "(a)", "#" + end if + end subroutine writeBCSurfacesASCII + end module tecplotIO diff --git a/src/solver/userSurfaceIntegrations.F90 b/src/solver/userSurfaceIntegrations.F90 index 760cca989..8124834f8 100644 --- a/src/solver/userSurfaceIntegrations.F90 +++ b/src/solver/userSurfaceIntegrations.F90 @@ -5,6 +5,227 @@ module userSurfaceIntegrations use userSurfaceIntegrationData contains + subroutine remapMesh(vars, conn, ptValid, nvars, npts, nElements, varsValid, connValid, nptsValid) + implicit none + + ! Inputs + integer(kind=intType), intent(in) :: npts, nElements, nvars ! Number of nodes and elements + real(kind=realType), dimension(nvars, npts), intent(in) :: vars ! Original nodes array (shape: 3 x npts) + integer(kind=intType), dimension(3, nElements), intent(in) :: conn ! Original connectivity array (shape: 3 x n_elements) + logical, dimension(npts), intent(in) :: ptValid ! Logical array indicating validity of nodes + + ! Outputs + real(kind=realType), dimension(:, :), allocatable, intent(out) :: varsValid ! Valid nodes array (shape: 3 x npts_valid) + integer(kind=intType), dimension(:, :), allocatable, intent(out) :: connValid ! Valid connectivity array (shape: 3 x n_elements_valid) + integer(kind=intType), intent(out) :: nptsValid ! Number of valid nodes + + ! Local variables + logical, dimension(:), allocatable :: nodeUsed ! Array to track if a node is used in a valid triangle + integer(kind=intType) :: i, j, k, countPtsValid, countConnsValid + integer(kind=intType), dimension(3, nElements) :: connTemp + + ! Allocate memory for node_used array + allocate (nodeUsed(npts)) + + ! Initialize node_used array + nodeUsed = .false. + + ! Count valid nodes + countPtsValid = 0 + do i = 1, npts + if (ptValid(i)) then + countPtsValid = countPtsValid + 1 + nodeUsed(i) = .true. + end if + end do + + ! Allocate memory for valid nodes array + allocate (varsValid(size(vars, 1), countPtsValid)) + + ! Count the number of valid elements + countConnsValid = 0 + do i = 1, nElements + if ((ptValid(conn(1, i)) .and. ptValid(conn(2, i)) .and. ptValid(conn(3, i)))) then + countConnsValid = countConnsValid + 1 + end if + end do + + ! Allocate memory for valid connectivity array + ! Copy valid nodes to new array + allocate (connValid(3, countConnsValid)) + + countPtsValid = 0 + do i = 1, npts + if (nodeUsed(i)) then + countPtsValid = countPtsValid + 1 + varsValid(:, countPtsValid) = vars(:, i) + + ! Update the connectivity array + do j = 1, nElements + if (conn(1, j) == i) then + connTemp(1, j) = countPtsValid + end if + if (conn(2, j) == i) then + connTemp(2, j) = countPtsValid + end if + if (conn(3, j) == i) then + connTemp(3, j) = countPtsValid + end if + end do + end if + end do + + ! Copy the valid connectivity array + countConnsValid = 0 + do i = 1, nElements + if ((ptValid(conn(1, i)) .and. ptValid(conn(2, i)) .and. ptValid(conn(3, i)))) then + countConnsValid = countConnsValid + 1 + connValid(:, countConnsValid) = connTemp(:, i) + end if + end do + + ! Deallocate node_used array + deallocate (nodeUsed) + + end subroutine remapMesh + + subroutine getUserSurfaceData(iSurf, sps, validVars, validConn) + + use constants + use block, only: flowDoms, nDom + use flowVarRefState, only: pRef, rhoRef + use communication, only: myid + use flowUtils, only: computePtot, computeTtot + use sorting, only: famInList + use utils, only: terminate + + ! Input Parameters + integer(kind=intType), intent(in) :: iSurf, sps + + ! Output Parameters + real(kind=realType), dimension(:, :), allocatable, intent(out) :: validVars + integer(kind=intType), dimension(:, :), allocatable, intent(out) :: validConn + + ! Working parameters + integer(kind=intType) :: i, j, k, jj, ierr, nn, iDim, nPts + real(kind=realType), dimension(:), allocatable :: recvBuffer1, recvBuffer2 + real(kind=realType), dimension(:, :), allocatable :: vars + real(kind=realType), dimension(:, :), allocatable :: solVars + logical, dimension(:), allocatable :: ptValid + type(userIntSurf), pointer :: surf + real(kind=realType) :: p, gamma, vmag + + if (nUserIntSurfs == 0) then + write (*, *) "Warning: No user-supplied surfaces to write" + return ! Nothing to do + end if + + ! Set the pointers for the required communication variables + domainLoop: do nn = 1, nDom + if (flowDoms(nn, 1, sps)%addGridVelocities) then + call terminate("userSurfaceIntegrations", "Cannot use user-supplied surface integrations"& + &"on with moving grids") + end if + + flowDoms(nn, 1, sps)%realCommVars(iRho)%var => flowDoms(nn, 1, sps)%w(:, :, :, iRho) + flowDoms(nn, 1, sps)%realCommVars(iVx)%var => flowDoms(nn, 1, sps)%w(:, :, :, iVx) + flowDoms(nn, 1, sps)%realCommVars(iVy)%var => flowDoms(nn, 1, sps)%w(:, :, :, iVy) + flowDoms(nn, 1, sps)%realCommVars(iVz)%var => flowDoms(nn, 1, sps)%w(:, :, :, iVz) + flowDoms(nn, 1, sps)%realCommVars(iZippFlowP)%var => flowDoms(nn, 1, sps)%P(:, :, :) + flowDoms(nn, 1, sps)%realCommVars(iZippFlowGamma)%var => flowDoms(nn, 1, sps)%gamma(:, :, :) + ! flowDoms(nn, 1, sps)%realCommVars(iZippFlowSface)%var => Not Implemented + + flowDoms(nn, 1, sps)%realCommVars(iZippFlowX)%var => flowDoms(nn, 1, sps)%x(:, :, :, 1) + flowDoms(nn, 1, sps)%realCommVars(iZippFlowY)%var => flowDoms(nn, 1, sps)%x(:, :, :, 2) + flowDoms(nn, 1, sps)%realCommVars(iZippFlowZ)%var => flowDoms(nn, 1, sps)%x(:, :, :, 3) + end do domainLoop + + ! Pointer for easier reading + surf => userIntSurfs(iSurf) + + nPts = size(surf%pts, 2) + + ! Communicate the face values and the nodal values + if (myid == 0) then + allocate (recvBuffer1(6 * nPts), recvBuffer2(3 * nPts)) + else + allocate (recvBuffer1(0), recvBuffer2(0)) + allocate (solVars(0, 0)) + end if + + call commUserIntegrationSurfaceVars(recvBuffer1, iRho, iZippFlowGamma, surf%flowComm) + call commUserIntegrationSurfaceVars(recvBuffer2, iZippFlowX, iZippFlowZ, surf%nodeComm) + + ! Only get the data on the root processor + if (myid == 0) then + + ! Allocate some temporary data + allocate (vars(npts, nZippFlowComm), ptValid(npts)) + ptValid = .True. + + ! We have to re-order the data according to the "inv" + ! array in each of the two comms. + do i = 1, nPts + + ! Flow Variables + j = surf%flowComm%inv(i) + vars(j, iRho:iZippFlowGamma) = recvBuffer1(6 * (i - 1) + iRho:6 * (i - 1) + iZippFlowGamma) + + if (.not. surf%flowComm%valid(i)) then + ptValid(j) = .False. + end if + + ! Sface is not implemented. To correctly do this, + ! interpolate the three components of 's', do the dot + ! product with the local normal to get the sFace value. + vars(j, iZippFlowSface) = zero + + ! Node Comm Values + j = surf%nodeComm%inv(i) + vars(j, iZippFlowX:iZippFlowZ) = recvBuffer2(3 * i - 2:3 * i) + + ! The additional pt-valid array + if (.not. surf%nodeComm%valid(i)) then + ptValid(j) = .False. + end if + end do + + ! We need to loop over the variables to compute and store the solution variables + allocate (solVars(11, npts)) + + do i = 1, npts + p = vars(i, iRhoE) + gamma = vars(i, iZippFlowGamma) + solVars(1, i) = vars(i, iZippFlowX) ! x-coordinate + solVars(2, i) = vars(i, iZippFlowY) ! y-coordinate + solVars(3, i) = vars(i, iZippFlowZ) ! z-coordinate + solVars(4, i) = vars(i, iVx) ! x-velocity + solVars(5, i) = vars(i, iVy) ! y-velocity + solVars(6, i) = vars(i, iVz) ! z-velocity + solVars(7, i) = vars(i, iRho) ! density + solVars(8, i) = vars(i, iRhoE) * pRef ! Static pressure + + ! Total pressure, stored in the 9th column + call computePtot(vars(i, iRho), vars(i, iVx), vars(i, iVy), vars(i, iVz), p, solVars(9, i)) + + ! Total temperature, stored in the 10th column + call computeTtot(vars(i, iRho), vars(i, iVx), vars(i, iVy), vars(i, iVz), p, solVars(10, i)) + + ! Compute the Mach number + vmag = sqrt(vars(i, iVx)**2 + vars(i, iVy)**2 + vars(i, iVz)**2) + solVars(11, i) = vmag / sqrt(gamma * p / vars(i, iRho)) ! Mach number + end do + + ! Remap the mesh to only include valid points + call remapMesh(solVars, surf%conn, ptValid, size(solVars, 1), npts, & + size(surf%conn, 2), validVars, validConn, nptsValid) + + deallocate (vars, ptValid) + end if + deallocate (recvBuffer1, recvBuffer2) + + end subroutine getUserSurfaceData + subroutine integrateUserSurfaces(localValues, famList, sps) use constants