Loading...
Searching...
No Matches
Public Types |
Public Member Functions |
Public Attributes |
Protected Member Functions |
Protected Attributes |
Static Protected Attributes |
Friends |
List of all members
FvStructuredSolver3D Class Reference
3D structured solver class More...
#include <fvstructuredsolver3d.h>
Inheritance diagram for FvStructuredSolver3D:
Collaboration diagram for FvStructuredSolver3D:
Public Types | |
| typedef void(FvStructuredSolver3D::* | fluxmethod) (MFloat *, MFloat *, MInt, MInt) |
 Public Types inherited from FvStructuredSolver< 3 > | |
| using | Timers = maia::structured::Timers_ |
Public Member Functions | |
| FvStructuredSolver3D (MInt, StructuredGrid< 3 > *, MBool *, const MPI_Comm comm) | |
| ~FvStructuredSolver3D () | |
| void | addGhostPointCoordinateValues () |
| Extrapolates and exchanges ghost point coordinates. More... | |
| void | extrapolateGhostPointCoordinatesBC () |
| void | initSolutionStep (MInt mode) |
| void | initialCondition () override |
| Computation of infinity values for the conservative and primitive variables. More... | |
| void | assignBndryCells () |
| void | nonReflectingBC () |
| void | applyBoundaryCondition () override |
| void | calcSurfaceMetrics () |
| void | initBndryCnds () |
| void | allocateSingularities () |
| void | applySandpaperTrip () |
| void | applySandpaperTripAirfoil () |
| void | initSandpaperTrip () |
| void | tripForceCoefficients (MFloat *, MFloat *, MFloat *, MInt, MInt) |
| void | tripFourierCoefficients (MFloat *, MInt, MFloat, MFloat) |
| void | computeDomainWidth () override |
| void | computeZonalConnections () |
| void | zonalExchange () override |
| void | zonalAllreduce () |
| void | zonalRealInterpolation () |
| void | zonalGather () |
| void | zonalSend () |
| void | zonalReceive () |
| void | zonalCheck () |
| void | zonalBufferAverage () |
| void | zonalScatter () |
| void | spanwiseAvgZonal (std::vector< MFloat * > &) override |
| void | distributeFluxToCells () |
| void | computeTimeStep () override |
| void | computeVorticity () override |
| void | computeLambda2Criterion () override |
| Function to compute the lambda_2 criterion. More... | |
| MFloat | computeTotalKineticEngergy () |
| MFloat | computeTotalPressure () |
| void | initializeNeighbourCellIds () |
| void | initFluxMethod () |
| void | computePrimitiveVariables () override |
| template<MFloat(FvStructuredSolver::*)(MInt) const = &FvStructuredSolver::dummy> | |
| void | computePrimitiveVariables_ () |
| MFloat | pressure (MInt) |
| MFloat | getPSI (MInt, MInt) |
| void | Ausm () |
| void | AusmLES (MFloat *QLeft, MFloat *QRight, const MInt dim, const MInt cellId) |
| AUSM central. More... | |
| void | AusmLES_PTHRC (MFloat *QLeft, MFloat *QRight, MInt dim, MInt cellId) |
| AUSM PTHRC. More... | |
| void | AusmDV (MFloat *QLeft, MFloat *QRight, MInt dim, MInt cellId) |
| template<MInt noVars> | |
| void | Muscl_AusmLES () |
| template<MInt noVars> | |
| void | Muscl_AusmLES_PTHRC () |
| template<MInt noVars> | |
| void | Muscl_AusmDV () |
| template<fluxmethod ausm, MInt noVars> | |
| void | MusclStretched_ () |
| void | Muscl (MInt timerId=-1) override |
| void | MusclRANS () |
| void | MusclVenkatakrishnan3D () |
| MUSCL with Venkatakrishan limiter. More... | |
| void | MusclAlbada () |
| void | MusclMinModLimiter () |
| void | computeCellLength () |
| void | computeReconstructionConstantsSVD () |
| void | VENKATAKRISHNAN_MOD_FCT (MFloat effNghbrDelta, MFloat srfcDelta, MFloat dxEpsSqr, MInt cellPos, MInt var, MFloatScratchSpace &minPhi) |
| Venkatakrishnan limiter, modified for better results. More... | |
| void | VENKATAKRISHNAN_FCT (MFloat effNghbrDelta, MFloat srfcDelta, MFloat dxEpsSqr, MInt cellPos, MInt var, MFloatScratchSpace &minPhi) |
| Standard Venkatakrishnan limiter. More... | |
| void | BARTH_JESPERSON_FCT (MFloat effNghbrDelta, MFloat srfcDelta, MFloat dxEpsSqr, MInt cellPos, MInt var, MFloatScratchSpace &minPhi) |
| Barth-Jesperson Limiter. More... | |
| MBool | maxResidual () |
| Computation of the maximum residual. More... | |
| MBool | rungeKuttaStep () |
| void | updateSpongeLayer () |
| void | addDisturbance () |
| void | viscousFlux () |
| void | viscousFluxRANS () |
| template<MBool twoEqRans = false> | |
| void | viscousFluxLES () |
| Viscous flux computation. More... | |
| void | moveGrid (const MBool isRestart, const MBool zeroPos) override |
| void | applyBodyForce (const MBool isRestart, const MBool zeroPos) override |
| void | initMovingGrid () override |
| void | initBodyForce () override |
| void | applyViscousBoundaryCondition () |
| void | applyInviscidBoundaryCondition () |
| virtual void | computeFrictionPressureCoef (MBool computePower) override |
| void | initFFTW (fftw_complex *uPhysField, fftw_complex *vPhysField, fftw_complex *wPhysField, MInt lx, MInt ly, MInt lz, MInt noPeakModes) |
| void | getFourierCoefficients (MFloat *k, MFloat k0, std::complex< MFloat > *fourierCoefficient) |
| Generates a single complex coefficient of Fourier series. More... | |
| MFloat | randnormal (MFloat mu, MFloat sigma) |
| Returns a normal distributed random-number with mu=mean and sigma=standard deviation. More... | |
| void | loadRestartBC2600 () |
| void | loadRestartBC2601 () |
| void | loadRestartSTG (MBool) |
| void | gather (const MBool, std::vector< std::unique_ptr< StructuredComm< 3 > > > &) override |
| void | scatter (const MBool, std::vector< std::unique_ptr< StructuredComm< 3 > > > &) override |
| void | waveGather () |
| void | waveScatter () |
| void | waveSend (std::vector< MPI_Request > &) |
| void | waveReceive (std::vector< MPI_Request > &) |
| void | waveExchange () override |
| void | spanwiseWaveReorder () override |
| void | gcFillGhostCells (std::vector< MFloat * > &) |
| void | gcExtrapolate (std::vector< MFloat * > &) |
| void | viscousFluxCorrection () |
| void | computeCumulativeAverage (MBool forceReset) override |
| void | initInterpolatedPoints () |
| Compute all the interpolation coefficients necessary for the interpolation output. More... | |
| void | saveInterpolatedPoints () override |
| void | saveNodalBoxes () override |
| void | savePointsToAsciiFile (MBool) override |
| Saves variables of given cells to ASCII file. More... | |
| void | initPointsToAsciiFile () override |
| virtual void | computePorousRHS (MBool) override |
| void | exchange6002 () |
| Public Member Functions inherited from FvStructuredSolver< 3 > | |
| FvStructuredSolver (MInt solverId, StructuredGrid< nDim > *, MBool *propertiesGroups, const MPI_Comm comm) | |
| Constructor of the structured solver. More... | |
| MBool | isActive () const |
| void | initializeFQField () |
| Counts the number of necessary FQ fields, allocates them and corrects the indexes of the FQ variable pointers. More... | |
| MInt | timer (const MInt timerId) const |
| MInt | noSolverTimers (const MBool allTimings) override |
| virtual void | writeHeaderAttributes (ParallelIoHdf5 *pio, MString fileType) |
| Overloaded version of writeHeaderAttributes that receives ParallelIoHdf5 object pointer instead of 'fileId'. More... | |
| virtual void | writePropertiesAsAttributes (ParallelIoHdf5 *pio, MString path) |
| Overloaded version of writePropertiesAsAttributes that receives ParallelIoHdf5 object pointer instead of 'fileId'. More... | |
| void | saveSolverSolution (MBool=false, const MBool=false) override |
| void | saveSolution (MBool) |
| Saves the soution to hdf5 file. More... | |
| void | savePartitions () |
| Saves the partitioned grid into an HDF5 file. Not used in production use but useful for debugging. More... | |
| void | saveBoxes () |
| virtual void | savePointsToAsciiFile (MBool) |
| virtual void | initPointsToAsciiFile () |
| virtual void | saveInterpolatedPoints () |
| virtual void | saveNodalBoxes () |
| virtual void | reIntAfterRestart (MBool) |
| MBool | prepareRestart (MBool, MBool &) override |
| Prepare the solvers for a grid-restart. More... | |
| void | writeRestartFile (MBool) override |
| void | writeRestartFile (const MBool, const MBool, const MString, MInt *) override |
| void | shiftCellValuesRestart (MBool) |
| void | loadRestartFile () |
| Load Restart File (primitive and conservative output) general formulation. More... | |
| MBool | loadBoxFile (MString, MString, MInt, MInt) |
| Load Box file general formulation. More... | |
| virtual void | loadRestartBC2600 () |
| virtual void | loadRestartBC2601 () |
| virtual void | loadRestartSTG (MBool) |
| void | saveForcesToAsciiFile (MBool) |
| Function to save the force coefficients and power to an ASCII file. More... | |
| void | saveAveragedVariables (MString, MInt, MFloat **) |
| Saves the averaged (mean) variables from postprocessing to an HDF5 file. More... | |
| void | saveProductionTerms (MString, MFloat **) |
| Writes the production terms into a given file. More... | |
| void | saveDissipation (MString, MFloat *) |
| Writes the dissipation into a given file. More... | |
| void | saveGradients (MString, MFloat **, MString *) |
| Writes the gradients into a given file. More... | |
| void | saveAverageRestart (MString, MInt, MFloat **, MFloat **, MFloat **, MFloat **) |
| Writes an restart file for postprocessing. More... | |
| void | getSolverTimings (std::vector< std::pair< MString, MFloat > > &solverTimings, const MBool allTimings) override |
| Get solver timings. More... | |
| void | getDomainDecompositionInformation (std::vector< std::pair< MString, MInt > > &domainInfo) override |
| Return decomposition information, i.e. number of local elements,... More... | |
| virtual MFloat | getCellLengthY (MInt, MInt, MInt) |
| virtual MFloat | getCellCoordinate (MInt, MInt) |
| MInt | noInternalCells () const override |
| Return the number of internal cells within this solver. More... | |
| void | saveAuxData () |
| void | saveForceCoefficient (ParallelIoHdf5 *parallelIoHdf5) |
| Saves force coefficients to an HDF5 file. More... | |
| void | computeAuxData () |
| void | computeAuxDataRoot () |
| virtual void | computeDomainWidth () |
| void | computeForceCoef () |
| Function to compute the force coefficient cl, split split into the viscous part cLv and the pressure part cLp. More... | |
| void | computeForceCoefRoot () |
| Function to compute the coefficient, split split into the viscous part cLv and the pressure part cLp The ROOT version is faster due to an MPI_Reduce instead of an MPI_Allreduce, but only root rank has data. More... | |
| virtual void | computeFrictionPressureCoef (MBool)=0 |
| virtual | ~FvStructuredSolver () |
| MFloat | computeRecConstSVD (const MInt ijk, const MInt noNghbrIds, MInt *nghbr, MInt ID, MInt sID, MFloatScratchSpace &tmpA, MFloatScratchSpace &tmpC, MFloatScratchSpace &weights, const MInt recDim) |
| AUX DATA ENDS /////////////////////////////////////////////////////////////. More... | |
| void | initializeFvStructuredSolver (MBool *propertiesGroups) |
| Structured Solver Constructor reads and allocate properties/variables. More... | |
| void | allocateAndInitBlockMemory () |
| void | allocateAuxDataMaps () |
| AUX DATA //////////////////////////////////////////////////////////////////. More... | |
| void | setRungeKuttaProperties () |
| This function reads the properties required for Runge Kutta time stepping. More... | |
| void | setSamplingProperties () |
| void | setNumericalProperties () |
| Reads and initializes properties associated with the numerical method. More... | |
| void | setInputOutputProperties () |
| Reads properties and initializes variables associated with input/output. More... | |
| void | setZonalProperties () |
| Set which zones are RANS and which are LES or if full LES or full RANS. More... | |
| void | allocateVariables () |
| void | setTestcaseProperties () |
| Reads and initializes properties associated with the Testcase. More... | |
| void | setMovingGridProperties () |
| Reads and initializes properties associated with the Moving Grid Methods. More... | |
| void | setBodyForceProperties () |
| Reads and initializes properties associated with the Moving Grid Methods. More... | |
| void | setPorousProperties () |
| Set properties for porous blocks. More... | |
| void | readAndSetSpongeLayerProperties () |
| void | setSTGProperties () |
| void | setProfileBCProperties () |
| void | createMPIGroups () |
| void | readAndSetAuxDataMap () |
| void | computePV () |
| MFloat | dummy (MInt) const |
| MFloat | pressure_twoEqRans (MInt cellId) const |
| MFloat | totalEnergy_twoEqRans (MInt cellId) const |
| void | partitionGrid () |
| virtual void | computePrimitiveVariables () |
| virtual void | computeConservativeVariables () |
| void | computeConservativeVariables_ () |
| void | saveVarToPrimitive (MInt, MInt, MFloat) |
| virtual void | gcFillGhostCells (std::vector< MFloat * > &) |
| void | computeSamplingInterval () |
| void | checkNans () |
| Checks whole domain for NaNs and adds the number of NaNs globally. More... | |
| void | setVolumeForce () |
| void | computeVolumeForces () |
| virtual void | computePorousRHS (MBool) |
| void | initPorous () |
| virtual void | applyBoundaryCondtition () |
| virtual void | moveGrid (MInt) |
| virtual void | moveGrid (MBool, MBool) |
| virtual void | initMovingGrid () |
| virtual void | initBodyForce () |
| virtual void | applyBodyForce (MBool, MBool) |
| virtual void | computeLambda2Criterion () |
| virtual void | computeVorticity () |
| void | exchange () |
| SVD STUFF ENDS /////////////////////////////////////////////////////////////. More... | |
| void | exchange (std::vector< std::unique_ptr< StructuredComm< nDim > > > &, std::vector< std::unique_ptr< StructuredComm< nDim > > > &) |
| Parallel exchange of primitive variables between partitions with MPI. More... | |
| void | send (const MBool, std::vector< std::unique_ptr< StructuredComm< nDim > > > &, std::vector< MPI_Request > &) |
| void | receive (const MBool, std::vector< std::unique_ptr< StructuredComm< nDim > > > &, std::vector< MPI_Request > &) |
| virtual void | gather (const MBool, std::vector< std::unique_ptr< StructuredComm< nDim > > > &) |
| virtual void | scatter (const MBool, std::vector< std::unique_ptr< StructuredComm< nDim > > > &) |
| MBool | isPeriodicComm (std::unique_ptr< StructuredComm< nDim > > &) |
| MBool | skipPeriodicDirection (std::unique_ptr< StructuredComm< nDim > > &) |
| virtual void | zonalExchange () |
| virtual void | spanwiseAvgZonal (std::vector< MFloat * > &) |
| virtual void | waveExchange () |
| virtual void | spanwiseWaveReorder () |
| void | setTimeStep () |
| void | setLimiterVisc () |
| void | fixTimeStepTravelingWave () |
| void | exchangeTimeStep () |
| void | initializeRungeKutta () |
| virtual void | computeTimeStep () |
| MBool | isInInterval (MInt) |
| MInt | getNoCells () |
| MInt | noVariables () const override |
| Return the number of variables. More... | |
| MInt | getNoActiveCells () |
| MInt * | getActiveCells () |
| MInt * | getOffsetCells () |
| MInt | getNoGhostLayers () |
| MInt | getWaveAvrgInterval () |
| MInt | getWaveStepOffset () |
| MInt * | getCellGrid () |
| MBool | isMovingGrid () |
| MInt | getGridMovingMethod () |
| MInt | getBodyForceMethod () |
| MBool | isStreamwiseTravelingWave () |
| MBool | isTravelingWave () |
| StructuredGrid< nDim > * | getGrid () |
| MFloat | getGamma () |
| MFloat | getSutherlandConstant () |
| MFloat | getRe0 () |
| MFloat | getMa () |
| MPI_Comm | getCommunicator () |
| virtual void | computeCumulativeAverage (MBool) |
| virtual void | loadSampleFile (MString) |
| virtual void | getSampleVariables (MInt, MFloat *) |
| MFloat | getPV (MInt var, MInt cellId) |
| void | setPV (MInt var, MInt cellId, MFloat value) |
| virtual MFloat | getSampleVorticity (MInt, MInt) |
| virtual MFloat | dvardxyz (MInt, MInt, MFloat *) |
| virtual MFloat | dvardx (MInt, MFloat *) |
| virtual void | loadAverageRestartFile (const char *, MFloat **, MFloat **, MFloat **, MFloat **) |
| virtual void | loadAveragedVariables (const char *) |
| void | convertRestartVariables (MFloat oldMa) |
| virtual void | convertRestartVariablesSTG (MFloat oldMa) |
| virtual bool | rungeKuttaStep ()=0 |
| virtual void | viscousFlux ()=0 |
| virtual void | Muscl (MInt=-1)=0 |
| virtual void | applyBoundaryCondition ()=0 |
| virtual void | initSolutionStep (MInt)=0 |
| virtual void | initialCondition ()=0 |
| void | tred2 (MFloatScratchSpace &A, MInt dim, MFloat *diag, MFloat *offdiag) |
| Householder Reduction according to Numercial Recipies in C: The Art of Scientific Computing. More... | |
| void | tqli2 (MFloat *diag, MFloat *offdiag, MInt dim) |
| Compute Eigenvalues with implicit shift according to Numercial Recipies in C: The Art of Scientific Computing. More... | |
| void | insertSort (MInt dim, MFloat *list) |
| Sorting function to sort list in ascending order. More... | |
| MFloat | pythag (MFloat a, MFloat b) |
| void | resetRHS () |
| Reset the right hand side to zero. More... | |
| void | rhs () |
| void | rhsBnd () |
| void | lhsBnd () |
| void | initSolver () override |
| virtual MBool | maxResidual () |
| MBool | solutionStep () override |
| void | preTimeStep () override |
| void | postTimeStep () override |
| : Performs the post time step More... | |
| void | finalizeInitSolver () override |
| void | cleanUp () override |
| virtual void | updateSpongeLayer ()=0 |
| MFloat | time () const override |
| Return the time. More... | |
| MInt | determineRestartTimeStep () const override |
| Load the restart time step from the restart file (useNonSpecifiedRestartFile enabled) More... | |
| MBool | hasRestartTimeStep () const override |
| void | init () |
| std::array< MInt, nDim > | beginP0 () |
| std::array< MInt, nDim > | beginP1 () |
| std::array< MInt, nDim > | beginP2 () |
| std::array< MInt, nDim > | endM2 () |
| std::array< MInt, nDim > | endM1 () |
| std::array< MInt, nDim > | endM0 () |
| Public Member Functions inherited from Solver | |
| MString | getIdentifier (const MBool useSolverId=false, const MString preString="", const MString postString="_") |
| virtual | ~Solver ()=default |
| virtual MInt | noInternalCells () const =0 |
| Return the number of internal cells within this solver. More... | |
| virtual MFloat | time () const =0 |
| Return the time. More... | |
| virtual MInt | noVariables () const |
| Return the number of variables. More... | |
| virtual void | getDimensionalizationParams (std::vector< std::pair< MFloat, MString > > &) const |
| Return the dimensionalization parameters of this solver. More... | |
| void | updateDomainInfo (const MInt domainId, const MInt noDomains, const MPI_Comm mpiComm, const MString &loc) |
| Set new domain information. More... | |
| virtual MFloat & | a_slope (const MInt, MInt const, const MInt) |
| virtual MBool | a_isBndryCell (const MInt) const |
| virtual MFloat & | a_FcellVolume (MInt) |
| virtual MInt | getCurrentTimeStep () const |
| virtual void | accessSampleVariables (MInt, MFloat *&) |
| virtual void | getSampleVariableNames (std::vector< MString > &NotUsed(varNames)) |
| virtual MBool | a_isBndryGhostCell (MInt) const |
| virtual void | saveCoarseSolution () |
| virtual void | getSolverSamplingProperties (std::vector< MInt > &NotUsed(samplingVarIds), std::vector< MInt > &NotUsed(noSamplingVars), std::vector< std::vector< MString > > &NotUsed(samplingVarNames), const MString NotUsed(featureName)="") |
| virtual void | initSolverSamplingVariables (const std::vector< MInt > &NotUsed(varIds), const std::vector< MInt > &NotUsed(noSamplingVars)) |
| virtual void | calcSamplingVariables (const std::vector< MInt > &NotUsed(varIds), const MBool NotUsed(exchange)) |
| virtual void | calcSamplingVarAtPoint (const MFloat *NotUsed(point), const MInt NotUsed(id), const MInt NotUsed(sampleVarId), MFloat *NotUsed(state), const MBool NotUsed(interpolate)=false) |
| virtual void | balance (const MInt *const NotUsed(noCellsToReceiveByDomain), const MInt *const NotUsed(noCellsToSendByDomain), const MInt *const NotUsed(targetDomainsByCell), const MInt NotUsed(oldNoCells)) |
| Perform load balancing. More... | |
| virtual MBool | hasSplitBalancing () const |
| Return if load balancing for solver is split into multiple methods or implemented in balance() More... | |
| virtual void | balancePre () |
| virtual void | balancePost () |
| virtual void | finalizeBalance () |
| virtual void | resetSolver () |
| Reset the solver/solver for load balancing. More... | |
| virtual void | cancelMpiRequests () |
| Cancel open mpi (receive) requests in the solver (e.g. due to interleaved execution) More... | |
| virtual void | setCellWeights (MFloat *) |
| Set cell weights. More... | |
| virtual MInt | noLoadTypes () const |
| virtual void | getDefaultWeights (MFloat *NotUsed(weights), std::vector< MString > &NotUsed(names)) const |
| virtual void | getLoadQuantities (MInt *const NotUsed(loadQuantities)) const |
| virtual MFloat | getCellLoad (const MInt NotUsed(cellId), const MFloat *const NotUsed(weights)) const |
| virtual void | limitWeights (MFloat *NotUsed(weights)) |
| virtual void | localToGlobalIds () |
| virtual void | globalToLocalIds () |
| virtual MInt | noCellDataDlb () const |
| Methods to inquire solver data information. More... | |
| virtual MInt | cellDataTypeDlb (const MInt NotUsed(dataId)) const |
| virtual MInt | cellDataSizeDlb (const MInt NotUsed(dataId), const MInt NotUsed(cellId)) |
| virtual void | getCellDataDlb (const MInt NotUsed(dataId), const MInt NotUsed(oldNoCells), const MInt *const NotUsed(bufferIdToCellId), MInt *const NotUsed(data)) |
| virtual void | getCellDataDlb (const MInt NotUsed(dataId), const MInt NotUsed(oldNoCells), const MInt *const NotUsed(bufferIdToCellId), MLong *const NotUsed(data)) |
| virtual void | getCellDataDlb (const MInt NotUsed(dataId), const MInt NotUsed(oldNoCells), const MInt *const NotUsed(bufferIdToCellId), MFloat *const NotUsed(data)) |
| virtual void | setCellDataDlb (const MInt NotUsed(dataId), const MInt *const NotUsed(data)) |
| virtual void | setCellDataDlb (const MInt NotUsed(dataId), const MLong *const NotUsed(data)) |
| virtual void | setCellDataDlb (const MInt NotUsed(dataId), const MFloat *const NotUsed(data)) |
| virtual void | getGlobalSolverVars (std::vector< MFloat > &NotUsed(globalFloatVars), std::vector< MInt > &NotUsed(globalIntVars)) |
| virtual void | setGlobalSolverVars (std::vector< MFloat > &NotUsed(globalFloatVars), std::vector< MInt > &NotUsed(globalIdVars)) |
| void | enableDlbTimers () |
| void | reEnableDlbTimers () |
| void | disableDlbTimers () |
| MBool | dlbTimersEnabled () |
| void | startLoadTimer (const MString name) |
| void | stopLoadTimer (const MString &name) |
| void | stopIdleTimer (const MString &name) |
| void | startIdleTimer (const MString &name) |
| MBool | isLoadTimerRunning () |
| virtual MInt | noSolverTimers (const MBool NotUsed(allTimings)) |
| virtual void | getSolverTimings (std::vector< std::pair< MString, MFloat > > &NotUsed(solverTimings), const MBool NotUsed(allTimings)) |
| virtual void | getDomainDecompositionInformation (std::vector< std::pair< MString, MInt > > &NotUsed(domainInfo)) |
| void | setDlbTimer (const MInt timerId) |
| virtual void | prepareAdaptation (std::vector< std::vector< MFloat > > &, std::vector< MFloat > &, std::vector< std::bitset< 64 > > &, std::vector< MInt > &) |
| virtual void | reinitAfterAdaptation () |
| virtual void | prepareAdaptation () |
| prepare adaptation for split adaptation before the adaptation loop More... | |
| virtual void | setSensors (std::vector< std::vector< MFloat > > &, std::vector< MFloat > &, std::vector< std::bitset< 64 > > &, std::vector< MInt > &) |
| set solver sensors for split adaptation within the adaptation loop More... | |
| virtual void | saveSensorData (const std::vector< std::vector< MFloat > > &, const MInt &, const MString &, const MInt *const) |
| virtual void | postAdaptation () |
| post adaptation for split adaptation within the adaptation loop More... | |
| virtual void | finalizeAdaptation () |
| finalize adaptation for split sadptation after the adaptation loop More... | |
| virtual void | refineCell (const MInt) |
| Refine the given cell. More... | |
| virtual void | removeChilds (const MInt) |
| Coarsen the given cell. More... | |
| virtual void | removeCell (const MInt) |
| Remove the given cell. More... | |
| virtual void | swapCells (const MInt, const MInt) |
| Swap the given cells. More... | |
| virtual void | swapProxy (const MInt, const MInt) |
| Swap the given cells. More... | |
| virtual MInt | cellOutside (const MFloat *, const MInt, const MInt) |
| Check whether cell is outside the fluid domain. More... | |
| virtual void | resizeGridMap () |
| Swap the given cells. More... | |
| virtual MBool | prepareRestart (MBool, MBool &) |
| Prepare the solvers for a grid-restart. More... | |
| virtual void | reIntAfterRestart (MBool) |
| MPI_Comm | mpiComm () const |
| Return the MPI communicator used by this solver. More... | |
| virtual MInt | domainId () const |
| Return the domainId (rank) More... | |
| virtual MInt | noDomains () const |
| virtual MBool | isActive () const |
| void | setSolverStatus (const MBool status) |
| MBool | getSolverStatus () |
| Get the solver status indicating if the solver is currently active in the execution recipe. More... | |
| MString | testcaseDir () const |
| Return the testcase directory. More... | |
| MString | outputDir () const |
| Return the directory for output files. More... | |
| MString | restartDir () const |
| Return the directory for restart files. More... | |
| MString | solverMethod () const |
| Return the solverMethod of this solver. More... | |
| MString | solverType () const |
| Return the solverType of this solver. More... | |
| MInt | restartInterval () const |
| Return the restart interval of this solver. More... | |
| MInt | restartTimeStep () const |
| Return the restart interval of this solver. More... | |
| MInt | solverId () const |
| Return the solverId. More... | |
| MBool | restartFile () |
| MInt | readSolverSamplingVarNames (std::vector< MString > &varNames, const MString featureName="") const |
| Read sampling variables names, store in vector and return the number of sampling variables. More... | |
| virtual MBool | hasRestartTimeStep () const |
| virtual MBool | forceAdaptation () |
| virtual void | preTimeStep ()=0 |
| virtual void | postTimeStep ()=0 |
| virtual void | initSolver ()=0 |
| virtual void | finalizeInitSolver ()=0 |
| virtual void | saveSolverSolution (const MBool NotUsed(forceOutput)=false, const MBool NotUsed(finalTimeStep)=false)=0 |
| virtual void | cleanUp ()=0 |
| virtual MBool | solutionStep () |
| virtual void | preSolutionStep (MInt) |
| virtual MBool | postSolutionStep () |
| virtual MBool | solverConverged () |
| virtual void | getInterpolatedVariables (MInt, const MFloat *, MFloat *) |
| virtual void | loadRestartFile () |
| virtual MInt | determineRestartTimeStep () const |
| virtual void | writeRestartFile (MBool) |
| virtual void | writeRestartFile (const MBool, const MBool, const MString, MInt *) |
| virtual void | setTimeStep () |
| virtual void | implicitTimeStep () |
| virtual void | prepareNextTimeStep () |
| Public Member Functions inherited from StructuredPostprocessing< nDim, FvStructuredSolver< nDim > > | |
| StructuredPostprocessing () | |
| Constructor for the postprocessing solver. More... | |
| ~StructuredPostprocessing () | |
| Destructor for the postprocessing solver. More... | |
| void | postprocessPreInit () |
| void | postprocessPreSolve () |
| void | postprocessPostSolve () |
| void | postprocessInSolve () |
Public Attributes | |
| std::vector< std::unique_ptr< StructuredZonalBC > > | m_zonalBC |
| void(FvStructuredSolver3D::* | reconstructSurfaceData )() |
| void(FvStructuredSolver3D::* | Venkatakrishnan_function )(MFloat, MFloat, MFloat, MInt, MInt, MFloatScratchSpace &) |
| void(FvStructuredSolver3D::* | viscFluxMethod )() |
| Public Attributes inherited from FvStructuredSolver< 3 > | |
| StructuredGrid< nDim > * | m_grid |
| MPI_Comm | m_StructuredComm |
| MInt | m_restartTimeStep |
| MString | m_outputFormat |
| MInt | m_lastOutputTimeStep |
| Public Attributes inherited from Solver | |
| std::set< MInt > | m_freeIndices |
| MBool | m_singleAdaptation = false |
| MBool | m_splitAdaptation = true |
| MBool | m_saveSensorData = false |
| Public Attributes inherited from StructuredPostprocessing< nDim, FvStructuredSolver< nDim > > | |
| MInt | m_restartTimeStep |
Protected Member Functions | |
| void | crossProduct (MFloat *, const MFloat *, const MFloat *) |
| MInt | cellIndex (const MInt, const MInt, const MInt) |
| MInt | pointIndex (const MInt, const MInt, const MInt) |
| MFloat | dist (MFloat *a, MFloat *b) |
| MInt | getPointIdFromCell (const MInt, const MInt, const MInt) |
| MInt | getPointIdfromPoint (const MInt, const MInt, const MInt, const MInt) |
| MInt | getCellIdfromCell (const MInt, const MInt, const MInt, const MInt) |
| MInt | surfId (MInt point, MInt isd, MInt dim) |
| void | manualInterpolationCorrection () |
| Manually correct errors made by the restart interpolation. More... | |
| void | interpolateFromDonor () |
| Interpolates the flow field from a given donor file. More... | |
| void | loadSampleFile (MString) override |
| Loads primitive variables from an HDF5 file. More... | |
| void | getSampleVariables (MInt, MFloat *) override |
| MFloat | getSampleVorticity (MInt, MInt) override |
| void | loadAverageRestartFile (const MChar *, MFloat **, MFloat **, MFloat **, MFloat **) override |
| Loads the postprocessing restart file to continue the postprocessing. More... | |
| void | loadAveragedVariables (const MChar *) override |
| Loads the averaged variables again to do further postprocessing. More... | |
| void | shiftAverageCellValuesRestart () |
| void | shiftAverageCellValues () |
| Shifts the averaged variables. More... | |
| MFloat | dvardxyz (MInt, MInt, MFloat *) override |
| MFloat | dvardx (MInt, MFloat *) override |
| MFloat | getCellLengthY (MInt, MInt, MInt) override |
| MFloat | getCellCoordinate (MInt, MInt) override |
| Protected Member Functions inherited from FvStructuredSolver< 3 > | |
| void | initTimers () |
| virtual void | initFsc () |
| Init for Falkner-Skan-Cooke flow. More... | |
| virtual MFloat | getFscPressure (MInt cellId) |
| virtual MFloat | getFscPressure (MFloat coordX) |
| virtual MFloat | getFscEta (MFloat coordX, MFloat coordY) |
| virtual void | getFscVelocity (MInt cellId, MFloat *const vel) |
| Load variables for the specified timeStep. More... | |
| virtual void | getFscVelocity (MFloat coordX, MFloat coordY, MFloat *const vel) |
| virtual void | initBlasius () |
| Init for Blasius boundary layer. More... | |
| virtual MFloat | getBlasiusEta (MFloat coordX, MFloat coordY) |
| virtual void | getBlasiusVelocity (MInt cellId, MFloat *const vel) |
| Load variables for the specified timeStep. More... | |
| virtual void | getBlasiusVelocity (MFloat coordX, MFloat coordY, MFloat *const vel) |
| Protected Member Functions inherited from Solver | |
| Solver (const MInt solverId, const MPI_Comm comm, const MBool isActive=true) | |
| MFloat | returnLoadRecord () const |
| MFloat | returnIdleRecord () const |
| Protected Member Functions inherited from StructuredPostprocessing< nDim, FvStructuredSolver< nDim > > | |
| void | initStructuredPostprocessing () |
| void | initAverageIn () |
| Initializes properties for averaging during solver run. More... | |
| void | initAverageVariables () |
| allocates memory for averageSolutions() and averageSolutionsInSolve() More... | |
| void | initTimeStepProperties () |
| Initializes timestep properties for postprocessing. More... | |
| void | initMovingAverage () |
| void | initProductionVariables () |
| void | initDissipationVariables () |
| void | averageSolutionsInSolve () |
| void | averageSolutions () |
| void | addAveragingSample () |
| Adds one sample to the summedVars. More... | |
| void | addTempWaveSample () |
| Adds for the travelling wave setups. More... | |
| void | saveAveragedSolution (MInt) |
| void | computeAveragedSolution () |
| Computes the mean variables from summed vars. More... | |
| void | computeAverageSkinFriction () |
| Computes skin friction of an averaged field. More... | |
| void | subtractPeriodicFluctuations () |
| void | subtractMean () |
| void | movingAverage () |
| void | movingAveragePost () |
| void | computeProductionTerms () |
| Computes the production terms from an averaged field. More... | |
| void | computeDissipationTerms () |
| Computes the production terms from an averaged field. More... | |
| void | decomposeCf () |
| void | decomposeCfDouble () |
| void | writeGradients () |
| void | loadAveragedSolution () |
| Loads an averaged file again. More... | |
| void | saveAverageRestart () |
| void | loadMeanFile (const MChar *fileName) |
| void | getSampleVariables (MInt cellId, const MFloat *&vars) |
| MInt | getNoPPVars () |
| Returns number of postprocessing variables. More... | |
| MInt | getNoVars () |
| MInt | getNoPPSquareVars () |
| Returns number of pp Square variables. More... | |
Static Protected Attributes | |
| static constexpr MInt | xsd = 0 |
| static constexpr MInt | ysd = 1 |
| static constexpr MInt | zsd = 2 |
| static constexpr const MInt | nDim = 3 |
| Static Protected Attributes inherited from StructuredPostprocessing< nDim, FvStructuredSolver< nDim > > | |
| static const MInt | xsd |
| static const MInt | ysd |
| static const MInt | zsd |
Friends | |
| template<MBool isRans> | |
| class | StructuredBndryCnd3D |
| class | FvStructuredSolver3DRans |
Additional Inherited Members | |
| Protected Types inherited from StructuredPostprocessing< nDim, FvStructuredSolver< nDim > > | |
| typedef void(StructuredPostprocessing::* | tpost) () |
| typedef std::vector< tpost > | tvecpost |
Detailed Description
Definition at line 32 of file fvstructuredsolver3d.h.
Member Typedef Documentation
◆ fluxmethod
| typedef void(FvStructuredSolver3D::* FvStructuredSolver3D::fluxmethod) (MFloat *, MFloat *, MInt, MInt) |
Definition at line 91 of file fvstructuredsolver3d.h.
Constructor & Destructor Documentation
◆ FvStructuredSolver3D()
| FvStructuredSolver3D::FvStructuredSolver3D | ( | MInt | solverId, |
| StructuredGrid< 3 > * | grid_, | ||
| MBool * | propertiesGroups, | ||
| const MPI_Comm | comm | ||
| ) |
Definition at line 17 of file fvstructuredsolver3d.cpp.
20 TRACE();
22
23 // count the no of necessary FQ fields and allocate
24 initializeFQField();
25
26 // compute the cell center coordinates from point coordinates
27 m_grid->computeCellCenterCoordinates();
28
31 } else {
33 }
34
35 allocateSingularities();
36
37 // allocate memory for aux data maps (cf,cp)
38 allocateAuxDataMaps();
39
40 // assign coordinates to all ghost points
42
43 // if we are RANS we should allocate a new RANS solver
46 }
47
48 initFluxMethod();
49
50 m_grid->computeCellCenterCoordinates();
52 m_grid->computeMetrics();
55 m_grid->computeJacobian();
57
58 // TODO_SS labels:FV,totest By now I am not sure if Code performs correctly for wrongly oriented meshes
63 if(m_cells->cellJac[cellId] < 0.0) mTerm(1, "Negative Jacobian found! Check first if code can cope this!");
64 }
65 }
66 }
67
69
71 computeZonalConnections();
72 }
73
74 // Assign handlers to the correct boundary conditions
75 assignBndryCells();
76
77 // Computation of modified wall distance in porous computation requires to set the porosity,
78 // Da-number etc. first; on the other hand we need to wait for initializeFQField to be called
80 initPorous();
81 // exchange6002();
82 }
83
86 m_structuredBndryCnd->computeWallDistances();
87 }
88 }
89
91 initInterpolatedPoints();
92 }
93
95 initPointsToAsciiFile();
96 }
97
99
101}
void initInterpolatedPoints()
Compute all the interpolation coefficients necessary for the interpolation output.
Definition: fvstructuredsolver3d.cpp:2063
void allocateSingularities()
Definition: fvstructuredsolver3d.cpp:9711
void initPointsToAsciiFile() override
Definition: fvstructuredsolver3d.cpp:2077
MInt cellIndex(const MInt, const MInt, const MInt)
Definition: fvstructuredsolver3d.cpp:4998
void addGhostPointCoordinateValues()
Extrapolates and exchanges ghost point coordinates.
Definition: fvstructuredsolver3d.cpp:4583
void assignBndryCells()
Definition: fvstructuredsolver3d.cpp:4569
std::unique_ptr< FvStructuredSolver3DRans > m_ransSolver
Definition: fvstructuredsolver3d.h:203
std::unique_ptr< StructuredBndryCnd< 3 > > m_structuredBndryCnd
Definition: fvstructuredsolver3d.h:183
void computeZonalConnections()
Definition: fvstructuredsolver3d.cpp:105
std::array< MInt, Timers::_count > m_timers
Definition: fvstructuredsolver.h:572
MInt m_pointsToAsciiOutputInterval
Definition: fvstructuredsolver.h:537
MInt m_intpPointsOutputInterval
Definition: fvstructuredsolver.h:493
MPI_Comm m_StructuredComm
Definition: fvstructuredsolver.h:324
void initializeFQField()
Counts the number of necessary FQ fields, allocates them and corrects the indexes of the FQ variable ...
Definition: fvstructuredsolver.cpp:277
void allocateAuxDataMaps()
AUX DATA //////////////////////////////////////////////////////////////////.
Definition: fvstructuredsolver.cpp:7804
void mTerm(const MInt errorCode, const MString &location, const MString &message)
Definition: functions.cpp:29
void printAllocatedMemory(const MLong oldAllocatedBytes, const MString &solverName, const MPI_Comm comm)
Prints currently allocated memory.
Definition: genericobject.cpp:26
◆ ~FvStructuredSolver3D()
| FvStructuredSolver3D::~FvStructuredSolver3D | ( | ) |
Definition at line 103 of file fvstructuredsolver3d.cpp.
103{ TRACE(); }
Member Function Documentation
◆ addDisturbance()
| void FvStructuredSolver3D::addDisturbance | ( | ) |
Definition at line 4524 of file fvstructuredsolver3d.cpp.
4524 {
4525 TRACE();
4526 cout << "enterin addDisturbance " << endl;
4531 MFloat pi2 = 8.0 * atan(1);
4532 MFloat period = 0.1;
4534
4535 // MFloat T=(m_cells->variables[CV->RHO_E][364])*(m_gamma-F1)/287.1500;
4536 // MFloat fluc_p = rhsdist*287.15000*T;
4537 // MFloat fluc_rhoE= (F1/(m_gamma-1))*fluc_p;
4538 // rhsdist=fluc_rhoE;
4540 cout << "sin " << amp1 << endl;
4541 MFloat centerloc = 3.1415;
4542 // go through all cells and adopt disturbance smoothly!!
4543 // m_cells->rightHandSide[CV->RHO_E][2191941]+=rhsdist;
4551 MFloat r = 0.025;
4557 /* MFloat x=m_cells->coordinates[0][cellId]-centreloc;
4558 MFloat y=m_cells->coordinates[1][cellId]-centreloc;
4559 MFloat z=m_cells->coordinates[2][cellId]-centreloc;
4560 m_cells->rightHandSide[CV->RHO_E][cellId]+=(F1+tanh(x/0.009))*(F1-tanh(x/0.009))*(F1+tanh(y/0.009))*(F1-tanh(y/0.009))*(F1+tanh(z/0.009))*(F1-tanh(z/0.009))*rhsdist;*/
4561 }
4562 }
4563 }
4564 }
4565
4566 cout << "leaving addDisturbance " << endl;
4567}
std::unique_ptr< MConservativeVariables< nDim > > CV
Definition: fvstructuredsolver.h:332
◆ addGhostPointCoordinateValues()
| void FvStructuredSolver3D::addGhostPointCoordinateValues | ( | ) |
Definition at line 4583 of file fvstructuredsolver3d.cpp.
4583 {
4584 TRACE();
4585 // > for debugging only: save cellId/pointId and domainId for all cells and points
4593 }
4594 }
4595 }
4596 }
4597 //< end debugging only
4598
4599 m_grid->extrapolateGhostPointCoordinates();
4600
4603
4605
4606 m_grid->computeCellCenterCoordinates();
4607
4608 // MUST be done after cell center computation!!!
4611
4614 }
4615
4617 m_grid->writePartitionedGrid();
4618 }
4619}
void extrapolateGhostPointCoordinatesBC()
Definition: fvstructuredsolver3d.cpp:4621
void computeReconstructionConstantsSVD()
Definition: fvstructuredsolver3d.cpp:9941
std::unique_ptr< StructuredFQVariables > FQ
Definition: fvstructuredsolver.h:469
MBool m_savePartitionOutput
Definition: fvstructuredsolver.h:669
std::vector< std::unique_ptr< StructuredComm< nDim > > > m_sndComm
Definition: fvstructuredsolver.h:728
std::vector< std::unique_ptr< StructuredComm< nDim > > > m_rcvComm
Definition: fvstructuredsolver.h:729
◆ allocateSingularities()
| void FvStructuredSolver3D::allocateSingularities | ( | ) |
Definition at line 9711 of file fvstructuredsolver3d.cpp.
9711 {
9716
9721 }
9722
9723 // 4 unknowns and 2*Nstar cells
9724 mAlloc(m_singularity[i].ReconstructionConstants, nDim + 1, m_singularity[i].totalCells * m_singularity[i].Nstar * 2,
9725 "ReconstructionConstants", 0.0, AT_);
9726 }
9727}
void mAlloc(T *&a, const MLong N, const MString &objectName, MString function)
allocates memory for one-dimensional array 'a' of size N
Definition: alloc.h:173
SingularInformation * m_singularity
Definition: fvstructuredsolver.h:835
◆ applyBodyForce()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 7863 of file fvstructuredsolver3d.cpp.
7863 {
7864 TRACE();
7867
7869 case 10: {
7870 // traveling wave case
7872 if(zeroPos) {
7873 t_offset = F0;
7874 }
7877
7878 MFloat fadeInFactor = 0.0;
7881 } else {
7882 fadeInFactor = 1.0;
7883 }
7884
7885 if(zeroPos) {
7886 fadeInFactor = F1;
7887 }
7888
7896
7897 MFloat transitionFactor = F0;
7899 transitionFactor = F0;
7903 transitionFactor = F1;
7905 transitionFactor = (1 + cos((x - m_waveOutBeginTransition) / transitionOutLength * pi)) * F1B2;
7906 } else {
7907 transitionFactor = F0;
7908 }
7909
7913
7918 }
7919 }
7920 }
7921
7922 break;
7923 }
7924 case 11: {
7925 // traveling wave case with predefined force field
7927 if(zeroPos) {
7928 t_offset = F0;
7929 }
7932
7933 MFloat fadeInFactor = 0.0;
7936 } else {
7937 fadeInFactor = 1.0;
7938 }
7939
7940 if(zeroPos) {
7941 fadeInFactor = F1;
7942 }
7943
7948
7951
7952 MFloat transitionFactor = F0;
7954 transitionFactor = F0;
7958 transitionFactor = F1;
7960 transitionFactor = (1 + cos((x - m_waveOutBeginTransition) / transitionOutLength * pi)) * F1B2;
7961 } else {
7962 transitionFactor = F0;
7963 }
7964
7967
7969 MFloat amp = 0.0;
7970
7978 amp = fstart + (pos - zzstart) / (zz - zzstart) * (f - fstart);
7979 } else if(pos > m_waveForceZ[jj + (m_grid->getMyBlockNoCells(0) - 1) * m_grid->getMyBlockNoCells(1)]) {
7986 amp = f + (pos - zz) / (zzend - zz) * (fend - f);
7987 } else {
7995 if(zz <= pos && pos < zzp) {
7996 amp = f + (pos - zz) / (zzp - zz) * (fp - f);
7997 break;
7998 }
7999 }
8000 }
8001
8003
8008 }
8009 }
8010 }
8011
8012 break;
8013 }
8014 default: {
8016 }
8017 }
8018}
MFloat m_waveEndTransition
Definition: fvstructuredsolver.h:367
MFloat m_waveDomainWidth
Definition: fvstructuredsolver.h:389
MFloat m_waveBeginTransition
Definition: fvstructuredsolver.h:366
MFloat * m_waveForceField
Definition: fvstructuredsolver.h:385
MFloat m_waveTemporalTransition
Definition: fvstructuredsolver.h:377
MInt m_bodyForceMethod
Definition: fvstructuredsolver.h:383
MFloat m_waveOutEndTransition
Definition: fvstructuredsolver.h:369
MFloat m_wavePenetrationHeight
Definition: fvstructuredsolver.h:382
MFloat m_movingGridTimeOffset
Definition: fvstructuredsolver.h:350
MFloat m_waveOutBeginTransition
Definition: fvstructuredsolver.h:368
◆ applyBoundaryCondition()
|
overridevirtual |
Implements FvStructuredSolver< 3 >.
Definition at line 2274 of file fvstructuredsolver3d.cpp.
2274 {
2275 TRACE();
2276 // treat Dirichlet and Neumann BC in one go!!!
2277 m_structuredBndryCnd->applyDirichletNeumannBC();
2278}
◆ applyInviscidBoundaryCondition()
| void FvStructuredSolver3D::applyInviscidBoundaryCondition | ( | ) |
Definition at line 9369 of file fvstructuredsolver3d.cpp.
9369{ TRACE(); }
◆ applySandpaperTrip()
| void FvStructuredSolver3D::applySandpaperTrip | ( | ) |
Definition at line 4077 of file fvstructuredsolver3d.cpp.
4077 {
4078 TRACE();
4079
4087
4089 m_tripTimeStep[ii] = tripTime;
4090
4091 // copy old values from H2 to H1
4094 }
4095
4096 // also copy the old mode coefficients
4099 }
4100
4101 // compute new fourier coefficients
4102 tripFourierCoefficients(&m_tripModesH2[offsetModes], m_tripNoModes, m_tripDomainWidth, m_tripCutoffZ[ii]);
4103 tripForceCoefficients(&m_tripModesH2[offsetModes], &m_tripH2[offset], m_tripCoords, m_tripNoCells, m_tripNoModes);
4104 }
4105
4110
4116
4117 MFloat force = 0.0;
4118 if(x > m_tripXOrigin[ii] - 2 * m_tripXLength[ii] && x < m_tripXOrigin[ii] + 2 * m_tripXLength[ii]
4119 && y > m_tripYOrigin[ii] - 2 * m_tripYHeight[ii] && y < m_tripYOrigin[ii] + 2 * m_tripYHeight[ii]) {
4122 * forceStrength;
4123 }
4124
4129 }
4130 }
4131 }
4132 }
4133}
void tripForceCoefficients(MFloat *, MFloat *, MFloat *, MInt, MInt)
Definition: fvstructuredsolver3d.cpp:4214
void tripFourierCoefficients(MFloat *, MInt, MFloat, MFloat)
Definition: fvstructuredsolver3d.cpp:4246
MFloat * m_tripDeltaTime
Definition: fvstructuredsolver.h:850
MFloat * m_tripMaxAmpSteady
Definition: fvstructuredsolver.h:847
MFloat m_tripDomainWidth
Definition: fvstructuredsolver.h:861
MFloat * m_tripMaxAmpFluc
Definition: fvstructuredsolver.h:848
constexpr std::underlying_type< FcCell >::type p(const FcCell property)
Converts property name to underlying integer value.
Definition: fccellproperties.h:32
Definition: geometrycontexttypes.h:23
◆ applySandpaperTripAirfoil()
| void FvStructuredSolver3D::applySandpaperTripAirfoil | ( | ) |
Definition at line 4135 of file fvstructuredsolver3d.cpp.
4135 {
4136 TRACE();
4137
4145
4147 m_tripTimeStep[ii] = tripTime;
4148
4149 // copy old values from H2 to H1
4152 }
4153
4154 // also copy the old mode coefficients
4157 }
4158
4159 // compute new fourier coefficients
4160 tripFourierCoefficients(&m_tripModesH2[offsetModes], m_tripNoModes, m_tripDomainWidth, m_tripCutoffZ[ii]);
4161 tripForceCoefficients(&m_tripModesH2[offsetModes], &m_tripH2[offset], m_tripCoords, m_tripNoCells, m_tripNoModes);
4162 }
4163
4168
4175
4176 MFloat force = 0.0;
4177 if(x > m_tripXOrigin[ii] - 2 * m_tripXLength[ii] && x < m_tripXOrigin[ii] + 2 * m_tripXLength[ii]
4178 && wallDist < 2 * m_tripYHeight[ii]) {
4179 MFloat factor = 0.0;
4181 factor = exp(-POW2((x - m_tripXOrigin[ii]) / m_tripXLength[ii]) - POW2(wallDist / m_tripYHeight[ii]));
4183 factor = exp(-POW2((x - m_tripXOrigin[ii]) / m_tripXLength[ii]) - POW2(wallDist / m_tripYHeight[ii]));
4184 }
4185
4186 force = forceStrength * factor;
4187 }
4188
4194
4198
4201
4204
4207 }
4208 }
4209 }
4210 }
4211}
MFloat * m_airfoilNormalVec
Definition: fvstructuredsolver.h:866
MFloat * m_airfoilWallDist
Definition: fvstructuredsolver.h:868
std::unique_ptr< MPrimitiveVariables< nDim > > PV
Definition: fvstructuredsolver.h:333
MInt m_airfoilNoWallPoints
Definition: fvstructuredsolver.h:867
◆ applyViscousBoundaryCondition()
| void FvStructuredSolver3D::applyViscousBoundaryCondition | ( | ) |
Definition at line 9371 of file fvstructuredsolver3d.cpp.
9371{ TRACE(); }
◆ assignBndryCells()
| void FvStructuredSolver3D::assignBndryCells | ( | ) |
Definition at line 4569 of file fvstructuredsolver3d.cpp.
◆ Ausm()
| void FvStructuredSolver3D::Ausm | ( | ) |
Definition at line 3689 of file fvstructuredsolver3d.cpp.
3689 {
3690 // Ausm routines have been moved and are called from inside Muscl (better performance)
3691}
◆ AusmDV()
Definition at line 2863 of file fvstructuredsolver3d.cpp.
2863 {
2867
2871
2873
2874 // left side
2881
2882 // right side
2889
2894
2895
2896 // compute lenght of metric vector for normalization
2899
2900 // scale by metric length to get velocity in the new basis (get normalized basis vectors)
2902
2903
2905
2906 MFloat AL = FRHOL * PL;
2907 MFloat AR = FRHOR * PR;
2908
2912
2913 AL = sqrt(gamma * AL);
2914 AR = sqrt(gamma * AR);
2915 AL = mMax(AL, AR);
2916 AR = AL;
2917
2920
2921 AL = AL * DGRAD;
2922 AR = AR * DGRAD;
2923
2926
2929
2934 MFloat SV = 2.0 * DGRAD / (m_cells->cellJac[I] + m_cells->cellJac[IP1]) * (FDV + (F1 - FDV) * getPSI(I, dim));
2938
2941
2942 MFloat FXMA = F1B2 * (XMAL1 + fabs(XMAL1));
2944 FXMA = F1B2 * (XMAR1 - fabs(XMAR1));
2945 const MFloat XMARM = ALPHAR * (-F1B4 * POW2(XMAR1 - F1) - FXMA) + FXMA + (mMin(-F1, XMAR) + F1);
2946
2950
2956
2959 UL = UL - UUL2;
2960 UR = UR - UUR2;
2963 VL = VL - UUL3;
2964 VR = VR - UUR3;
2967 WL = WL - UUL4;
2968 WR = WR - UUR4;
2969
2971
2972 flux[CV->RHO_U][I] = RHOU2 * (UL + UR) + ARHOU2 * (UL - UR) + PLR * surf0 + RHOUL * UUL2 + RHOUR * UUR2;
2973 flux[CV->RHO_V][I] = RHOU2 * (VL + VR) + ARHOU2 * (VL - VR) + PLR * surf1 + RHOUL * UUL3 + RHOUR * UUR3;
2974 flux[CV->RHO_W][I] = RHOU2 * (WL + WR) + ARHOU2 * (WL - WR) + PLR * surf2 + RHOUL * UUL4 + RHOUR * UUR4;
2975 flux[CV->RHO_E][I] = RHOU2 * (e0 + e1) + ARHOU2 * (e0 - e1) + PLR * dxdtau;
2976 flux[CV->RHO][I] = RHOU;
2977}
◆ AusmLES()
|
inline |
Can be used for moving grids, dxt term is included
Definition at line 2692 of file fvstructuredsolver3d.cpp.
2693 {
2694 // MFloat pFactor[3]={F0,F0,F0};
2698
2702
2704
2705 // calculate pressure
2711
2717
2718 // compute lenght of metric vector for normalization
2721
2722 // scale by metric length to get velocity in the new basis (get normalized basis vectors)
2724
2725
2727
2728
2729 // speed of sound
2732
2735
2738
2741
2744
2747
2748 const MFloat RHOU = F1B2 * (MALR * (RHO_AL + RHO_AR) + fabs(MALR) * (RHO_AL - RHO_AR)) * metricLength;
2750 // multiply by metric length to take surface area into account
2752
2753 //==>fluxes:
2755 flux[CV->RHO_U][I] = RHOU2 * (UL + UR) + AbsRHO_U2 * (UL - UR) + PLR * surf0;
2756 flux[CV->RHO_V][I] = RHOU2 * (VL + VR) + AbsRHO_U2 * (VL - VR) + PLR * surf1;
2757 flux[CV->RHO_W][I] = RHOU2 * (WL + WR) + AbsRHO_U2 * (WL - WR) + PLR * surf2;
2758 flux[CV->RHO_E][I] = RHOU2 * (e0 + e1) + AbsRHO_U2 * (e0 - e1) + PLR * dxdtau;
2759 flux[CV->RHO][I] = RHOU;
2760}
◆ AusmLES_PTHRC()
|
inline |
Same AUSM scheme as AusmLES with additional damping controlled by the 4th order pressure derivative. Pressure needs to computed beforehand.
Definition at line 2771 of file fvstructuredsolver3d.cpp.
2771 {
2774
2778
2780
2782
2783 // calculate pressure
2789
2795
2796 // compute lenght of metric vector for normalization
2799
2800 // scale by metric length to get velocity in the new basis (get normalized basis vectors)
2802
2803
2805
2806 // speed of sound
2809
2812
2814
2815 // 4th order pressure damping
2820
2825
2830
2834
2838
2840
2843
2846
2849
2850 const MFloat RHOU = F1B2 * (MALR * (RHO_AL + RHO_AR) + fabs(MALR) * (RHO_AL - RHO_AR)) * metricLength;
2852 // multiply by metric length to take surface area into account
2854
2856 flux[CV->RHO_U][I] = RHOU2 * (UL + UR) + AbsRHO_U2 * (UL - UR) + PLR * surf0;
2857 flux[CV->RHO_V][I] = RHOU2 * (VL + VR) + AbsRHO_U2 * (VL - VR) + PLR * surf1;
2858 flux[CV->RHO_W][I] = RHOU2 * (WL + WR) + AbsRHO_U2 * (WL - WR) + PLR * surf2;
2859 flux[CV->RHO_E][I] = RHOU2 * (e0 + e1) + AbsRHO_U2 * (e0 - e1) + PLR * dxdtau;
2860 flux[CV->RHO][I] = RHOU;
2861}
MInt getCellIdfromCell(const MInt, const MInt, const MInt, const MInt)
Definition: fvstructuredsolver3d.cpp:5002
MFloat m_fgammaMinusOne
Definition: fvstructuredsolver.h:444
◆ BARTH_JESPERSON_FCT()
| void FvStructuredSolver3D::BARTH_JESPERSON_FCT | ( | MFloat | effNghbrDelta, |
| MFloat | srfcDelta, | ||
| MFloat | dxEpsSqr, | ||
| MInt | cellPos, | ||
| MInt | var, | ||
| MFloatScratchSpace & | minPhi | ||
| ) |
- Date
- 2015
Definition at line 2677 of file fvstructuredsolver3d.cpp.
◆ calcSurfaceMetrics()
| void FvStructuredSolver3D::calcSurfaceMetrics | ( | ) |
◆ cellIndex()
Definition at line 4998 of file fvstructuredsolver3d.cpp.
◆ computeCellLength()
| void FvStructuredSolver3D::computeCellLength | ( | ) |
Definition at line 866 of file fvstructuredsolver3d.cpp.
866 {
867 // this function can be moved into the MusclSchemeStreched later but for testing it is easier
868 // REMEMBER: FOR MOVINg GRIDS THIS NEEDS TO BE CALLED EACH TIME
869
870 MInt P1 = -1, P2 = -1, P3 = -1, P4 = -1, P5 = -1, P6 = -1, P7 = -1, P8 = -1;
871 MFloat f1x = F0, f1y = F0, f1z = F0, f2x = F0, f2y = F0, f2z = F0;
876 P1 = getPointIdFromCell(i, j, k);
877 P2 = getPointIdfromPoint(P1, 1, 0, 0);
878 P3 = getPointIdfromPoint(P1, 1, 1, 0);
879 P4 = getPointIdfromPoint(P1, 1, 0, 1);
880 P5 = getPointIdfromPoint(P1, 1, 1, 1);
881 P6 = getPointIdfromPoint(P1, 0, 1, 0);
882 P7 = getPointIdfromPoint(P1, 0, 1, 1);
883 P8 = getPointIdfromPoint(P1, 0, 0, 1);
884 //----------Idirection
885 // face 1
886 f1x = F1B4
888 + m_grid->m_coordinates[0][P8]);
889 f1y = F1B4
891 + m_grid->m_coordinates[1][P8]);
892 f1z = F1B4
894 + m_grid->m_coordinates[2][P8]);
895 // face 2
896 f2x = F1B4
898 + m_grid->m_coordinates[0][P5]);
899 f2y = F1B4
901 + m_grid->m_coordinates[1][P5]);
902 f2z = F1B4
904 + m_grid->m_coordinates[2][P5]);
906 //----------Jdirection
907 // face 1
908 f1x = F1B4
910 + m_grid->m_coordinates[0][P8]);
911 f1y = F1B4
913 + m_grid->m_coordinates[1][P8]);
914 f1z = F1B4
916 + m_grid->m_coordinates[2][P8]);
917 // face 2
918 f2x = F1B4
920 + m_grid->m_coordinates[0][P7]);
921 f2y = F1B4
923 + m_grid->m_coordinates[1][P7]);
924 f2z = F1B4
926 + m_grid->m_coordinates[2][P7]);
928 //----------Kdirection
929 // face 1
930 f1x = F1B4
932 + m_grid->m_coordinates[0][P6]);
933 f1y = F1B4
935 + m_grid->m_coordinates[1][P6]);
936 f1z = F1B4
938 + m_grid->m_coordinates[2][P6]);
939 // face 2
940 f2x = F1B4
942 + m_grid->m_coordinates[0][P8]);
943 f2y = F1B4
945 + m_grid->m_coordinates[1][P8]);
946 f2z = F1B4
948 + m_grid->m_coordinates[2][P8]);
950 }
951 }
952 }
953}
MInt getPointIdFromCell(const MInt, const MInt, const MInt)
Definition: fvstructuredsolver3d.cpp:5788
MInt getPointIdfromPoint(const MInt, const MInt, const MInt, const MInt)
Definition: fvstructuredsolver3d.cpp:5792
◆ computeCumulativeAverage()
|
overridevirtual |
do the time average of the flow variables for LES and store them.
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 4267 of file fvstructuredsolver3d.cpp.
4267 {
4268 TRACE();
4277 }
4278 } else {
4279 const MFloat timeFacRho = m_Ma * (1.0 / m_zonalAveragingFactor) * sqrt(PV->TInfinity) * m_physicalTimeStep;
4280 const MFloat timeFacU = m_Ma * (1.0 / m_zonalAveragingFactor) * sqrt(PV->TInfinity) * m_physicalTimeStep;
4281 const MFloat timeFacV = m_Ma * (1.0 / m_zonalAveragingFactor) * sqrt(PV->TInfinity) * m_physicalTimeStep;
4282 const MFloat timeFacW = m_Ma * (1.0 / m_zonalAveragingFactor) * sqrt(PV->TInfinity) * m_physicalTimeStep;
4283 const MFloat timeFacE = m_Ma * (1.0 / m_zonalAveragingFactor) * sqrt(PV->TInfinity) * m_physicalTimeStep;
4284
4287 // Exponential averaging
4289 timeFacRho * m_cells->pvariables[PV->RHO][cellId] + (1.0 - timeFacRho) * m_cells->fq[FQ->AVG_RHO][cellId];
4291 timeFacU * m_cells->pvariables[PV->U][cellId] + (1.0 - timeFacU) * m_cells->fq[FQ->AVG_U][cellId];
4293 timeFacV * m_cells->pvariables[PV->V][cellId] + (1.0 - timeFacV) * m_cells->fq[FQ->AVG_V][cellId];
4295 timeFacW * m_cells->pvariables[PV->W][cellId] + (1.0 - timeFacW) * m_cells->fq[FQ->AVG_W][cellId];
4297 timeFacE * m_cells->pvariables[PV->P][cellId] + (1.0 - timeFacE) * m_cells->fq[FQ->AVG_P][cellId];
4298 }
4299 }
4300 }
4301}
MFloat m_physicalTimeStep
Definition: fvstructuredsolver.h:565
MFloat m_zonalAveragingFactor
Definition: fvstructuredsolver.h:720
◆ computeDomainWidth()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 3693 of file fvstructuredsolver3d.cpp.
3693 {
3694 // if we only compute cd,cl for a part of the domain (m_auxDataCoordinateLimits == true)
3695 // only compute the average with the width of this section
3698 } else {
3699 MFloat minCoordinate = F0, maxCoordinate = F0, minCoordinateGlobal = F0, maxCoordinateGlobal = F0;
3700 MInt lowPoint = -1, highPoint = -1;
3701
3708 } else {
3709 highPoint = getPointIdFromCell(m_noGhostLayers, m_noGhostLayers, m_nCells[0] - m_noGhostLayers);
3711 }
3712
3715 MPI_Allreduce(&minCoordinate, &minCoordinateGlobal, 1, MPI_DOUBLE, MPI_MIN, m_StructuredComm, AT_, "minCoordinate",
3716 "minCoordinateGlobal");
3717 MPI_Allreduce(&maxCoordinate, &maxCoordinateGlobal, 1, MPI_DOUBLE, MPI_MAX, m_StructuredComm, AT_, "maxCoordinate",
3718 "maxCoordinateGlobal");
3719 m_globalDomainWidth = fabs(maxCoordinateGlobal - minCoordinateGlobal);
3721 }
3722}
MBool m_auxDataCoordinateLimits
Definition: fvstructuredsolver.h:489
MInt m_forceAveragingDir
Definition: fvstructuredsolver.h:478
MFloat * m_auxDataLimits
Definition: fvstructuredsolver.h:490
MFloat m_globalDomainWidth
Definition: fvstructuredsolver.h:488
int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allreduce
Definition: mpioverride.cpp:952
◆ computeFrictionPressureCoef()
|
inlineoverridevirtual |
Implements FvStructuredSolver< 3 >.
Definition at line 138 of file fvstructuredsolver3d.h.
◆ computeLambda2Criterion()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 9162 of file fvstructuredsolver3d.cpp.
9162 {
9163 TRACE();
9165 MFloat d[3] = {F0, F0, F0};
9166 MFloat e[3] = {F0, F0, F0};
9167 J.fill(F0);
9168
9169 // MInt IMJK, IJMK, IMJMK, IJKM, IJMKM, IMJKM;
9174
9175 // compute the lambda2 criterion.
9187
9200
9213
9226
9227 // Compute the matrix (S^2+Omega^2)/2
9228 J(0, 0) = POW2(dudx) + dvdx * dudy + dwdx * dudz;
9229 J(0, 1) = F1B2 * (dudx * dudy + dudx * dvdx + dvdy * dudy + dvdy * dvdx + dwdy * dudz + dvdz * dwdx);
9230 J(0, 2) = F1B2 * (dudx * dudz + dudx * dwdx + dudy * dvdz + dvdx * dwdy + dwdz * dudz + dwdz * dwdx);
9231 J(1, 0) = J(0, 1);
9232 J(1, 1) = POW2(dvdy) + dvdx * dudy + dvdz * dwdy;
9233 J(1, 2) = F1B2 * (dudz * dvdx + dwdx * dudy + dvdy * dvdz + dvdy * dwdy + dwdz * dvdz + dwdz * dwdy);
9234 J(2, 0) = J(0, 2);
9235 J(2, 1) = J(1, 2);
9236 J(2, 2) = POW2(dwdz) + dwdx * dudz + dvdz * dwdy;
9237
9238 // perform householder tridiagonalization of
9239 // symmetric real matrix
9241 // compute eigenvalues
9243 // sort eigenvalues
9245
9246 lambda2[cellId] = d[1];
9247 }
9248 }
9249 }
9250}
void tred2(MFloatScratchSpace &A, MInt dim, MFloat *diag, MFloat *offdiag)
Householder Reduction according to Numercial Recipies in C: The Art of Scientific Computing.
Definition: fvstructuredsolver.cpp:6859
void tqli2(MFloat *diag, MFloat *offdiag, MInt dim)
Compute Eigenvalues with implicit shift according to Numercial Recipies in C: The Art of Scientific C...
Definition: fvstructuredsolver.cpp:6934
void insertSort(MInt dim, MFloat *list)
Sorting function to sort list in ascending order.
Definition: fvstructuredsolver.cpp:6998
◆ computePorousRHS()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 5501 of file fvstructuredsolver3d.cpp.
5501 {
5502 TRACE();
5503
5505
5511
5512 if(isRans) {
5514 } else { // LES
5519 MFloat velAbs = F0;
5522 }
5523 velAbs = sqrt(velAbs);
5528 -(rRe0 * rDa * por[IJK] * muLam[IJK] + porPOW2 * sqrt(rDa) * cf[IJK] * velAbs * pvars[PV->RHO][IJK])
5529 * pvars[dim /*PV->VV[dim] <-- why is this giving linking error???*/][IJK] * m_cells->cellJac[IJK];
5530 }
5531 }
5532 }
5533 }
5534 }
5535}
◆ computePrimitiveVariables()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 9703 of file fvstructuredsolver3d.cpp.
9703 {
9705 computePrimitiveVariables_<&FvStructuredSolver::pressure_twoEqRans>();
9706 else
9707 computePrimitiveVariables_();
9708}
void computePrimitiveVariables_()
Definition: fvstructuredsolver3d.cpp:9669
◆ computePrimitiveVariables_()
| void FvStructuredSolver3D::computePrimitiveVariables_ |
Definition at line 9669 of file fvstructuredsolver3d.cpp.
9669 {
9671
9674
9680 MFloat velPOW2 = F0;
9683 velPOW2 += POW2(pvars[vel][cellId]);
9684 }
9685
9686 // density and pressure:
9689 gammaMinusOne
9690 * (cvars[CV->RHO_E][cellId] - F1B2 * pvars[PV->RHO][cellId] * velPOW2 + (this->*pressure_func)(cellId));
9691
9693 // TODO_SS labels:FV,totest Does it make sense to forbid negative values. BCs sometimes negate the neighbor
9694 // value
9697 }
9698 }
9699 }
9700 }
9701}
MInt m_noRansEquations
Definition: fvstructuredsolver.h:698
◆ computeReconstructionConstantsSVD()
| void FvStructuredSolver3D::computeReconstructionConstantsSVD | ( | ) |
Definition at line 9941 of file fvstructuredsolver3d.cpp.
9941 {
9942 MInt nghbr[30], dim = 0;
9944 m_orderOfReconstruction = 1;
9946 MInt maxNoSingularityRecNghbrIds = 14;
9950 MFloat counter = F0;
9951 MFloat avg = F0;
9952 MFloat maxc = F0;
9953
9956 MInt totalCells = 1;
9958
9959 //(p)reset the reconstruction constants
9963 }
9964 }
9965
9968 if(len1[j] != 0) totalCells *= len1[j];
9969 }
9970
9973 dim = n;
9976 } else {
9979 }
9980 }
9981
9985 MInt count = 0;
9987 temp[dim] = 1;
9988
9989 nghbr[count++] = cellIndex(ii, jj, kk);
9990 nghbr[count++] = cellIndex(ii + temp[0], jj + temp[1], kk + temp[2]);
9991
9992 // the coordinates of the corner where the viscousflux should be corrected.
9993 MInt ijk = getPointIdFromCell(ii + m_singularity[i].Viscous[0], jj + m_singularity[i].Viscous[1],
9994 kk + m_singularity[i].Viscous[2]);
9995 ijk = getPointIdfromPoint(ijk, 1, 1, 1);
9996
9999 nghbr[count++] = cellIndex(ii + change[0], jj + change[1], kk + change[2]);
10000 nghbr[count++] = cellIndex(ii + temp[0] + change[0], jj + temp[1] + change[1], kk + temp[2] + change[2]);
10001 }
10002
10004 cerr << "Something wrong with the singularities in the LS coeffiecient computation" << endl;
10005 }
10006
10007 // weighted Least square
10008 weights.fill(F0);
10009
10010 // Compute weights with RBF (take mean distance as R0)
10012 MInt nghbrId = nghbr[n];
10013 MFloat dxdx = F0;
10016 }
10017
10018 weights[n] = 1 / dxdx; // RBF( dxdx, POW2( dist) );
10019 }
10020
10021 MInt id2 = ii - start[0] + ((jj - start[1]) + (kk - start[2]) * len1[1]) * len1[0];
10023
10025 avg += condNum;
10026 maxc = mMax(maxc, condNum);
10027 counter += F1;
10028 if(condNum < F0 || condNum > 1e7 || std::isnan(condNum)) {
10030 << " with large condition number: " << condNum << " num of neighbor" << count << "x" << recDim << " "
10031 << " coords " << m_grid->m_coordinates[0][ijk] << ", " << m_grid->m_coordinates[1][ijk] << ", "
10032 << m_grid->m_coordinates[2][ijk] << endl;
10033 }
10034 }
10035 }
10036 }
10037 }
10038 }
10039}
MFloat computeRecConstSVD(const MInt ijk, const MInt noNghbrIds, MInt *nghbr, MInt ID, MInt sID, MFloatScratchSpace &tmpA, MFloatScratchSpace &tmpC, MFloatScratchSpace &weights, const MInt recDim)
AUX DATA ENDS /////////////////////////////////////////////////////////////.
Definition: fvstructuredsolver.cpp:8020
MInt m_orderOfReconstruction
Definition: fvstructuredsolver.h:435
MFloat ** ReconstructionConstants
Definition: structuredgrid.h:27
◆ computeTimeStep()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 4348 of file fvstructuredsolver3d.cpp.
4348 {
4349 TRACE();
4350
4351 m_timeStep = 1000.0;
4356
4361
4363
4364 // compute the speed of sound
4366
4367 // no need for simplified metrics, since information is already contained
4368 // in cell metrics
4369 const MFloat lenXi = sqrt(POW2(metric[0][cellId]) + POW2(metric[1][cellId]) + POW2(metric[2][cellId]));
4370 const MFloat lenEt = sqrt(POW2(metric[3][cellId]) + POW2(metric[4][cellId]) + POW2(metric[5][cellId]));
4371 const MFloat lenZe = sqrt(POW2(metric[6][cellId]) + POW2(metric[7][cellId]) + POW2(metric[8][cellId]));
4372
4373 // contravariant velocities
4374 MFloat U_c = F0;
4375 MFloat V_c = F0;
4376 MFloat W_c = F0;
4377
4382 }
4383
4384 // subtract grid velocity
4385 U_c -= dxtx[cellId];
4386 V_c -= dxty[cellId];
4387 W_c -= dxtz[cellId];
4388
4389 U_c = fabs(U_c);
4390 V_c = fabs(V_c);
4391 W_c = fabs(W_c);
4392
4393 // has area information in it due to metric terms
4395
4396 // divide volume information (jacobian) through area to get characteristic length for CFL
4398
4401 m_timeStep = F1;
4402 m_timeRef = F1;
4403 } else {
4405 }
4406 }
4407 }
4408 }
4409}
◆ computeTotalKineticEngergy()
| MFloat FvStructuredSolver3D::computeTotalKineticEngergy | ( | ) |
Definition at line 5749 of file fvstructuredsolver3d.cpp.
5749 {
5750 TRACE();
5751 MFloat localEnergy = F0;
5756 localEnergy += (POW2(m_cells->pvariables[PV->U][cellId]) + POW2(m_cells->pvariables[PV->V][cellId])
5759 }
5760 }
5761 }
5762 return localEnergy;
5763}
◆ computeTotalPressure()
| MFloat FvStructuredSolver3D::computeTotalPressure | ( | ) |
Definition at line 5765 of file fvstructuredsolver3d.cpp.
5765 {
5766 TRACE();
5767 MFloat localPressure = F0;
5773 }
5774 }
5775 }
5776 return localPressure;
5777}
◆ computeVorticity()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 9093 of file fvstructuredsolver3d.cpp.
9093 {
9094 TRACE();
9102
9103
9114
9118
9122
9126
9133
9140
9147
9148 vortx[IJK] = F1B2 * (dwdy - dvdz) / jac[IJK];
9149 vorty[IJK] = F1B2 * (dudz - dwdx) / jac[IJK];
9150 vortz[IJK] = F1B2 * (dvdx - dudy) / jac[IJK];
9151 }
9152 }
9153 }
9154}
◆ computeZonalConnections()
| void FvStructuredSolver3D::computeZonalConnections | ( | ) |
Definition at line 105 of file fvstructuredsolver3d.cpp.
105 {
107 // do a spanwise average of the pressure in
108 // all LES zones
110 }
111
116 m_zonalBC.resize(noZonalBCMaps);
117
122 cout << "///////////////////////////////////////////////////////////////////" << endl
123 << "Starting zonal bc creation for " << noZonalBCMaps << " zonal bc maps" << endl;
124 }
125
128 cout << "///////////////////////////////////////////////////////////////////" << endl;
129 cout << "//////////// ZONAL BC " << id << " ///////////////////////////////////////////" << endl;
130 cout << "///////////////////////////////////////////////////////////////////" << endl;
131 }
132
135
137 cout << "Zonal BC" << id << ": Initialization ..." << endl;
138 }
139
140 // set up the data of the zonal domain and noCellsBC
142 MInt hasRcvDomain = 0;
148
151 m_windowInfo->checkZonalBCMaps(m_windowInfo->m_zonalBCMaps[id], m_windowInfo->physicalBCMap[bcId]);
152
153 MInt stgIP = 0;
156 stgIP = 1; // in this case use one more cell, required for STG
158 }
159
162 m_zonalBC[id]->noCellsLocalBC += (end[0] + stgIP - start[0]) * (end[1] - start[1]) * (end[2] - start[2]);
169 hasRcvDomain = 1;
170 break;
171 }
172 }
173
175 cout << "Zonal BC" << id << ": Initialization ... Finished!" << endl;
176 cout << "Zonal BC" << id << ": Collecting local coordinates ... " << endl;
177 }
178
179 MFloatScratchSpace localCoordinatesBC(nDim, std::max(1, m_zonalBC[id]->noCellsLocalBC), AT_, "localCoordinatesBC");
180 MIntScratchSpace localReceiverIds(std::max(1, m_zonalBC[id]->noCellsLocalBC), AT_, "localReceiverIds");
181 MIntScratchSpace localMapCellIds(std::max(1, m_zonalBC[id]->noCellsLocalBC), AT_, "m_localMapCellIds");
182
184 MInt bcCellId = 0;
189
191 localMapCellIds[bcCellId] = bcCellId;
192 } else {
193 localMapCellIds[bcCellId] = cellId;
194 }
195
198 }
199
200 localReceiverIds[bcCellId] = domainId();
201 ++bcCellId;
202 }
203 }
204 }
205 }
206
208 cout << "Zonal BC" << id << ": Collecting local coordinates ... Finished!" << endl;
209 cout << "Zonal BC" << id << ": Exchanging information about local coordinates ..." << endl;
210 }
211
212 MPI_Allreduce(&m_zonalBC[id]->noCellsLocalBC, &m_zonalBC[id]->noCellsGlobalBC, 1, MPI_INT, MPI_SUM,
214 MPI_Allgather(&m_zonalBC[id]->noCellsLocalBC, 1, MPI_INT, &noZonalCells[0], 1, MPI_INT, m_StructuredComm, AT_,
215 "m_zonalBC[id]->noCellsLocalBC", "noZonalCells[0]");
216
218 cout << "Zonal BC" << id << ": Exchanging information about local coordinates ... Finished!" << endl;
219 }
220
225 MInt noCells = 0;
227 noCells += noZonalCells[i];
228 }
229 cout << "Zonal BC" << id << ": Gathering coordinates for " << noCells << " cells on all partitions..." << endl;
230 }
231
233 mAlloc(m_zonalBC[id]->coordinatesGlobalBC, 3, std::max(1, m_zonalBC[id]->noCellsGlobalBC), "coordinatesGlobalBC",
234 0.0, AT_);
235 mAlloc(m_zonalBC[id]->globalLocalMapCellIds, std::max(1, m_zonalBC[id]->noCellsGlobalBC), "globalLocalMapCellIds",
236 -1, AT_);
237 mAlloc(m_zonalBC[id]->globalReceiverIds, std::max(1, m_zonalBC[id]->noCellsGlobalBC), "globalReceiverIds", -1, AT_);
238
239 noCellsBCOffsets[0] = 0;
241 noCellsBCOffsets[i] = noCellsBCOffsets[i - 1] + noZonalCells[i - 1];
242 }
243
244
245 // Distribute all zonal cell information (all coordinates, cellIds, domainIds) to all domains
248 &m_zonalBC[id]->coordinatesGlobalBC[dim][0], &noZonalCells[0], &noCellsBCOffsets[0], MPI_DOUBLE,
250 }
251 MPI_Allgatherv(&localReceiverIds[0], noZonalCells[domainId()], MPI_INT, &m_zonalBC[id]->globalReceiverIds[0],
253 "m_zonalBC[id]->globalReceiverIds[0]");
254 MPI_Allgatherv(&localMapCellIds[0], noZonalCells[domainId()], MPI_INT, &m_zonalBC[id]->globalLocalMapCellIds[0],
256 "m_zonalBC[id]->globalLocalMapCellIds[0]");
257
259 MInt noCells = 0;
261 noCells += noZonalCells[i];
262 }
263 cout << "Zonal BC" << id << ": Gathering coordinates for " << noCells << " cells on all partitions... Finished!"
264 << endl;
265 cout << "Zonal BC" << id << ": Searching for suitable interpolation partners..." << endl;
266 }
272
274 hasInterpolationPartnerDomain = false;
275 }
276
279
282 hasPartnerLocalBC.begin(),
283 hasInterpolationPartnerDomain);
284
286 cout << "Zonal BC" << id << ": Searching for suitable interpolation partners... Finished!" << endl;
287 cout << "Zonal BC" << id << ": Creating sending and receiving information..." << endl;
288 }
289
293 // to find out which domains will be rcvDomains
294
296 MPI_Allreduce(&hasRcvDomain, &m_zonalBC[id]->noGlobalRcvDomains, 1, MPI_INT, MPI_SUM, m_StructuredComm, AT_,
297 "hasRcvDomain", "m_zonalBC[id]->noGlobalRcvDomains");
298
299 mAlloc(m_zonalBC[id]->globalRcvDomainIds, std::max(1, m_zonalBC[id]->noGlobalRcvDomains), "m_globalRcvDomainIds",
300 -1, AT_);
301
302 MInt pos = 0;
304 if(noZonalCells[j] > 0) {
306 ++pos;
307 }
308 }
309
310 // to find out which domains will be sndDomains
313 if(hasPartnerLocalBC[cellIdBC]) {
314 hasSndDomain = true;
315 break;
316 }
317 }
318
320 MPI_Allreduce(&hasSndDomain, &m_zonalBC[id]->noGlobalSndDomains, 1, MPI_INT, MPI_SUM, m_StructuredComm, AT_,
321 "hasSndDomain", "m_zonalBC[id]->noGlobalSndDomains");
322
323
325 MPI_Allgather(&hasSndDomain, 1, MPI_INT, &hasSndDomainInfo[0], 1, MPI_INT, m_StructuredComm, AT_, "hasSndDomain",
326 "hasSndDomainInfo[0]");
327
328 mAlloc(m_zonalBC[id]->globalSndDomainIds, std::max(1, m_zonalBC[id]->noGlobalSndDomains), "m_globalSndDomainIds",
329 -1, AT_);
330
331 pos = 0;
333 if(hasSndDomainInfo[j] > 0) {
335 pos++;
336 }
337 }
341 // compute no of cells which our domain can interpolate and has to send
344 if(hasPartnerLocalBC[cellIdBC]) {
346 }
347 }
348
349 mAlloc(m_zonalBC[id]->localCommReceiverIds, std::max(1, m_zonalBC[id]->noBufferSndSize), "localCommReceiverIds",
350 AT_);
351 mAlloc(m_zonalBC[id]->localBufferMapCellIds, std::max(1, m_zonalBC[id]->noBufferSndSize), "localCommReceiverIds",
352 AT_);
353 mAlloc(m_zonalBC[id]->localBufferIndexCellIds, std::max(1, m_zonalBC[id]->noBufferSndSize), "localCommReceiverIds",
354 AT_);
355
356 // collect the receiver domainIds for each cell that our domain can interpolate
357 pos = 0;
359 if(hasPartnerLocalBC[cellIdBC]) {
361 pos++;
362 }
363 }
364
365 // count to how many different domains our domain needs to send data
373 break;
374 }
375 }
376 }
377 }
378 }
379
381 MIntScratchSpace localBufferSndSize(std::max(1, m_zonalBC[id]->noSndNghbrDomains), AT_, "localBufferSndSize");
382
383 // save each receiving domainId in localRcvId
384 pos = 0;
390 localRcvId[pos] = d;
391 pos++;
392 break;
393 }
394 }
395 }
396 }
397 }
398
399
400 // save the cellId and mapCellId (the cellId of the receiving cell on the nghbr Domain)
401 pos = 0;
403 if(hasPartnerLocalBC[cellIdBC]) {
406 m_zonalBC[id]->localBufferIndexCellIds[pos] = cellIdBC; // this is needed to get InterpolatedVars in zonalGather
407 pos++;
408 }
409 }
410
411 // compute the buffer size for each domain our domain needs to send data to
415 localBufferSndSize[i]++;
416 }
417 }
418 }
419
420 // get m_zonalBC[id]->noRcvNghbrDomains and SndId
422
424 MInt bufferRcvSize = 0;
427 if(hasPartnerLocalBC[cellIdBC]) {
428 bufferRcvSize++;
429 }
430 }
431 }
432
433 // get the bufferSndSize
435 }
436
438 MPI_Alltoall(&sndBufferRcvSize[0], 1, MPI_INT, &rcvBufferRcvSize[0], 1, MPI_INT, m_StructuredComm, AT_,
439 "sndBufferRcvSize[0]", "rcvBufferRcvSize[0]");
440
441 // get m_zonalBC[id]->noRcvNghbrDomains, localSndId, localbufferRcvSize.
442 // These are from which domains each rcvDoamin will receive
443 // and how much size it is.
447 if(rcvBufferRcvSize[i] > 0) { // get noRcvNghbrDoamins
449 }
450 }
451 }
452
454 MIntScratchSpace localBufferRcvSize(std::max(1, m_zonalBC[id]->noRcvNghbrDomains), AT_, "localBufferRcvSize");
455
456 // get localSndId from which domains each rcvDoamin will receive data
457 pos = 0;
460 if(rcvBufferRcvSize[i] > 0) { // sndId is the info the receiver Domains have from which domains I will receive
461 localSndId[pos] = i;
462 pos++;
463 }
464 }
465 }
466 // get localbufferRcvSize
467 pos = 0;
470 if(rcvBufferRcvSize[i] > 0) {
471 localBufferRcvSize[pos] = rcvBufferRcvSize[i];
472 pos++;
473 }
474 }
475 }
476
480
481 mAlloc(m_zonalBC[id]->sndRequest, std::max(1, m_zonalBC[id]->noSndNghbrDomains), "sndRequest", AT_);
482 mAlloc(m_zonalBC[id]->sndStatus, std::max(1, m_zonalBC[id]->noSndNghbrDomains), "sndStatus", AT_);
483 mAlloc(m_zonalBC[id]->rcvRequest, std::max(1, m_zonalBC[id]->noRcvNghbrDomains), "rcvRequest", AT_);
484 mAlloc(m_zonalBC[id]->rcvStatus, std::max(1, m_zonalBC[id]->noRcvNghbrDomains), "rcvStatus", AT_);
485
486
488 unique_ptr<StructuredZonalComm> sndCommPtr =
489 make_unique<StructuredZonalComm>(localBufferSndSize[i], localRcvId[i], m_zonalBC[id]->noZonalVariables);
492 }
493
495 unique_ptr<StructuredZonalComm> rcvCommPtr =
496 make_unique<StructuredZonalComm>(localBufferRcvSize[i], localSndId[i], m_zonalBC[id]->noZonalVariables);
499 }
500
501
505
507 pos = 0;
508
512 pos++;
513 }
514 }
515
516 MInt err = MPI_Isend((void*)m_zonalBC[id]->sndComm[noSnd]->mapCellId, m_zonalBC[id]->sndComm[noSnd]->bufferSize,
519 if(err)
521 << " throws error " << endl;
522 }
523
525 MInt err =
526 MPI_Irecv((void*)&m_zonalBC[id]->rcvComm[noRcv]->mapCellId[0], m_zonalBC[id]->rcvComm[noRcv]->bufferSize,
529 if(err)
531 << " throws error " << endl;
532 }
533
534 MPI_Waitall(m_zonalBC[id]->noSndNghbrDomains, m_zonalBC[id]->sndRequest, m_zonalBC[id]->sndStatus, AT_);
535 MPI_Waitall(m_zonalBC[id]->noRcvNghbrDomains, m_zonalBC[id]->rcvRequest, m_zonalBC[id]->rcvStatus, AT_);
536
538
540 cout << "Zonal BC" << id << ": Creating sending and receiving information... Finished!" << endl;
541 cout << "Zonal BC" << id << ": Completed!" << endl;
542 }
543 }
544}
std::vector< std::unique_ptr< StructuredZonalBC > > m_zonalBC
Definition: fvstructuredsolver3d.h:69
std::unique_ptr< FvStructuredSolverWindowInfo< nDim > > m_windowInfo
Definition: fvstructuredsolver.h:419
std::vector< MFloat * > m_zonalSpanwiseAvgVars
Definition: fvstructuredsolver.h:703
int MPI_Isend(const void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Isend
Definition: mpioverride.cpp:481
int MPI_Barrier(MPI_Comm comm, const MString &name)
same as MPI_Barrier
Definition: mpioverride.cpp:874
int MPI_Allgatherv(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, const int recvcounts[], const int displs[], MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allgatherv
Definition: mpioverride.cpp:1315
int MPI_Irecv(void *buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Irecv
Definition: mpioverride.cpp:606
int MPI_Waitall(int count, MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Waitall
Definition: mpioverride.cpp:766
int MPI_Alltoall(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Alltoall
Definition: mpioverride.cpp:1356
int MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allgather
Definition: mpioverride.cpp:1274
◆ crossProduct()
◆ dist()
◆ distributeFluxToCells()
| void FvStructuredSolver3D::distributeFluxToCells | ( | ) |
Definition at line 4265 of file fvstructuredsolver3d.cpp.
4265{ TRACE(); }
◆ dvardx()
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 9435 of file fvstructuredsolver3d.cpp.
9435 {
9437
9444
9448
9450 (dvardxi * m_cells->cellMetrics[xsd * 3 + xsd][IJK] + dvardet * m_cells->cellMetrics[ysd * 3 + xsd][IJK]
9452
9454}
◆ dvardxyz()
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 9414 of file fvstructuredsolver3d.cpp.
9414 {
9416
9423
9427
9429 (dvardxi * m_cells->cellMetrics[xsd * 3 + dir][IJK] + dvardet * m_cells->cellMetrics[ysd * 3 + dir][IJK]
9431
9433}
◆ exchange6002()
| void FvStructuredSolver3D::exchange6002 | ( | ) |
GATHER & SEND
RECEIVE
WAIT
SCATTER
Definition at line 10043 of file fvstructuredsolver3d.cpp.
10043 {
10045 // if (m_rans && m_ransMethod!=RANS_KEPSILON)
10046 // mTerm(1, "Porous RANS computation is only supported by k-epsilon model!");
10047 // if (!m_porous)
10048 // mTerm(1, "bc6002 requires the property porous to be set to true!");
10049
10050 // 0) Check if initBc6002 is called for the first time
10051 // for(MInt bcId_=0; bcId_ < bcId; ++bcId_) {
10052 // if (m_physicalBCMap[bcId_]->BC==6002) return;
10053 // }
10054
10055 // Determine normal vectors and save for later use
10058 && m_structuredBndryCnd->m_physicalBCMap[bcId_]->Nstar == -1) {
10061
10064
10066 {0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1},
10067 {0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0}};
10068
10070
10079 const MInt incp[nDim] = {IJKP[normalDir], IJKP[firstTangentialDir], IJKP[secondTangentialDir]};
10081
10085
10089
10090 // compute four surrounding points of surface centroid
10092 const MInt pp1 = getPointIdfromPoint(ijk, pp[normalDir][0], pp[normalDir][1], pp[normalDir][2]);
10093 const MInt pp2 = getPointIdfromPoint(ijk, pp[normalDir][3], pp[normalDir][4], pp[normalDir][5]);
10094 const MInt pp3 = getPointIdfromPoint(ijk, pp[normalDir][6], pp[normalDir][7], pp[normalDir][8]);
10096 helper[normalDir] = 1;
10098 getPointIdfromPoint(ijk, n * helper[0], n * helper[1], n * helper[2]); // point lying outside domain
10099
10100 // compute the velocity of the surface centroid
10109 }
10110
10111 // compute normal vector of surface
10112 crossProduct(normalVec, firstVec, secondVec);
10113 const MFloat normalLength = sqrt(POW2(normalVec[0]) + POW2(normalVec[1]) + POW2(normalVec[2]));
10114
10115 MFloat sgn = (std::inner_product(&normalVec[0], &normalVec[0] + nDim, &normalVec_[0], 0.0) < 0.0) ? -1 : 1;
10117
10119 normalVec[dim] /= normalLength;
10121 }
10122 }
10123 }
10124 }
10125 }
10126
10127 // Determine normal vector at singularities
10131 // only correct for bc 6000 not for bc 4000-5000
10132 if(singularity.BC == -6000) {
10135 if(singularity.BCsingular[n] == -6002) takeIt = true;
10136 }
10137
10138 if(takeIt) {
10141 if(singularity.end[n] - singularity.start[n] > 1) {
10142 // dim=n;
10143 // start[n]=singularity.start[n]+1;
10144 start[n] = singularity.start[n] + 1;
10145 end[n] = singularity.end[n] - 1;
10146 } else {
10147 start[n] = singularity.start[n];
10148 end[n] = singularity.end[n];
10149 }
10150 }
10151
10162 }
10163 }
10164 }
10165 }
10166 }
10167 }
10168 }
10169
10170 // TODO_SS labels:FV,toenhance The exchange of all the normals is an overhead, because it is only needed at
10171 // singularity points
10175 // MInt sendSizeTotal = 0;
10176 std::vector<MInt> receiveSizes;
10178 // TODO_SS labels:FV right now exchange at all 6000er not only 6002 because of singularities
10179 if(snd->bcId == -6000 || snd->bcId == -6002) {
10180 // Gather
10181
10182 MInt size = 1;
10184 size *= snd->endInfoCells[dim] - snd->startInfoCells[dim];
10185 // std::vector<MFloat> snd->cellBuffer(size*(1+nDim));
10186 // sendSizeTotal += size;
10187
10188 MInt pos = 0;
10193 // TODO_SS labels:FV Latter only allocate FQ->POROSITY for m_blockType==porous
10194 // if (m_blockType=="fluid")
10195 // snd->cellBuffer[pos++] = 1;
10196 // else
10200 }
10201 ++pos;
10202 }
10203 }
10204 }
10205
10206 // Send
10208 MInt err = MPI_Isend((void*)&snd->cellBuffer[0], size * (nDim + 1), MPI_DOUBLE, snd->nghbrId, tag,
10211
10212
10213 // Determine size of receive buffer
10214 // TODO_SS labels:FV right now exchange at all 6000er not only 6002 because of singularities
10215 size = 1;
10217 size *= snd->endInfoCells[dim] - snd->startInfoCells[dim];
10218 receiveSizes.push_back(size);
10219 }
10220 }
10221
10225 std::vector<MFloat> bufferRcv(std::accumulate(receiveSizes.begin(), receiveSizes.end(), 0) * (nDim + 1));
10226 MInt offset = 0;
10227 MInt cnt = 0;
10229 // TODO_SS labels:FV right now exchange at all 6000er not only 6002 because of singularities
10230 if(rcv->bcId == -6002 || rcv->bcId == -6000) {
10233 MInt err = MPI_Irecv(&bufferRcv[offset], rcvSize * (1 + nDim), MPI_DOUBLE, rcv->nghbrId, tag, m_StructuredComm,
10234 &rcv->mpi_request, AT_, "(void*)&rcvSize");
10236
10237 offset += rcvSize * (1 + nDim);
10238 ++cnt;
10239 }
10240 }
10241
10246 if(snd->bcId == -6002 || snd->bcId == -6000) {
10247 MPI_Wait(&(snd->mpi_request), &(snd->mpi_status), AT_);
10248 }
10249 }
10250
10251
10253 if(rcv->bcId == -6002 || rcv->bcId == -6000) {
10254 MPI_Wait(&(rcv->mpi_request), &(rcv->mpi_status), AT_);
10255 }
10256 }
10257
10258
10263 offset = 0;
10264 cnt = 0;
10266 if(rcv->bcId == -6002 || rcv->bcId == -6000) {
10267 // TODO_SS labels:FV right now exchange at all 6000er not only 6002 because of singularities
10268
10269 MInt k2, j2, i2, id2;
10275 MInt totalCells = 1;
10277
10279 len1[j] = rcv->endInfoCells[j] - rcv->startInfoCells[j];
10280 if(len1[j] != 0) totalCells *= len1[j];
10281 // added check the step for RCV part !!!!!!!!important
10282 step2[rcv->orderInfo[j]] = rcv->stepInfo[j];
10283 }
10284
10285 // Sanity check
10286 ASSERT(totalCells == receiveSizes[cnt], "");
10287
10289 start2[j] = 0;
10290 end2[j] = len1[j] - 1;
10291 len2[rcv->orderInfo[j]] = len1[j];
10292 if(step2[j] < 0) {
10294 start2[j] = end2[j];
10295 end2[j] = dummy;
10296 }
10297 }
10298
10299 MInt pos = 0;
10300 k2 = start2[2];
10302 j2 = start2[1];
10304 i2 = start2[0];
10306 start1[rcv->orderInfo[0]] = i2;
10307 start1[rcv->orderInfo[1]] = j2;
10308 start1[rcv->orderInfo[2]] = k2;
10309
10310 id2 = start1[0] + (start1[1] + start1[2] * len2[1]) * len2[0];
10314 normals_temp[m_noCells * d + cellId] = bufferRcv[offset * (nDim + 1) + (d + 1) * totalCells + id2];
10315 // m_cells->fq[FQ->NORMAL[d]][cellId] = bufferRcv[offset*(nDim+1)+(d+1)*totalCells+id2];
10316 }
10317
10318 i2 += step2[0];
10319 pos++;
10320 }
10321 j2 += step2[1];
10322 }
10323 k2 += step2[2];
10324 }
10325
10326 offset += totalCells;
10327 ++cnt;
10328 }
10329 }
10330
10331 // Determine normal vector at singularities
10334 // only correct for bc 6000 not for bc 4000-5000
10335 if(singularity.BC == -6000) {
10338 if(singularity.BCsingular[n] == -6002) takeIt = true;
10339 }
10340
10341 if(takeIt) {
10344 if(singularity.end[n] - singularity.start[n] > 1) {
10345 // dim=n;
10346 // start[n]=singularity.start[n]+1;
10347 start[n] = singularity.start[n] + 1;
10348 end[n] = singularity.end[n] - 1;
10349 } else {
10350 start[n] = singularity.start[n];
10351 end[n] = singularity.end[n];
10352 }
10353 }
10354
10356
10364 }
10369 temp[d] += normals_temp[m_noCells * d + nghbr];
10370 }
10371
10372 MFloat l = 0;
10376 }
10377 l = sqrt(l);
10380 }
10381
10387 }
10388 }
10389 }
10390 }
10391 }
10392 }
10393 }
10394
10396#if 0
10397 // 1) Gather
10398 std::vector<MFloat> sendBuffer;
10399 std::vector<MInt> sendcounts;
10400 std::vector<MInt> snghbrs;
10401 std::vector<MInt> tags;
10403 if (snd->bcId==-6002) {
10404 snghbrs.push_back(snd->nghbrId);
10405 tags.push_back(domainId()+(snd->tagHelper)*m_solver->noDomains());
10406
10407 MInt size = 1;
10409 size *= snd->endInfoCells[dim] - snd->startInfoCells[dim];
10410 sendcounts.push_back(size);
10411
10412 sendBuffer.resize(pos+size);
10416 // TODO_SS labels:FV Latter only allocate FQ->POROSITY for m_blockType==porous
10417// if (m_blockType=="fluid")
10418// sendBuffer[pos++] = 1;
10419// else
10421 }
10422 }
10423 }
10424 }
10425
10427
10428 // 2) Send & receive
10429 std::vector<MInt> recvcounts(noNeighborDomains);
10430 std::vector<MInt> rdispls(noNeighborDomains);
10431 std::vector<MFloat> recvBuffer = maia::mpi::mpiExchangePointToPoint(&sendBuffer[0],
10432 &snghbrs[0],
10433 noNeighborDomains,
10434 &sendcounts[0],
10435 &snghbrs[0],
10436 noNeighborDomains,
10437 m_StructuredComm,
10438 m_solver->domainId(),
10439 1,
10440 recvcounts.data(),
10441 rdispls.data());
10442 // 2.5) Send & receive the tags
10443 std::vector<MInt> sendcounts2(noNeighborDomains, 1);
10444 std::vector<MInt> recvTags = maia::mpi::mpiExchangePointToPoint(&tags[0],
10445 &snghbrs[0],
10446 noNeighborDomains,
10447 &sendcounts2[0],
10448 &snghbrs[0],
10449 noNeighborDomains,
10450 m_StructuredComm,
10451 m_solver->domainId(),
10452 1);
10453
10454 // 3) Scatter
10456 if (rcv->bcId==-6002) {
10458 MInt n;
10459 for (n = 0; n < noNeighborDomains; ++n) {
10460 if (tag==recvTags[n])
10461 break;
10462 }
10464
10467
10469 MInt j2, i2, id2;
10475 MInt totalCells=1;
10477
10479 len1[j]=rcv->endInfoCells[j] - rcv->startInfoCells[j];
10480 if(len1[j]!=0) totalCells*=len1[j];
10481 //added check the step for RCV part !!!!!!!!important
10482 step2[rcv->orderInfo[j]]=rcv->stepInfo[j];
10483 }
10484
10485 //TODO_SS labels:FV,totest check if this assert makes sense
10486 ASSERT(noReceivedElements==totalCells, "noReceivedElements == totalCells");
10487
10489 start2[j]=0;
10490 end2[j]=len1[j]-1;
10491 len2[rcv->orderInfo[j]]=len1[j];
10492 if(step2[j]<0) {
10494 start2[j]=end2[j];
10495 end2[j]=dummy;
10496 }
10497 }
10498
10499
10500 MInt pos=0;
10501 j2=start2[1];
10503 i2=start2[0];
10505 start1[rcv->orderInfo[0]]=i2;
10506 start1[rcv->orderInfo[1]]=j2;
10507
10508 id2=start1[0]+start1[1]*len2[0];
10511
10512 i2+=step2[0];
10513 pos++;
10514 }
10515 j2+=step2[1];
10516 }
10517 }
10518 }
10519#endif
10521}
void crossProduct(MFloat *, const MFloat *, const MFloat *)
Definition: fvstructuredsolver3d.cpp:5782
int MPI_Wait(MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Wait
Definition: mpioverride.cpp:732
std::vector< T > mpiExchangePointToPoint(const T *const sendBuffer, const MInt *const snghbrs, const MInt nosnghbrs, const MInt *const sendcounts, const MInt *const rnghbrs, const MInt nornghbrs, const MPI_Comm &mpi_comm, const MInt domainId, const MInt dataBlockSize, MInt *const recvcounts_=nullptr, MInt *const rdispls_=nullptr)
Definition: mpiexchange.h:1057
◆ extrapolateGhostPointCoordinatesBC()
| void FvStructuredSolver3D::extrapolateGhostPointCoordinatesBC | ( | ) |
Definition at line 4621 of file fvstructuredsolver3d.cpp.
4621 {
4623 // all the periodic BCs are NOT included.
4624 // also skip the channel bc
4626 || m_structuredBndryCnd->m_physicalBCMap[bcId]->BC == 2402
4627 || (m_structuredBndryCnd->m_physicalBCMap[bcId]->BC >= 6000
4628 && m_structuredBndryCnd->m_physicalBCMap[bcId]->BC < 6010)) {
4629 continue;
4630 }
4631
4636 MInt pos[3], fix[3], mirror[3], ijk[3], extendijk[3];
4637 MInt pointId, FixPointId, MirrorPointId;
4638
4639 extendijk[0] = 1;
4640 extendijk[1] = 1;
4641 extendijk[2] = 1;
4642 extendijk[index] = 0;
4643
4644 for(ijk[2] = start[2]; ijk[2] < end[2] + extendijk[2]; ++ijk[2]) {
4645 for(ijk[1] = start[1]; ijk[1] < end[1] + extendijk[1]; ++ijk[1]) {
4646 for(ijk[0] = start[0]; ijk[0] < end[0] + extendijk[0]; ++ijk[0]) {
4648 if(index == m) {
4649 if(step == 1) {
4650 pos[m] = ijk[m] + 1;
4651 fix[m] = start[m];
4652 mirror[m] = 2 * fix[m] - pos[m];
4653 } else {
4654 pos[m] = ijk[m];
4655 fix[m] = end[m];
4656 mirror[m] = 2 * fix[m] - pos[m];
4657 }
4658 } else {
4659 pos[m] = ijk[m];
4660 fix[m] = ijk[m];
4661 mirror[m] = ijk[m];
4662 }
4663 } // m
4664
4665 pointId = pointIndex(pos[0], pos[1], pos[2]);
4666 FixPointId = pointIndex(fix[0], fix[1], fix[2]);
4667 MirrorPointId = pointIndex(mirror[0], mirror[1], mirror[2]);
4668
4670 m_grid->m_coordinates[dim][pointId] =
4672 }
4673 } // ijk
4674 }
4675 }
4676 } // bcid
4677}
MInt pointIndex(const MInt, const MInt, const MInt)
Definition: fvstructuredsolver3d.cpp:5007
IdType index(const FloatType *const x, const IdType level)
Return Hilbert index for given location and level in 2D or 3D.
Definition: hilbert.h:165
◆ gather()
|
override |
Definition at line 4684 of file fvstructuredsolver3d.cpp.
4685 {
4686 for(auto& snd : sndComm) {
4689
4690 std::array<MInt, nDim> begin{snd->startInfoCells[0], snd->startInfoCells[1], snd->startInfoCells[2]};
4691 std::array<MInt, nDim> end{snd->endInfoCells[0], snd->endInfoCells[1], snd->endInfoCells[2]};
4692 std::array<MInt, nDim> size{snd->endInfoCells[0] - snd->startInfoCells[0],
4693 snd->endInfoCells[1] - snd->startInfoCells[1],
4694 snd->endInfoCells[2] - snd->startInfoCells[2]};
4696
4697 // Aliasing the unique_pointer in snd to a raw point is needed for PSTL on NVHPC
4698 auto* cellBuffer = snd->cellBuffer.get();
4699 auto* variables = &(snd->variables[0]);
4703 const MInt bufferId = totalSize * var + (i - begin[0]) + ((j - begin[1]) + (k - begin[2]) * size[1]) * size[0];
4704 cellBuffer[bufferId] = variables[var][cellId];
4705 });
4706 }
4707 }
4708}
MBool isPeriodicComm(std::unique_ptr< StructuredComm< nDim > > &)
Definition: fvstructuredsolver.cpp:8197
MBool skipPeriodicDirection(std::unique_ptr< StructuredComm< nDim > > &)
Definition: fvstructuredsolver.cpp:8191
◆ gcExtrapolate()
| void FvStructuredSolver3D::gcExtrapolate | ( | std::vector< MFloat * > & | variables | ) |
Definition at line 4910 of file fvstructuredsolver3d.cpp.
4910 {
4911 TRACE();
4912 // i-direction
4921
4923 variables[varPos][cellId] = (F2 * variables[varPos][cellIdAdj1] - variables[varPos][cellIdAdj2]);
4924 }
4925
4929
4931 variables[varPos][cellId] = (F2 * variables[varPos][cellIdAdj1] - variables[varPos][cellIdAdj2]);
4932 }
4933 }
4934 }
4935 }
4936
4937 // j-direction
4944
4946 variables[varPos][cellId] = (F2 * variables[varPos][cellIdAdj1] - variables[varPos][cellIdAdj2]);
4947 }
4948
4952
4954 variables[varPos][cellId] = (F2 * variables[varPos][cellIdAdj1] - variables[varPos][cellIdAdj2]);
4955 }
4956 }
4957 }
4958 }
4959
4960 // k-direction
4967
4969 variables[varPos][cellId] = (F2 * variables[varPos][cellIdAdj1] - variables[varPos][cellIdAdj2]);
4970 }
4971
4975
4977 variables[varPos][cellId] = (F2 * variables[varPos][cellIdAdj1] - variables[varPos][cellIdAdj2]);
4978 }
4979 }
4980 }
4981 }
4982}
◆ gcFillGhostCells()
|
virtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 4898 of file fvstructuredsolver3d.cpp.
4898 {
4899 gcExtrapolate(variables);
4900
4901 std::vector<std::unique_ptr<StructuredComm<nDim>>> gcSndComm;
4902 std::vector<std::unique_ptr<StructuredComm<nDim>>> gcRcvComm;
4903
4905
4906 m_windowInfo->createCommunicationExchangeFlags(gcSndComm, gcRcvComm, (MInt)variables.size(), varPtr);
4907 exchange(gcSndComm, gcRcvComm);
4908}
void gcExtrapolate(std::vector< MFloat * > &)
Definition: fvstructuredsolver3d.cpp:4910
void exchange()
SVD STUFF ENDS /////////////////////////////////////////////////////////////.
Definition: fvstructuredsolver.cpp:8061
◆ getCellCoordinate()
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 975 of file fvstructuredsolver3d.cpp.
◆ getCellIdfromCell()
|
inlineprotected |
Definition at line 5002 of file fvstructuredsolver3d.cpp.
◆ getCellLengthY()
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 955 of file fvstructuredsolver3d.cpp.
◆ getFourierCoefficients()
| void FvStructuredSolver3D::getFourierCoefficients | ( | MFloat * | k, |
| MFloat | k0, | ||
| std::complex< MFloat > * | fourierCoefficient | ||
| ) |
Original Implementation by Georg Eitel Amor for a given wavenumber k, and a certain energy spectrum (see Appendix of Orszag, 1969)
Definition at line 9880 of file fvstructuredsolver3d.cpp.
9880 {
9881 TRACE();
9882 MFloat r[6], s[6], kAbs, energy;
9883 complex<MFloat> uHat, vHat, wHat;
9884 // complex<MFloat>* fourierCoefficient;
9885
9886 kAbs = sqrt(k[0] * k[0] + k[1] * k[1] + k[2] * k[2]);
9887
9888 // the zero-frequency component is always set to zero, so there is no offset
9890 // fourierCoefficient = new complex<MFloat>[3];
9891
9892 fourierCoefficient[0] = complex<MFloat>(0, 0);
9893 fourierCoefficient[1] = complex<MFloat>(0, 0);
9894 fourierCoefficient[2] = complex<MFloat>(0, 0);
9895 } else {
9896 // energy = (kAbs/k0)*(kAbs/k0)*(kAbs/k0)*(kAbs/k0) * exp(-2.0*(kAbs/k0)*(kAbs/k0));
9897 // energy = pow(kAbs/k0,8.0) * exp(-4.0*(kAbs/k0)*(kAbs/k0)); // set spectral distribution
9898 energy = pow(kAbs / k0, 4.0) * exp(-2.0 * (kAbs / k0) * (kAbs / k0)); // set spectral distribution
9899 energy *= exp(2.0) * 0.499 * (m_Ma * LBCS)
9900 * (m_Ma * LBCS); // set maximal fluctuation amplitude to 20% of the freestream velocity (for 128^3: 0.88)
9901
9902 // determine Fourier coefficients:
9903 // r and s are Independant random vector fields with independant
9904 // components (zero mean and rms according to energy spectrum).
9905 // Each vector has three components for k and another three for -k.
9906
9908 r[i] = randnormal(0.0, PI * sqrt(energy) / (SQRT2 * kAbs));
9909 // r[i] = randNumGen.randNorm(0.0, PI*sqrt(energy)/(SQRT2*kAbs));
9910 s[i] = randnormal(0.0, PI * sqrt(energy) / (SQRT2 * kAbs));
9911 // s[i] = randNumGen.randNorm(0.0, PI*sqrt(energy)/(SQRT2*kAbs));
9912 }
9913
9914 uHat = (1.0 - k[0] * k[0] / (kAbs * kAbs)) * complex<MFloat>(r[0] + r[3], s[0] - s[3])
9915 - k[0] * k[1] / (kAbs * kAbs) * complex<MFloat>(r[1] + r[4], s[1] - s[4])
9916 - k[0] * k[2] / (kAbs * kAbs) * complex<MFloat>(r[2] + r[5], s[2] - s[5]);
9917
9918
9919 vHat = -k[1] * k[0] / (kAbs * kAbs) * complex<MFloat>(r[0] + r[3], s[0] - s[3])
9920 + (1.0 - k[1] * k[1] / (kAbs * kAbs)) * complex<MFloat>(r[1] + r[4], s[1] - s[4])
9921 - k[1] * k[2] / (kAbs * kAbs) * complex<MFloat>(r[2] + r[5], s[2] - s[5]);
9922
9923
9924 wHat = -k[2] * k[0] / (kAbs * kAbs) * complex<MFloat>(r[0] + r[3], s[0] - s[3])
9925 - k[2] * k[1] / (kAbs * kAbs) * complex<MFloat>(r[1] + r[4], s[1] - s[4])
9926 + (1.0 - k[2] * k[2] / (kAbs * kAbs)) * complex<MFloat>(r[2] + r[5], s[2] - s[5]);
9927
9928 // fourierCoefficient = new complex<MFloat>[3];
9929
9930 // fourierCoefficient[0] = complex<MFloat>(sqrt(2*energy)/SQRT2,sqrt(2*energy)/SQRT2);// uHat;
9931 // fourierCoefficient[1] = complex<MFloat>(sqrt(2*energy)/SQRT2,sqrt(2*energy)/SQRT2);//vHat;
9932 // fourierCoefficient[2] = complex<MFloat>(sqrt(2*energy)/SQRT2,sqrt(2*energy)/SQRT2);//wHat;
9933
9934 fourierCoefficient[0] = uHat;
9935 fourierCoefficient[1] = vHat;
9936 fourierCoefficient[2] = wHat;
9937 // return fourierCoefficient;
9938 }
9939}
MFloat randnormal(MFloat mu, MFloat sigma)
Returns a normal distributed random-number with mu=mean and sigma=standard deviation.
Definition: fvstructuredsolver3d.cpp:2233
◆ getPointIdFromCell()
|
inlineprotected |
Definition at line 5788 of file fvstructuredsolver3d.cpp.
◆ getPointIdfromPoint()
|
inlineprotected |
Definition at line 5792 of file fvstructuredsolver3d.cpp.
◆ getPSI()
Definition at line 10524 of file fvstructuredsolver3d.cpp.
◆ getSampleVariables()
|
overrideprotectedvirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 9373 of file fvstructuredsolver3d.cpp.
9373 {
9379}
◆ getSampleVorticity()
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 9381 of file fvstructuredsolver3d.cpp.
◆ initBndryCnds()
| void FvStructuredSolver3D::initBndryCnds | ( | ) |
Definition at line 4574 of file fvstructuredsolver3d.cpp.
◆ initBodyForce()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 8044 of file fvstructuredsolver3d.cpp.
◆ initFFTW()
| void FvStructuredSolver3D::initFFTW | ( | fftw_complex * | uPhysField, |
| fftw_complex * | vPhysField, | ||
| fftw_complex * | wPhysField, | ||
| MInt | lx, | ||
| MInt | ly, | ||
| MInt | lz, | ||
| MInt | noPeakModes | ||
| ) |
Definition at line 9732 of file fvstructuredsolver3d.cpp.
9733 {
9734 TRACE();
9735 // first check for odd numbers of cells in each direction
9736 if(lx % 2 != 0 || ly % 2 != 0 || lz % 2 != 0) {
9737 stringstream errorMessage;
9738 errorMessage << " FFTInit: Domain size must NOT be an odd number!: (lx)x(ly)x(lz)-> " << lx << "x" << ly << "x"
9739 << lz << endl;
9740 mTerm(1, AT_, errorMessage.str());
9741 }
9742
9745
9746 MFloat waveVector[3], k0;
9747
9748 complex<MFloat>* fourierCoefficient = new complex<MFloat>[3];
9749
9750 fftw_complex *uHatField, *vHatField, *wHatField;
9751
9752 fftw_plan planU, planV, planW;
9753
9754 // 1) Allocation of Fourier coefficients
9755 uHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
9756 vHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
9757 wHatField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
9758
9759 // 2) Creation of the plans for the FFTW
9760 planU = fftw_plan_dft_3d(lx, ly, lz, uHatField, uPhysField, FFTW_BACKWARD, FFTW_MEASURE);
9761 planV = fftw_plan_dft_3d(lx, ly, lz, vHatField, vPhysField, FFTW_BACKWARD, FFTW_MEASURE);
9762 planW = fftw_plan_dft_3d(lx, ly, lz, wHatField, wPhysField, FFTW_BACKWARD, FFTW_MEASURE);
9763
9768 // u-component
9769 uHatField[cellId][0] = 0.0;
9770 uHatField[cellId][1] = 0.0;
9771 uPhysField[cellId][0] = 0.0;
9772 uPhysField[cellId][1] = 0.0;
9773
9774 // v-component
9775 vHatField[cellId][0] = 0.0;
9776 vHatField[cellId][1] = 0.0;
9777 vPhysField[cellId][0] = 0.0;
9778 vPhysField[cellId][1] = 0.0;
9779
9780 // w-component
9781 wHatField[cellId][0] = 0.0;
9782 wHatField[cellId][1] = 0.0;
9783 wPhysField[cellId][0] = 0.0;
9784 wPhysField[cellId][1] = 0.0;
9785 }
9786 }
9787 }
9788
9789 // IMPORTANT NOTICE ON USE OF FFTW:
9790 // comment from Georg Eitel Amor:
9791 // FFTW stores the coefficients for positive wavenumbers in the first half of the array,
9792 // and those for negative wavenumbers in reverse order in the second half.
9793 // [0, 1, ... , N/2-1, N/2, ... , N-1]
9794 // - the entry for zero-wavenumber is at position 0
9795 // - the k-th entry and the (N-k)th entry correspond to wavenumbers with opposite sign
9796 // - the entry at position N/2 corresponds to the Nyquist wavenumber and appears only once
9797
9798 // peak wave number of energy spectrum
9799 k0 = 2.0 * PI / (lx / noPeakModes);
9800
9804 // wave-vector: k(p,q,r) = (2 \pi p / lx, 2 \pi q / ly, 2 \pi r / lz)
9805 waveVector[0] = (p)*2.0 * PI / lx;
9806 waveVector[1] = (q)*2.0 * PI / ly;
9807 waveVector[2] = (r)*2.0 * PI / lz;
9808
9809 getFourierCoefficients(waveVector, k0, fourierCoefficient);
9810
9811 // 1. Positive frequencies:
9812 uHatField[r + lz * (q + ly * p)][0] = real(fourierCoefficient[0]);
9813 uHatField[r + lz * (q + ly * p)][1] = imag(fourierCoefficient[0]);
9814
9815 vHatField[r + lz * (q + ly * p)][0] = real(fourierCoefficient[1]);
9816 vHatField[r + lz * (q + ly * p)][1] = imag(fourierCoefficient[1]);
9817
9818 wHatField[r + lz * (q + ly * p)][0] = real(fourierCoefficient[2]);
9819 wHatField[r + lz * (q + ly * p)][1] = imag(fourierCoefficient[2]);
9820
9821 // 2. Negative frequencies:
9822 if(p > 1 && q > 1 && r > 1) {
9823 if(p < lx / 2 && q < ly / 2 && r < lz / 2) {
9824 // since the physical velocity field is real, the coefficients for negative frequencies
9825 // are the complex conjugate of those for positive frequencies
9826 uHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = uHatField[r + lz * (q + ly * p)][0];
9827 uHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -uHatField[r + lz * (q + ly * p)][1];
9828
9829 vHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = vHatField[r + lz * (q + ly * p)][0];
9830 vHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -vHatField[r + lz * (q + ly * p)][1];
9831
9832 wHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][0] = wHatField[r + lz * (q + ly * p)][0];
9833 wHatField[(lz - r) + lz * ((ly - q) + ly * (lx - p))][1] = -wHatField[r + lz * (q + ly * p)][1];
9834 }
9835 }
9836 }
9837 }
9838 }
9839
9840 // Do Fourier transform (backward, see plan definition)
9841 // Definition in one dimension:
9842 // u(x) = \sum_{j=0}^{lx-1} \hat{u}_j exp(i 2 \pi j x / lx)
9843
9844 fftw_execute(planU);
9845 fftw_execute(planV);
9846 fftw_execute(planW);
9847
9848
9849 // normalize (this preserves the norm of the basis functions)
9856
9860 }
9861 }
9862 }
9863
9864 fftw_destroy_plan(planU);
9865 fftw_destroy_plan(planV);
9866 fftw_destroy_plan(planW);
9867 fftw_free(uHatField);
9868 fftw_free(vHatField);
9869 fftw_free(wHatField);
9870}
void getFourierCoefficients(MFloat *k, MFloat k0, std::complex< MFloat > *fourierCoefficient)
Generates a single complex coefficient of Fourier series.
Definition: fvstructuredsolver3d.cpp:9880
◆ initFluxMethod()
| void FvStructuredSolver3D::initFluxMethod | ( | ) |
Definition at line 674 of file fvstructuredsolver3d.cpp.
◆ initialCondition()
|
overridevirtual |
Initialization ot the entire flow field
- Date
- 16.02.2011
Implements FvStructuredSolver< 3 >.
Definition at line 1112 of file fvstructuredsolver3d.cpp.
1112 {
1113 TRACE();
1115 MFloat UT;
1116 MFloat pressureCH = F0;
1117
1127
1128 // compute conservative variables
1137
1138 // internal Reynolds number Re0 = Re / ( rho8*M*sqrt(T8)/T8^F072)
1140
1141 // reference enthalpies (needed for combustion computations)
1143
1144 // reference time (convection time)
1146
1147 m_deltaP = F0;
1148 // pressure loss per unit length dp = rho_00 u_tau^2 L / D ) here: D=1.0, L=1;
1149 // channel: dp = rho_00 u_tau^2 L / D )
1150 // m_deltaP = POW2( m_Ma * m_ReTau * sqrt(PV->TInfinity) / m_Re ) * CV->rhoInfinity / m_referenceLength;x
1151 // result is obtained by making deltap dimensionless with a_0^2 and rho_0
1152 // pipe: dp = lambda * L/D * rho/2 * u^2, lambda = 0.3164 Re^(-1/4) (Blasius)
1153
1157 CV->ransInfinity[0] = chi * (lamVisc);
1159 }
1160
1175
1177 cout << "////////////////////////////////////////////////" << endl;
1178 cout << "////////// Initial Condition summary ///////////" << endl;
1179 cout << "////////////////////////////////////////////////" << endl;
1193 }
1194
1196 // inflow condition
1197 // ----------------
1199 case 0: {
1200 // parallel inflow field
1202 // go through every cell
1206 }
1207
1209
1212 }
1213 }
1214
1215 break;
1216 }
1217 case 43: {
1218 // parallel inflow field with pressure peak in the middle of the domain
1220 // go through every cell
1224 }
1225
1227
1228 MFloat radius =
1229 sqrt(POW2(m_cells->coordinates[0][cellid] - 0.5) + POW2(m_cells->coordinates[1][cellid] - 0.5));
1230
1231 // impose pressure peak in the middle of the domain
1232 if(radius <= 0.05) {
1233 MFloat pAmp = 0.005;
1236 }
1237 }
1238 break;
1239 }
1240 case 333: {
1241 // parallel inflow field
1243 // go through every cell
1247 }
1248
1250
1251 // impose pressure peak in the middle of the domain
1253 MFloat pAmp = 0.005;
1257 }
1258 }
1259 break;
1260 }
1261 case 314: // stagnating flow field
1262 {
1267 }
1268
1270 }
1271 cout << "I.C. stagnating flow field was applied! " << endl;
1272 break;
1273 }
1274 case 315: // Poiseuille flow
1275 {
1278
1282
1285 }
1286 // all velocity components are 0 except u, Poiseuille distribution:
1289
1290 // the pressure is defined through the axial pressure gradient:
1292 PV->PInfinity - F3 * (x + 15.0) * SUTHERLANDLAW(PV->TInfinity) * PV->UInfinity * POW2(F2 / y_max) / m_Re0;
1293
1294 // compressible Poiseuille Temperature distribution
1295 // T = T0 - (m_gamma-F1)*m_Pr * POW2(PV->UInfinity*F3/F2) * (F1B2*(F1+ pow( (y-y_max/F2)/(y_max/F2), F4 )
1296 // )-POW2((y-y_max/F2)/(y_max/F2)));
1297 m_cells->pvariables[PV->RHO][cellid] = CV->rhoInfinity; // m_gamma * m_cells->variables[ PV->P ][cellid] /
1298 // T;//
1299 }
1300
1301 break;
1302 }
1303 case 101: // TAYLOR_GREEN_VORTEX
1304 {
1305 // domain boundaries are all 2*pi
1306 // rho=1.0;
1307 // u=A*SIN(x)*COS(y)*COS(z)
1308 // v=- A*COS(x)*SIN(y)*COS(z)
1309 // w=0.0
1310 // p=A*A*rho*( 1./(Ms*Ms*kappa) + 1./16.*(COS(2*x)*COS(2.*z)+ 2.*COS(2.*y) +2.*COS(2.*x) +COS(2*y)*COS(2.*z)))
1311 // Ms =0.1 maximum Mach number
1312 // A= speed magnitude set to 1.0
1315 MFloat x = F0;
1317 MFloat z = F0;
1330 PV->PInfinity
1332 }
1333 }
1334 }
1335 break;
1336 }
1337 case 1234: {
1338 // laminar channel flow
1340 m_channelHeight = 2.0;
1341 m_channelLength = 6.2831;
1342 // m_deltaP=32.0*SUTHERLANDLAW(PV->TInfinity)*m_Ma*sqrt(PV->TInfinity)*m_channelLength/(m_Re0*m_channelHeight);
1343 m_deltaP =
1344 -12.0 * PV->UInfinity * SUTHERLANDLAW(PV->TInfinity) * m_channelLength / (POW2(m_channelHeight) * m_Re0);
1345
1348 MFloat u =
1349 PV->UInfinity; // POW2(m_channelHeight)*m_deltaP*m_Re0*m_referenceLength*m_referenceLength/m_channelLength;
1354
1355 // channel height is in j direction
1356 // so prescribe a mean profile in u(y)
1358 pressureCH = m_channelPresInlet + (m_channelPresOutlet - m_channelPresInlet) * x / m_channelLength;
1359 // density:
1366 }
1367 }
1368 }
1369
1370 break;
1371 }
1372 case 1233: {
1373 // turbulent channel with perturbations
1374 // calculate the Pressure loss;
1375
1376 // for the law of the wall
1383 MFloat yplus = F0;
1386
1388
1390 // mean velocity profile
1392 m_channelPresOutlet =
1393 PV->PInfinity - CV->rhoInfinity * POW2(uTau) * F2 * (xINIT + m_channelLength) / m_channelHeight;
1397 m_log << "deltaP: " << deltaP << " cfTheoretisch: " << cfTheo << " cdTheo: " << cdTheo << endl;
1398
1399
1403 // channel height is in j direction
1404 // so prescribe a mean profile in u(y)
1406 // channel starts at x=0 or we need to prescribe an offset
1409 pressureCH =
1410 m_channelPresInlet + ((m_channelPresOutlet - m_channelPresInlet) * (x - xINIT) / m_channelLength);
1411
1413 MFloat velFactor = 2.0;
1415 yplus = prefactor * y;
1416 } else {
1418 }
1419
1420 // C1 etc are defined in maiaconstants.h
1421 if(yplus <= 5.0) {
1423 } else if(yplus <= 30 && yplus > 5.0) {
1425 } else if(yplus > 30) {
1427 }
1428
1433 }
1434 }
1435 }
1436 // create the fluctuations:
1437
1438
1439 // this way of creating the
1440 // fluctuations is only possible for
1441 // quite small mounts of cells
1442 // if too large this approach will not work anymore
1443 // 125000000 is a random integer number which has to be tested
1444 // else the approach with white noise will be employed
1447
1448 if(totalBlockCells <= 12500000) {
1449 fftw_complex *uPhysField, *vPhysField, *wPhysField;
1450
1451 // we need the total number of points
1452 MInt lx = m_grid->getMyBlockNoCells(2), ly = m_grid->getMyBlockNoCells(1), lz = m_grid->getMyBlockNoCells(0);
1453
1454 // field of velocities from positve frequencies
1455 uPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
1456 vPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
1457 wPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
1458
1459 // no real parallel computation is done
1460 // should be implemented later !!!!
1461
1467 MInt m_noPeakModes = 100;
1468 initFFTW(uPhysField, vPhysField, wPhysField, lx, ly, lz, m_noPeakModes);
1469 // copy values into the sendRCVbuffer
1470
1475 }
1476 }
1477 MPI_Bcast(&sendRecvBufferU[0], lx * ly * lz, MPI_DOUBLE, 0, m_StructuredComm, AT_, "sendRecvBufferU[0]");
1478 MPI_Bcast(&sendRecvBufferV[0], lx * ly * lz, MPI_DOUBLE, 0, m_StructuredComm, AT_, "sendRecvBufferV[0]");
1479 MPI_Bcast(&sendRecvBufferW[0], lx * ly * lz, MPI_DOUBLE, 0, m_StructuredComm, AT_, "sendRecvBufferW[0]");
1480
1486 }
1487 }
1488 } else {
1489 // create velocity field
1490 MInt m_noPeakModes = 100;
1491 initFFTW(uPhysField, vPhysField, wPhysField, lx, ly, lz, m_noPeakModes);
1492 }
1493 // now we need to distribute the
1494 } else {
1495 MFloat amp = 0.15;
1506 }
1507 }
1508 }
1509 }
1510
1511 break;
1512 }
1513 case 1236: {
1514 // pipe with perturbations
1515 // calculate the Pressure loss;
1516 // for the law of the wall
1517 // const MFloat C1 = m_channelC1;
1518 // const MFloat C2 = m_channelC2;
1523 // MFloat prefactor=m_ReTau;//uTau*m_Re0/SUTHERLANDLAW(PV->TInfinity);
1524
1529
1536 m_log << "=========== Turb. Pipe Flow Initial Condition Summary Finished =========== " << endl;
1537
1541 // IMPORTANT PIPE LENGTH HAS TO GO INTO X-DIRECTION
1542 // so prescribe a mean profile in u(r)
1543 // centerline is assumed at y=0.0,z=0.0
1544
1546 // determine the radius
1547 MFloat r = sqrt(POW2(m_cells->coordinates[1][cellId]) + POW2(m_cells->coordinates[2][cellId]));
1548
1550 // determine the pressure drop
1551 pressureCH = m_deltaP / m_channelLength * (x - m_channelInflowPlaneCoordinate) + PV->PInfinity;
1553 // important: viscous sublayer is not build explicitly. The flow has to build it by itself
1560 }
1561 }
1562 }
1563
1564 // create the fluctuations:
1565 // this way of creating the
1566 // fluctuations is only possible for
1567 // quite small mounts of cells
1568 // if too large this approach will not work anymore
1569 // 125000000 is a random integer number which has to be tested
1570 // else the approach with white noise will be employed
1573
1574 if(totalBlockCells <= 12500000) {
1575 fftw_complex *uPhysField, *vPhysField, *wPhysField;
1576
1577 // we need the total number of points
1578 MInt lx = m_grid->getMyBlockNoCells(2), ly = m_grid->getMyBlockNoCells(1), lz = m_grid->getMyBlockNoCells(0);
1579
1580 // field of velocities from positve frequencies
1581 uPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
1582 vPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
1583 wPhysField = (fftw_complex*)fftw_malloc(lx * ly * lz * sizeof(fftw_complex));
1584
1585 // no real parallel computation is done
1586 // should be implemented later !!!!
1587
1593 MInt m_noPeakModes = 100;
1594 initFFTW(uPhysField, vPhysField, wPhysField, lx, ly, lz, m_noPeakModes);
1595 // copy values into the sendRCVbuffer
1596
1601 }
1602 }
1603 MPI_Bcast(&sendRecvBufferU[0], lx * ly * lz, MPI_DOUBLE, 0, m_StructuredComm, AT_, "sendRecvBufferU[0]");
1604 MPI_Bcast(&sendRecvBufferV[0], lx * ly * lz, MPI_DOUBLE, 0, m_StructuredComm, AT_, "sendRecvBufferV[0]");
1605 MPI_Bcast(&sendRecvBufferW[0], lx * ly * lz, MPI_DOUBLE, 0, m_StructuredComm, AT_, "sendRecvBufferW[0]");
1606
1612 }
1613 }
1614 } else {
1615 // create velocity field
1616 MInt m_noPeakModes = 100;
1617 initFFTW(uPhysField, vPhysField, wPhysField, lx, ly, lz, m_noPeakModes);
1618 }
1619 // now we need to distribute the
1620 } else {
1621 MFloat amp = 0.15;
1632 }
1633 }
1634 }
1635 }
1636
1637 break;
1638 }
1639 case 79092: {
1640 // approximate mean turbulent boundary layer profile
1658
1661 } else if(yPlus <= 10.0) {
1663 } else if(yPlus <= 30 && yPlus > 10.0) {
1667 uTau
1668 * ((F1 / kappa) * log(max(yPlus, epss)) + C1 + 2 * (PI1 / kappa) * (3 * y * y - 2 * y * y * y));
1669 }
1670
1675
1678 }
1679 }
1680 }
1681 }
1682
1683 break;
1684 }
1685
1686 case 79091: // Turbulent plate
1687 {
1696
1703 const MFloat yplus = m_cells->coordinates[1][cellId] * sqrt(cf / 2.) * CV->rhoUInfinity / mu * m_Re0;
1705
1706 // 1-7th profile
1707 // log-law + wake
1710 } else if(yplus < 10) {
1712 } else {
1714 utau
1715 * ((1. / K) * log(max(yplus, epss)) + C1 + 2 * PI1 / K * (3 * eta * eta - 2 * eta * eta * eta)),
1716 PV->UInfinity);
1717 }
1718
1723
1726 }
1727 }
1728 }
1729 }
1730
1731 break;
1732 }
1733 case 11: // point source in the middle
1734 {
1735 // contain an initial perturbation in the middle
1737 // go through every cell
1741 }
1742
1743 MFloat amp = 0.0001;
1747 MFloat r = 0.025;
1753 }
1754 break;
1755 }
1756 /* TESTCASE from C. Bogey, C. Bailly, Three-dimensional non-reflective boundary conditions
1757 * for acoustic simulations: far field formulation and validation test cases
1758 * Acta acustica united with acustica Volume 88 (2002), 463-471
1759 */
1760 case 111: {
1761 MFloat amp = 0.01;
1762 MFloat alpha = log(2) / 9;
1764 // go through every cell
1774 }
1775 }
1776 break;
1777 }
1778 /* TESTCASE from C. Bogey, C. Bailly, Three-dimensional non-reflective boundary conditions
1779 * for acoustic simulations: far field formulation and validation test cases
1780 * Acta acustica united with acustica Volume 88 (2002), 463-471 the vortex
1781 */
1782 case 112: {
1783 MFloat amp = 0.003;
1787 MFloat r = F0;
1788 MFloat phi = F0;
1789 MFloat v = F0;
1791 // go through every cell
1801 phi = atan(z / y);
1807 }
1808 break;
1809 }
1810 case 113: {
1811 MFloat amp = 0.003;
1815 MFloat r = F0;
1816 MFloat phi = F0;
1817 MFloat v = F0;
1819 // go through every cell
1829 phi = atan(z / y);
1835 }
1836 break;
1837 }
1838 case 2: // shear flow is prescribed
1839 {
1841 MFloat x = F0;
1848 if(x < 0.5) {
1852 } else {
1856 }
1858 }
1859 }
1860 }
1861 break;
1862 }
1863 case 900: { // jet Freund
1876 }
1877 }
1878 }
1879 break;
1880 }
1881 case 4001: { // test periodic rotation boundary conditions
1886 // MFloat x= m_cells->coordinates[0][cellId];
1891 MFloat rmax = 10.0;
1897 }
1898 }
1899 }
1900 break;
1901 }
1902 case 42: {
1907
1910 } else {
1912 }
1913
1915 }
1916 break;
1917 }
1918 case 44: {
1919 MFloat amp = 0.003;
1923 MFloat r = F0;
1924 MFloat phi = F0;
1925 MFloat v = F0;
1927 // go through every cell
1937 phi = atan(z / y);
1943 }
1944 break;
1945 }
1946 case 777: {
1947 // fsc
1951 ofstream fscf;
1952 // write pressure to file
1953 fscf.open("pressure.dat", ios::trunc);
1954 if(fscf) {
1955 fscf << "#x_maia x_fsc p" << endl;
1959 }
1960 fscf.close();
1961 }
1962 // write velocity to file
1963 fscf.open("velocity_x0.dat", ios::trunc);
1964 if(fscf) {
1965 fscf << "#y eta u v w" << endl;
1967 // coord
1970 // velocity
1975 fscf << endl;
1976 }
1977
1978 fscf.close();
1979 }
1980 }
1981
1984 getFscVelocity(cellid, vel);
1987 }
1990 }
1991
1992 break;
1993 }
1994 case 999: {
1995 // blasius laminar bl
1997 if(!m_useBlasius) mTerm(1, "property Blasius not set. Refer to the description of the property");
1999 ofstream blasiusf;
2000 // write velocity to file
2001 blasiusf.open("velocity_x0.dat", ios::trunc);
2002 if(blasiusf) {
2003 blasiusf << "#y eta u v" << endl;
2004 MFloat d0 = 0.0;
2005 MFloat d1 = 0.0;
2006 MFloat d2 = 0.0;
2008
2010 // coord
2013 // velocity
2018 blasiusf << endl;
2019
2021 d0 = y;
2022 d0Set = true;
2023 }
2024
2026 d1 += (1 - vel[0] / PV->UInfinity) * 0.05;
2028 }
2029 }
2030 blasiusf.close();
2031
2033 cout << "d0: " << d0 << " d1: " << d1 << " d2: " << d2 << endl;
2034 }
2035 }
2036
2039 getBlasiusVelocity(cellid, vel);
2042 }
2045 }
2046
2047 break;
2048 }
2049 default: {
2050 // put the parallel flow field input in here
2051 // force output that no specific initial condition was chosen
2053 break;
2054 }
2055 }
2056 }
2057}
void initFFTW(fftw_complex *uPhysField, fftw_complex *vPhysField, fftw_complex *wPhysField, MInt lx, MInt ly, MInt lz, MInt noPeakModes)
Definition: fvstructuredsolver3d.cpp:9732
virtual MFloat getBlasiusEta(MFloat coordX, MFloat coordY)
Definition: fvstructuredsolver.cpp:7726
MInt m_initialCondition
Definition: fvstructuredsolver.h:660
MFloat m_referenceLength
Definition: fvstructuredsolver.h:430
MFloat m_channelInflowPlaneCoordinate
Definition: fvstructuredsolver.h:761
MInt m_blasius_noPoints
Definition: fvstructuredsolver.h:902
virtual void getBlasiusVelocity(MInt cellId, MFloat *const vel)
Load variables for the specified timeStep.
Definition: fvstructuredsolver.cpp:7721
MFloat m_channelPresOutlet
Definition: fvstructuredsolver.h:757
virtual MFloat getFscEta(MFloat coordX, MFloat coordY)
Definition: fvstructuredsolver.cpp:7586
virtual MFloat getFscPressure(MInt cellId)
Definition: fvstructuredsolver.cpp:7560
MFloat m_channelPresInlet
Definition: fvstructuredsolver.h:756
virtual void getFscVelocity(MInt cellId, MFloat *const vel)
Load variables for the specified timeStep.
Definition: fvstructuredsolver.cpp:7581
int MPI_Bcast(void *buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm, const MString &name, const MString &varname)
same as MPI_Bcast
Definition: mpioverride.cpp:1114
◆ initializeNeighbourCellIds()
| void FvStructuredSolver3D::initializeNeighbourCellIds | ( | ) |
◆ initInterpolatedPoints()
| void FvStructuredSolver3D::initInterpolatedPoints | ( | ) |
Definition at line 2063 of file fvstructuredsolver3d.cpp.
2063 {
2064 TRACE();
2066 m_pointInterpolation =
2067 make_unique<StructuredInterpolation<3>>(m_nCells, m_cells->coordinates, m_cells->pvariables, m_StructuredComm);
2068 // allocate domains points of the lines
2070 m_intpPointsHasPartnerLocal);
2071 MPI_Allreduce(m_intpPointsHasPartnerLocal, m_intpPointsHasPartnerGlobal, m_intpPointsNoPointsTotal, MPI_INT,
2072 MPI_SUM, m_StructuredComm, AT_, "m_intpPointsHasPartnerLocal", "m_intpPointsHasPartnerGlobal");
2073 }
2074}
std::unique_ptr< StructuredInterpolation< 3 > > m_pointInterpolation
Definition: fvstructuredsolver3d.h:208
MInt m_intpPointsNoPointsTotal
Definition: fvstructuredsolver.h:496
MInt * m_intpPointsHasPartnerGlobal
Definition: fvstructuredsolver.h:504
MInt * m_intpPointsHasPartnerLocal
Definition: fvstructuredsolver.h:505
MFloat ** m_intpPointsCoordinates
Definition: fvstructuredsolver.h:503
◆ initMovingGrid()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 6492 of file fvstructuredsolver3d.cpp.
◆ initPointsToAsciiFile()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 2077 of file fvstructuredsolver3d.cpp.
2077 {
2078 TRACE();
2080 m_pointsToAsciiInterpolation =
2081 make_unique<StructuredInterpolation<3>>(m_nCells, m_cells->coordinates, m_cells->pvariables, m_StructuredComm);
2082
2083 // allocate domains points of the lines
2084 m_pointsToAsciiInterpolation->prepareInterpolation(m_pointsToAsciiNoPoints, m_pointsToAsciiCoordinates,
2086 MPI_Allreduce(m_pointsToAsciiHasPartnerLocal, m_pointsToAsciiHasPartnerGlobal, m_pointsToAsciiNoPoints, MPI_INT,
2087 MPI_SUM, m_StructuredComm, AT_, "m_pointsToAsciiHasPartnerLocal", "m_pointsToAsciiHasPartnerGlobal");
2088
2093 << " x: " << m_pointsToAsciiCoordinates[0][pointId] << " y: " << m_pointsToAsciiCoordinates[1][pointId]
2095 }
2096 }
2097}
MFloat ** m_pointsToAsciiCoordinates
Definition: fvstructuredsolver.h:541
MInt * m_pointsToAsciiHasPartnerGlobal
Definition: fvstructuredsolver.h:545
MInt m_pointsToAsciiNoPoints
Definition: fvstructuredsolver.h:538
MInt * m_pointsToAsciiHasPartnerLocal
Definition: fvstructuredsolver.h:546
std::unique_ptr< StructuredInterpolation< nDim > > m_pointsToAsciiInterpolation
Definition: fvstructuredsolver.h:547
◆ initSandpaperTrip()
| void FvStructuredSolver3D::initSandpaperTrip | ( | ) |
Definition at line 3741 of file fvstructuredsolver3d.cpp.
◆ initSolutionStep()
|
virtual |
Implements FvStructuredSolver< 3 >.
Definition at line 978 of file fvstructuredsolver3d.cpp.
978 {
979 TRACE();
980
981 std::ignore = mode;
982
983 // Compute infinity values from property file
984 // and (if no restart) fill cells according
985 // to the initialCondition property
986 initialCondition();
987
989 interpolateFromDonor();
991 loadRestartFile();
992 }
993
994 // timestep will be computed and
995 // set if globalTimeStep == 0 or
996 // constantTimeStep = 0
997 setTimeStep();
998
999 // initialize moving grid
1000 // functions and move grid
1001 // to correct position
1005 initMovingGrid();
1008 }
1009
1011 initBodyForce();
1012 }
1013
1015 initSandpaperTrip();
1016 }
1017
1018 // Get the correct values
1019 // in the exchange ghostcells
1020 exchange();
1021
1023 m_ransSolver->computeTurbViscosity();
1024 }
1025
1029 zonalExchange();
1030 }
1031
1032 // Call the init function of each BC
1033 initBndryCnds();
1034
1035 // Apply boundary conditions
1036 // and fill the non-exchange ghostcells
1037 applyBoundaryCondition();
1038
1040 m_ransSolver->computeTurbViscosity();
1041 }
1042
1043 // Check for NaNs
1044 checkNans();
1045
1046 computeConservativeVariables();
1047}
void initSandpaperTrip()
Definition: fvstructuredsolver3d.cpp:3724
void initBodyForce() override
Definition: fvstructuredsolver3d.cpp:8020
void zonalExchange() override
Definition: fvstructuredsolver3d.cpp:546
void spanwiseAvgZonal(std::vector< MFloat * > &) override
Definition: fvstructuredsolver3d.cpp:4303
void applyBoundaryCondition() override
Definition: fvstructuredsolver3d.cpp:2274
void interpolateFromDonor()
Interpolates the flow field from a given donor file.
Definition: fvstructuredsolver3d.cpp:1058
void initialCondition() override
Computation of infinity values for the conservative and primitive variables.
Definition: fvstructuredsolver3d.cpp:1112
void computeCumulativeAverage(MBool forceReset) override
Definition: fvstructuredsolver3d.cpp:4267
void initMovingGrid() override
Definition: fvstructuredsolver3d.cpp:6449
void checkNans()
Checks whole domain for NaNs and adds the number of NaNs globally.
Definition: fvstructuredsolver.cpp:7029
void loadRestartFile()
Load Restart File (primitive and conservative output) general formulation.
Definition: fvstructuredsolver.cpp:5630
virtual void computeConservativeVariables()
Definition: fvstructuredsolver.cpp:4456
MBool m_useSandpaperTrip
Definition: fvstructuredsolver.h:839
MBool m_restartInterpolation
Definition: fvstructuredsolver.h:662
◆ interpolateFromDonor()
|
protected |
Instead of a restart from a restart file this method will interpolate from a donor solution onto the flow field
Definition at line 1058 of file fvstructuredsolver3d.cpp.
1058 {
1061
1065 }
1066
1068 } else {
1071 donorConservative = Context::getSolverProperty<MBool>("donorConservative", false, AT_);
1072 }
1073
1074 if(donorConservative) {
1077 }
1078 computePrimitiveVariables();
1079 } else {
1082 }
1083 }
1084
1085 // interpolate NU_T for STG if this is an initial start (otherwise not needed)
1088 }
1089 }
1090
1094 }
1095
1098 }
1099
1100 manualInterpolationCorrection();
1101}
static MBool propertyExists(const MString &name, MInt solver=m_noSolvers)
This function checks if a property exists in general.
Definition: context.cpp:494
void computePrimitiveVariables() override
Definition: fvstructuredsolver3d.cpp:9703
void manualInterpolationCorrection()
Manually correct errors made by the restart interpolation.
Definition: fvstructuredsolver3d.cpp:2189
std::unique_ptr< StructuredInterpolation< 3 > > m_structuredInterpolation
Definition: fvstructuredsolver3d.h:207
MBool m_stgInitialStartup
Definition: fvstructuredsolver.h:593
MInt m_spongeLayerType
Definition: fvstructuredsolver.h:621
MString * m_variableNames
Definition: fvstructuredsolver.h:410
MString * m_pvariableNames
Definition: fvstructuredsolver.h:411
◆ loadAveragedVariables()
|
overrideprotectedvirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 9546 of file fvstructuredsolver3d.cpp.
9546 {
9547 TRACE();
9549
9550 // open the file
9551 MString fileNameStr = fileName;
9553
9554 stringstream blockNumber;
9555 blockNumber << m_blockId;
9557 blockPathStr += blockNumber.str();
9560
9562 pio.readArray(m_summedVars[var], blockPathStr, m_avgVariableNames[var], nDim, ioOffset, ioSize);
9563 }
9564
9568 }
9569 }
9570
9572 shiftAverageCellValues();
9573}
void shiftAverageCellValues()
Shifts the averaged variables.
Definition: fvstructuredsolver3d.cpp:9636
Definition: parallelio_hdf5.h:19
MFloat ** m_summedVars
Definition: fvstructuredpostprocessing.h:93
MBool m_averagingFavre
Definition: fvstructuredpostprocessing.h:119
MString * m_avgVariableNames
Definition: fvstructuredpostprocessing.h:121
MInt getNoVars()
Definition: fvstructuredpostprocessing.cpp:2424
MString * m_avgFavreNames
Definition: fvstructuredpostprocessing.h:122
MInt getNoPPVars()
Returns number of postprocessing variables.
Definition: fvstructuredpostprocessing.cpp:2381
MFloat ** m_favre
Definition: fvstructuredpostprocessing.h:98
◆ loadAverageRestartFile()
|
overrideprotectedvirtual |
- Date
- 06.01.2016
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 9464 of file fvstructuredsolver3d.cpp.
9465 {
9466 TRACE();
9469 cout << "Opening file: " << fileName << endl;
9470 }
9471
9472 // open the file
9473 MString fileNameStr = fileName;
9475
9476 // now read in the data!
9480 stringstream blockNumber;
9481 blockNumber << m_blockId;
9483 blockPathStr += blockNumber.str();
9484 MInt offset = 0;
9487
9488
9494 offset = noVariables();
9495
9502 }
9503
9508 }
9509
9516
9518
9523 }
9524
9529 }
9530
9535 }
9536
9538 shiftAverageCellValuesRestart();
9539}
void shiftAverageCellValuesRestart()
Definition: fvstructuredsolver3d.cpp:9580
MInt noVariables() const override
Definition: fvstructuredsolver.h:211
MInt m_noSamples
Definition: fvstructuredpostprocessing.h:143
MBool m_kurtosis
Definition: fvstructuredpostprocessing.h:117
MBool m_skewness
Definition: fvstructuredpostprocessing.h:116
MBool m_averageVorticity
Definition: fvstructuredpostprocessing.h:118
◆ loadRestartBC2600()
|
virtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 8617 of file fvstructuredsolver3d.cpp.
8617 {
8620 cout << "Loading BC2600 values..." << endl;
8621 }
8623 ParallelIo::size_type bcCells[3] = {m_grid->getMyBlockNoCells(0), m_grid->getMyBlockNoCells(1), m_noGhostLayers};
8624 MInt noCellsBC = bcCells[0] * bcCells[1] * bcCells[2];
8625 ParallelIo::size_type bcOffset[3] = {0, 0, 0};
8628 stringstream restartFileName;
8629 MString restartFile = Context::getSolverProperty<MString>("restartVariablesFileName", m_solverId, AT_);
8631
8633 stringstream pathStr;
8635 std::string path = pathStr.str();
8636
8638 pio.readArray(&tmpRestartVars[var * noCellsBC], path, m_pvariableNames[var], nDim, bcOffset, bcCells);
8639 }
8640 }
8641
8642 MPI_Bcast(&tmpRestartVars(0, 0), noCellsBC * m_maxNoVariables, MPI_DOUBLE, 0, m_commBC2600, AT_,
8643 "tmpRestartVars(0");
8644
8645 MInt startGC[3] = {0, 0, 0};
8646 MInt endGC[3] = {0, 0, 0};
8647
8649 startGC[1] = m_noGhostLayers;
8650 }
8652 startGC[0] = m_noGhostLayers;
8653 }
8655 endGC[1] = m_noGhostLayers;
8656 }
8658 endGC[0] = m_noGhostLayers;
8659 }
8660
8665 MInt globalI = i;
8668 MInt cellIdBC = globalI + (globalJ + globalK * bcCells[1]) * bcCells[2];
8669
8670 // load values from restart field
8673 }
8674 }
8675 }
8676 }
8677
8678
8684
8688 }
8689
8693
8697 }
8698 }
8699 }
8700
8701
8702 // Fix diagonal cells at end of domain
8707 i + (m_noGhostLayers + m_bc2600noActiveCells[1] - 2 + k * m_bc2600noCells[1]) * m_bc2600noCells[2];
8709 i + (m_noGhostLayers + m_bc2600noActiveCells[1] - 1 + k * m_bc2600noCells[1]) * m_bc2600noCells[2];
8711 i + (m_noGhostLayers + m_bc2600noActiveCells[1] + k * m_bc2600noCells[1]) * m_bc2600noCells[2];
8717
8718 const MFloat slope = (m_bc2600Variables[var][cellIdA1] - m_bc2600Variables[var][cellIdA2]) / distA1A2;
8724 }
8725 }
8726 }
8727 }
8728 }
8729
8731 cout << "Loading BC2600 values... SUCCESSFUL!" << endl;
8732 }
8733 }
8734}
MInt * m_bc2600noOffsetCells
Definition: fvstructuredsolver.h:783
MInt * m_bc2600noActiveCells
Definition: fvstructuredsolver.h:782
MFloat ** m_bc2600Variables
Definition: fvstructuredsolver.h:779
MInt m_commBC2600MyRank
Definition: fvstructuredsolver.h:708
MBool m_bc2600InitialStartup
Definition: fvstructuredsolver.h:778
◆ loadRestartBC2601()
|
virtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 8737 of file fvstructuredsolver3d.cpp.
8737 {
8738 // if we have the prescribing boundary,
8739 // put the values from the restart file into the ghostcells
8742 cout << "Loading restart values 2601" << endl;
8743 }
8744 ParallelIo::size_type bcCells[3] = {0, 0, 0};
8745 ParallelIo::size_type bcOffset[3] = {0, 0, 0};
8746
8747 bcCells[0] = m_grid->getMyBlockNoCells(0);
8748 bcCells[1] = m_noGhostLayers;
8749 bcCells[2] = m_grid->getMyBlockNoCells(2);
8750
8751 MInt noCellsBC = bcCells[0] * bcCells[1] * bcCells[2];
8753
8755 stringstream restartFileName;
8756 MString restartFile = Context::getSolverProperty<MString>("restartVariablesFileName", m_solverId, AT_);
8758
8760 stringstream pathStr;
8762 std::string path = pathStr.str();
8763
8765 pio.readArray(&tmpRestartVars[var * noCellsBC], path, m_pvariableNames[var].c_str(), nDim, bcOffset, bcCells);
8766 }
8767 }
8768
8769 MPI_Bcast(&tmpRestartVars[0], noCellsBC * PV->noVariables, MPI_DOUBLE, 0, m_StructuredComm, AT_,
8770 "tmpRestartVars[0]");
8771
8773 MInt startGC[3] = {0, 0, 0};
8774 MInt endGC[3] = {0, 0, 0};
8775
8777 startGC[2] = m_noGhostLayers;
8778 }
8780 startGC[0] = m_noGhostLayers;
8781 }
8783 endGC[2] = m_noGhostLayers;
8784 }
8786 endGC[0] = m_noGhostLayers;
8787 }
8788
8794 MInt globalJ = j;
8796 MInt cellIdBC = globalI + (globalJ + globalK * m_noGhostLayers) * m_grid->getMyBlockNoCells(2);
8797
8798 // load values from restart field
8801 }
8802 }
8803 }
8804 }
8805
8806 // Fix diagonal cells at start of domain
8819 }
8820 }
8821 }
8822 }
8823
8824 // Fix diagonal cells at end of domain
8837 }
8838 }
8839 }
8840 }
8841 }
8842 }
8843}
MBool m_bc2601InitialStartup
Definition: fvstructuredsolver.h:789
◆ loadRestartSTG()
|
virtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 8845 of file fvstructuredsolver3d.cpp.
8845 {
8848 cout << "Restart file has conservative variables, converting STG variables to primitive!" << endl;
8849 }
8850
8851 stringstream restartFileName;
8852 MString restartFile = Context::getSolverProperty<MString>("restartVariablesFileName", m_solverId, AT_);
8855
8856 stringstream blockNumber;
8857 blockNumber << m_blockId;
8859 blockPathStr += blockNumber.str();
8860
8861 MInt restartNoEddies = 0;
8862 pio.getAttribute(&restartNoEddies, "stgNRAN", "");
8863
8865 m_log << "STG: NRAN in restart file (" << restartNoEddies << ") not the same as given in property file ("
8868 } else {
8870 }
8871
8872 // if this is an initialStartup the new
8873 // eddies need to be created in any case
8878 // also load nu_t into fq field
8879 pio.readArray(m_cells->fq[FQ->NU_T], blockPathStr, FQ->fqNames[FQ->NU_T].c_str(), nDim, offset, size);
8881 } else {
8882 // has to be set manually for restart from RANS profile
8884 ParallelIo::size_type VBStart = 0;
8885
8886
8892 cout << "NRAN in property differs from NRAN in restart file"
8893 << " NOT READING EDDIES FROM RESTART!" << endl;
8894 }
8895 } else {
8897
8898 if(pio.hasDataset("FQeddies", stgGlobalPathStr)) {
8900 cout << "FQeddies field is at new position /stgGlobal/FQeddies" << endl;
8901 }
8902 } else {
8904 cout << "FQeddies field is NOT at new position! Using old path within solver..." << endl;
8905 }
8906 stgGlobalPathStr = blockPathStr;
8908 cout << "Loading FQeddies from path " << stgGlobalPathStr << endl;
8909 }
8910 }
8911
8912 // do this only serially
8914 cout << "Loading STG Eddies..." << endl;
8915 }
8919 }
8920
8923 cout << "Loading STG Eddies... SUCCESSFUL!" << endl;
8924 }
8925 }
8926
8927
8931
8933 ParallelIo::size_type bcActiveCells[3] = {m_grid->getMyBlockNoCells(0), m_grid->getMyBlockNoCells(1), 3};
8936 ParallelIo::size_type bcOffset[3] = {0, 0, 0};
8939
8941 cout << "Loading STG Datasets..." << endl;
8942 }
8947 stringstream fieldName;
8948 stringstream stgPath;
8949 stgPath << blockPathStr << "/stg";
8950 fieldName << "stgFQ" << var;
8951 pioLocal.readArray(&tmpRestartVars[var * noCellsBC], stgPath.str(), fieldName.str(), nDim, bcOffset,
8952 bcActiveCells);
8953 }
8954
8955 // reorder the cells in the global array
8964
8966 tmpRestartVars[var * noCellsBC + cellId] = tmpRestartVars[var * noCellsBC + cellId_org];
8967 tmpRestartVars[var * noCellsBC + cellId_org] = F0;
8968 }
8969 }
8970 }
8971 }
8972
8973 // apply periodic bc
8980
8985
8987 tmpRestartVars[var * noCellsBC + gcId0] = tmpRestartVars[var * noCellsBC + acId0];
8988 tmpRestartVars[var * noCellsBC + gcId1] = tmpRestartVars[var * noCellsBC + acId1];
8989 tmpRestartVars[var * noCellsBC + gcId2] = tmpRestartVars[var * noCellsBC + acId2];
8990 tmpRestartVars[var * noCellsBC + gcId3] = tmpRestartVars[var * noCellsBC + acId3];
8991 }
8992 }
8993 }
8994
8995 // extrapolation at the bottoms
9002
9004 tmpRestartVars[var * noCellsBC + gc1] =
9005 2 * tmpRestartVars[var * noCellsBC + ac1] - tmpRestartVars[var * noCellsBC + ac2];
9006 tmpRestartVars[var * noCellsBC + gc2] =
9007 2 * tmpRestartVars[var * noCellsBC + gc1] - tmpRestartVars[var * noCellsBC + ac1];
9008 }
9009 }
9010 }
9011
9012 // extrapolation at the top
9019
9020
9022 tmpRestartVars[var * noCellsBC + gc1] =
9023 2 * tmpRestartVars[var * noCellsBC + ac1] - tmpRestartVars[var * noCellsBC + ac2];
9024 tmpRestartVars[var * noCellsBC + gc2] =
9025 2 * tmpRestartVars[var * noCellsBC + gc1] - tmpRestartVars[var * noCellsBC + ac1];
9026 }
9027 }
9028 }
9029 }
9030
9031
9032 // now broadcast to everyone
9033 MPI_Bcast(&tmpRestartVars[0], noCellsBC * m_stgNoVariables, MPI_DOUBLE, 0, m_commStg, AT_, "tmpRestartVars[0]");
9035 cout << "Loading STG Datasets... SUCCESSFUL!" << endl;
9036 }
9037
9045 MInt globalI = i;
9048 MInt cellIdBCGlobal = globalI + (globalJ + globalK * bcCells[1]) * bcCells[2];
9050
9051 // load values from restart field
9054 }
9055
9056 if(!isPrimitiveOutput) {
9062
9068
9074 }
9075
9076 if(i < 2) {
9082 }
9083 }
9084 }
9085 }
9086 }
9087 }
9088
9089 }
9090}
MInt m_stgNoEddieProperties
Definition: fvstructuredsolver.h:594
MBool m_stgCreateNewEddies
Definition: fvstructuredsolver.h:602
◆ loadSampleFile()
|
overrideprotectedvirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 9393 of file fvstructuredsolver3d.cpp.
9393 { // loading files for averaging (pre- and postsolve)
9394 TRACE();
9395
9396 // open the file
9398
9399 stringstream blockNumber;
9400 blockNumber << m_blockId;
9402 blockPathStr += blockNumber.str();
9405
9407 pio.readArray(m_cells->pvariables[var], blockPathStr, m_pvariableNames[var], nDim, offset, size);
9408 }
9409
9411}
void shiftCellValuesRestart(MBool)
◆ manualInterpolationCorrection()
|
protected |
In case some cells did not get correct values in the interpolation process (due to no matching donor cells) Put your routines that will correct these cells here.
Definition at line 2190 of file fvstructuredsolver3d.cpp.
◆ maxResidual()
|
virtual |
This function computes the maxResidual using Res = deltaT/(CFL*VolOfCell) * |RHS|
with deltaT depending on local or global time stepping is used. checks if the computed max density residual is below the convergence criterion and returns boolean variable
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 5554 of file fvstructuredsolver3d.cpp.
5554 {
5555 TRACE();
5556
5558 MFloat epsilon = pow(10.0, -10.0);
5559 m_avrgResidual = F0;
5561 MFloat tmpResidual = F0;
5562 MFloat maxResidual1 = F0;
5563 MInt maxResIndex[3];
5564 // MInt localCounter=F0;
5565 MFloat maxResidualOrg = F0;
5566 MFloat localMaxResidual = F0;
5567 MFloat localAvrgResidual = F0;
5568 MFloat accumAvrgResidual = F0;
5569 MFloat globalMaxResidual = F0;
5570 m_workload = F0;
5571 // MInt accumCounter=0;
5573 maxResIndex[dim] = F0;
5574 }
5575
5580 // cerr << cellId << endl;
5581 tmpResidual = m_timeStep / (m_cfl * m_cells->cellJac[cellId]) * fabs(m_cells->rightHandSide[CV->RHO][cellId]);
5582 m_avrgResidual += tmpResidual;
5583
5584 if(tmpResidual > maxResidual1) {
5585 maxResIndex[0] = i - m_noGhostLayers;
5586 maxResIndex[1] = j - m_noGhostLayers;
5587 maxResIndex[2] = k - m_noGhostLayers;
5588 maxResidual1 = tmpResidual;
5589 }
5590 }
5591 }
5592 }
5593
5594 // localCounter = counter;
5595 localMaxResidual = maxResidual1;
5596 localAvrgResidual = m_avrgResidual;
5597 // reset average Residual
5598 m_avrgResidual = F0;
5599
5600 MFloat localTotalEnergy = F0;
5601 MFloat globalTotalEnergy = F0;
5602 MFloat globalPressure = F0;
5603 MFloat localPressure = F0;
5605 localTotalEnergy = computeTotalKineticEngergy();
5606 MPI_Allreduce(&localTotalEnergy, &globalTotalEnergy, 1, MPI_DOUBLE, MPI_SUM, m_StructuredComm, AT_,
5607 "localTotalEnergy", "globalTotalEnergy");
5608 globalTotalEnergy /= ((16.0 * pow(4 * atan(1), 3.0))); // divided by the overall volume
5609 localPressure = computeTotalPressure();
5610 MPI_Allreduce(&localPressure, &globalPressure, 1, MPI_DOUBLE, MPI_SUM, m_StructuredComm, AT_, "localPressure",
5611 "globalPressure");
5612 globalPressure /= ((16.0 * pow(4 * atan(1), 3.0))); // divided by the overall volume
5613 }
5614
5615 // if( noDomains()>1 )
5616 // {
5617 // MPI_Allreduce(m_residualSnd, &m_residualRcv, 1, m_mpiStruct, m_resOp, m_StructuredComm, AT_, "m_residualSnd",
5618 // "m_residualRcv" );
5619 MPI_Allreduce(&localAvrgResidual, &accumAvrgResidual, 1, MPI_DOUBLE, MPI_SUM, m_StructuredComm, AT_,
5620 "localAvrgResidual", "accumAvrgResidual");
5621 MPI_Allreduce(&localMaxResidual, &globalMaxResidual, 1, MPI_DOUBLE, MPI_MAX, m_StructuredComm, AT_,
5622 "localMaxResidual", "globalMaxResidual");
5623
5625 maxResidualOrg = globalMaxResidual;
5626 // globalMaxResidual=globalMaxResidual;//m_residualRcv.maxRes;
5627 // for(MInt i=0; i<3; i++)
5628 //{
5629 // maxResIndex[i]=m_residualRcv.maxCellIndex[i];
5630 //}
5631
5632 // }
5633
5634 // cout << "m_avrgResidual = " << m_avrgResidual<< " | totalCells " << m_totalGridCells <<endl;
5636 // write first residuals;
5642 m_eps)) { // so only domainId with the max writes out ==> no need to communicate the max index[i]
5643
5644 // write out values into residual file
5645 FILE* f_residual;
5646 f_residual = fopen("./Residual", "a+");
5649 fprintf(f_residual, "#iter, physTime, time, dT, wLoad, avrgRes, maxRes, blockId, i, j, k ");
5650 fclose(f_residual);
5651 }
5652 } else {
5654 // write out values into residual file
5655 FILE* f_residual;
5656 f_residual = fopen("./Residual", "a+");
5659 fprintf(f_residual,
5660 "#iter, physTime, time, dT, wLoad, avrgRes, maxRes, blockId, i, j, k, k_mean, p_mean, pProbe ");
5661 }
5662 }
5663 }
5664
5665 // normalize residuals
5666 globalMaxResidual = globalMaxResidual / m_firstMaxResidual;
5668
5669 // question if "( m_avrgResidual >= F0 || m_avrgResidual < F0 ) {} else {"
5670 // is better to capture the also inf???
5671
5673 cerr << "Solution diverged, average residual is nan " << endl;
5677 }
5678
5679 // convergence Check
5680
5684 }
5685
5686 // need again special treatment for TGV due to symmetrie many processors would write out (!!!!IMPORTANT NO CORRECT
5687 // INDEX WILL BE WRITTEN OUT )
5689 // processor with the highest Residual writes out!!! saves communication;
5691 m_eps)) { // so only domainId with the max writes out ==> no need to communicate the max index[i]
5692 // write out values into residual file
5693 FILE* f_residual;
5694 f_residual = fopen("./Residual", "a+");
5701 fprintf(f_residual, " %1.10e", globalMaxResidual);
5706 fprintf(f_residual, "\n");
5707 fclose(f_residual);
5708 }
5709 } else {
5711 MFloat dissip = F0;
5713 m_kineticEOld = 1.2474901617e-03;
5714 } else {
5716 m_kineticEOld = globalTotalEnergy;
5717 }
5718 // write out values into residual file
5719 // compute the dissipation rate
5720
5721 FILE* f_residual;
5722 f_residual = fopen("./Residual", "a+");
5729 fprintf(f_residual, " %1.10e", globalMaxResidual);
5734 fprintf(f_residual, " %1.10e", globalTotalEnergy); // kinetic Energy
5735 fprintf(f_residual, " %1.10e", globalPressure); // averaged pressure
5736 fprintf(f_residual, " %1.10e", dissip); // dissipation rate
5737 fprintf(f_residual, "\n");
5738 fclose(f_residual);
5739 }
5740 }
5741
5743 return true;
5744 } else {
5745 return false;
5746 }
5747}
MFloat computeTotalPressure()
Definition: fvstructuredsolver3d.cpp:5765
MFloat computeTotalKineticEngergy()
Definition: fvstructuredsolver3d.cpp:5749
MFloat m_firstAvrgResidual
Definition: fvstructuredsolver.h:647
MFloat m_convergenceCriterion
Definition: fvstructuredsolver.h:656
void saveSolverSolution(MBool=false, const MBool=false) override
Definition: fvstructuredsolver.cpp:5073
MFloat m_firstMaxResidual
Definition: fvstructuredsolver.h:646
MInt m_residualInterval
The number of timesteps before writing the next residual.
Definition: solver.h:87
◆ moveGrid()
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 5798 of file fvstructuredsolver3d.cpp.
5798 {
5799 TRACE();
5802
5804 case 1: // oscillating grid in x-direction
5805 {
5806 // we need some relaxation function in wall-normal direction
5807 // here we take a linear function: ( y_max - y ) / y_max
5811
5817
5818 m_grid->m_coordinates[0][pointId] =
5819 x + (1 - m_grid->m_coordinates[1][pointId] / y_max) * amp * sin(2 * pi * t * frequency);
5820 m_grid->m_velocity[0][pointId] = (1 - m_grid->m_coordinates[1][pointId] / y_max) * amp * 2 * pi * frequency
5821 * cos(2 * pi * t * frequency);
5822 }
5823 }
5824 }
5825
5826 break;
5827 }
5828 case 2: // channel with moving indentation
5829 {
5833 const MFloat y_max = F1; // can be used to scale (old mesh had y_max = m_mgOrgCoordinates[0][32] = 1/30)
5834
5835 // values from Ralph and Pedley:
5843
5849
5850 MFloat g = F0;
5851
5853 g = (1.0 + tanhl(beta * (m_grid->m_coordinates[0][pointId] - (x2 + x3) / 2.0l) / y_max)) / 2.0;
5855 g = F1;
5857 g = (1.0 - tanhl(beta * (m_grid->m_coordinates[0][pointId] - (x4 + x5) / 2.0l) / y_max)) / 2.0;
5858 }
5859
5862 }
5863 }
5864 }
5865
5866 break;
5867 }
5868 case 3: // piston moving in x-direction
5869 {
5875
5878 m_grid->m_acceleration[0][pointId] = F0;
5879 }
5880 }
5881 }
5882
5883 break;
5884 }
5885 case 4: // inner grid movement
5886 {
5891
5892 // for Square:
5897
5901
5906
5909
5910 MFloat g = F0;
5914 }
5915 m_grid->m_coordinates[0][pointId] = x * (F1 - h * g * (F1 - x));
5916 m_grid->m_velocity[0][pointId] = x * (-hvel * g * (F1 - x));
5917
5918 g = F0;
5920 g = ((x < x3) ? (1.0l + tanhl(beta * (x - (x2 + x3) / 2.0l) / y_max)) / 2.0l
5921 : ((x < x4) ? 1.0l : (1.0l - tanhl(beta * (x - (x4 + x5) / 2.0l) / y_max)) / 2.0l));
5922 }
5925 }
5926 }
5927 }
5928
5929 break;
5930 }
5931 case 9:
5932 case 10: {
5933 // traveling wave case
5937 if(zeroPos) {
5938 t_offset = F0;
5939 }
5943
5944 MFloat fadeInFactor = 0.0;
5945 MFloat fadeInFactorPrime = 0.0;
5946 MFloat fadeInFactorPrimePrime = 0.0;
5949 fadeInFactorPrime = (pi / m_waveTemporalTransition) * sin(t_offset / m_waveTemporalTransition * pi) * F1B2;
5950 fadeInFactorPrimePrime =
5952 } else {
5953 fadeInFactor = 1.0;
5954 fadeInFactorPrime = 0.0;
5955 fadeInFactorPrimePrime = 0.0;
5956 }
5957
5958 if(zeroPos) {
5959 fadeInFactor = F1;
5960 }
5961
5962
5970
5971 MFloat transitionFactor = F0;
5973 transitionFactor = F0;
5977 transitionFactor = F1;
5979 transitionFactor = (1 + cos((xInit - m_waveOutBeginTransition) / transitionOutLength * pi)) * F1B2;
5980 } else {
5981 transitionFactor = F0;
5982 }
5983
5984 MFloat yTransitionFactor = F1;
5986 yTransitionFactor = F1;
5988 yTransitionFactor = (1 + cos((yInit - m_waveYBeginTransition) / yTransitionLength * pi)) * F1B2;
5989 } else {
5990 yTransitionFactor = F0;
5991 }
5992
5994
5997 const MFloat funcPrime = transitionFactor * yTransitionFactor * (2 * PI * m_waveSpeed / m_waveLength)
5998 * m_waveAmplitude
6003
6004 m_grid->m_coordinates[1][pointId] = func * fadeInFactor + yInit;
6005 m_grid->m_velocity[1][pointId] = func * fadeInFactorPrime + funcPrime * fadeInFactor;
6006 m_grid->m_acceleration[1][pointId] =
6007 funcPrimePrime * fadeInFactor + 2 * funcPrime * fadeInFactorPrime + func * fadeInFactorPrimePrime;
6008 }
6009 }
6010 }
6011
6012 break;
6013 }
6014 case 11: {
6015 // traveling wave channel (Tomiyama & Fukagata 2013)
6017 if(zeroPos) {
6018 t_offset = F0;
6019 }
6020
6021 MFloat fadeInFactor = 0;
6023
6024 if(t_offset < timeRelaxation) {
6025 fadeInFactor = (1.0 - cos(t_offset / timeRelaxation * pi)) * F1B2;
6026 } else {
6027 fadeInFactor = 1.0;
6028 }
6029
6030 if(zeroPos) {
6031 fadeInFactor = F1;
6032 }
6033
6040 MFloat yRelaxation = F0;
6041
6042 if(yInit <= F0) {
6043 yRelaxation = F1;
6044 } else if(yInit > F0 && yInit < F1) {
6045 yRelaxation = F1 - yInit;
6046 } else if(yInit > F1 && yInit < F2) {
6047 yRelaxation = yInit - F1;
6048 } else {
6049 yRelaxation = F1;
6050 }
6051
6052 if(yInit <= F1) {
6053 m_grid->m_coordinates[1][pointId] =
6054 yInit
6055 + fadeInFactor
6056 * (m_waveAmplitude * yRelaxation
6058 m_grid->m_velocity[1][pointId] = fadeInFactor
6061 } else {
6062 m_grid->m_coordinates[1][pointId] =
6063 yInit
6064 - fadeInFactor
6065 * (m_waveAmplitude * yRelaxation
6067 m_grid->m_velocity[1][pointId] = -fadeInFactor
6070 }
6071 }
6072 }
6073 }
6074
6075 break;
6076 }
6077 case 12: {
6078 // streamwise traveling wave case
6080 if(zeroPos) {
6081 t_offset = F0;
6082 }
6086
6087 MFloat fadeInFactor = 0.0;
6088 MFloat fadeInFactorPrime = 0.0;
6089 MFloat fadeInFactorPrimePrime = 0.0;
6092 fadeInFactorPrime = (pi / m_waveTemporalTransition) * sin(t_offset / m_waveTemporalTransition * pi) * F1B2;
6093 fadeInFactorPrimePrime =
6095 } else {
6096 fadeInFactor = 1.0;
6097 fadeInFactorPrime = 0.0;
6098 fadeInFactorPrimePrime = 0.0;
6099 }
6100
6101 if(zeroPos) {
6102 fadeInFactor = F1;
6103 }
6104
6105
6112
6113 MFloat transitionFactor = F0;
6115 transitionFactor = F0;
6119 transitionFactor = F1;
6121 transitionFactor = (1 + cos((xInit - m_waveOutBeginTransition) / transitionOutLength * pi)) * F1B2;
6122 } else {
6123 transitionFactor = F0;
6124 }
6125
6126 MFloat yTransitionFactor = F1;
6128 yTransitionFactor = F1;
6130 yTransitionFactor = (1 + cos((yInit - m_waveYBeginTransition) / yTransitionLength * pi)) * F1B2;
6131 } else {
6132 yTransitionFactor = F0;
6133 }
6134
6137 const MFloat funcPrime = transitionFactor * yTransitionFactor * (2 * PI * m_waveSpeed / m_waveLength)
6138 * m_waveAmplitude
6143
6144 m_grid->m_coordinates[1][pointId] = func * fadeInFactor + yInit;
6145 m_grid->m_velocity[1][pointId] = func * fadeInFactorPrime + funcPrime * fadeInFactor;
6146 m_grid->m_acceleration[1][pointId] =
6147 funcPrimePrime * fadeInFactor + 2 * funcPrime * fadeInFactorPrime + func * fadeInFactorPrimePrime;
6148 }
6149 }
6150 }
6151
6152 break;
6153 }
6154 case 13: {
6155 // traveling wave on airfoil
6157 if(zeroPos) {
6158 t_offset = F0;
6159 }
6163
6164 MFloat fadeInFactor = 0;
6166
6167 if(t_offset < timeRelaxation) {
6168 fadeInFactor = (1.0 - cos(t_offset / timeRelaxation * pi)) * F1B2;
6169 } else {
6170 fadeInFactor = 1.0;
6171 }
6172
6173 if(zeroPos) {
6174 fadeInFactor = F1;
6175 }
6176
6178
6179 if(myBlockId == 0) {
6181 MFloat transitionFactor = F0;
6186
6190
6192 transitionFactor = F0;
6196 transitionFactor = F1;
6198 transitionFactor = (1 + cos((xWall - m_waveOutBeginTransition) / transitionOutLength * pi)) * F1B2;
6199 } else {
6200 transitionFactor = F0;
6201 }
6202
6203 // linear increase of amplitude
6207
6208 if(wallNormalVec[1] < 0.0) {
6210 inc = 1.0 + heightGrad * (xWall - m_waveEndTransition);
6211 amp0 = m_waveAmplitudePressure;
6212 }
6213
6220
6221
6222 MFloat yTransitionFactor = F1;
6225 yTransitionFactor = F1;
6227 yTransitionFactor = (1 + cos((normalDist - m_waveYBeginTransition) / yTransitionLength * pi)) * F1B2;
6228 } else {
6229 yTransitionFactor = F0;
6230 }
6231
6233 inc * fadeInFactor
6234 * (amp0 * yTransitionFactor * transitionFactor
6237 inc * fadeInFactor
6240
6241 m_grid->m_coordinates[0][pIJK] = xInit + normalDisplacement * wallNormalVec[0];
6242 m_grid->m_coordinates[1][pIJK] = yInit + normalDisplacement * wallNormalVec[1];
6243
6244 m_grid->m_velocity[0][pIJK] = normalVel * wallNormalVec[0];
6245 m_grid->m_velocity[1][pIJK] = normalVel * wallNormalVec[1];
6246 }
6247 }
6248 }
6249 }
6250 break;
6251 }
6252 case 14: {
6253 // oscillating cylinder
6255 if(zeroPos) {
6256 t_offset = F0;
6257 }
6258
6259 MFloat fadeInFactor = 0;
6261
6262 if(t_offset < timeRelaxation) {
6263 fadeInFactor = (1.0 - cos(t_offset / timeRelaxation * pi)) * F1B2;
6264 } else {
6265 fadeInFactor = 1.0;
6266 }
6267
6268 if(zeroPos) {
6269 fadeInFactor = F1;
6270 }
6271
6278 // const MFloat z = m_grid->m_initCoordinates[2][pIJK];
6279
6281
6282 MFloat spaceTransition = F1;
6283
6284 if(r < 1.0) {
6285 spaceTransition = F0;
6286 } else if(r >= 1.0 && r <= 41.0) {
6287 spaceTransition = fabs(r - 1.0) / 40.0;
6288 } else {
6289 spaceTransition = F1;
6290 }
6291
6292 m_grid->m_coordinates[1][pIJK] =
6293 y * spaceTransition
6294 + (1.0 - spaceTransition) * (y + fadeInFactor * m_oscAmplitude * sin(2 * PI * m_oscFreq * t_offset));
6295 m_grid->m_velocity[1][pIJK] =
6296 (1.0 - spaceTransition)
6298 m_grid->m_acceleration[1][pIJK] =
6299 (1.0 - spaceTransition)
6300 * (-fadeInFactor * m_oscAmplitude * POW2(2 * PI * m_oscFreq) * sin(2 * PI * m_oscFreq * t_offset));
6301 }
6302 }
6303 }
6304
6305 break;
6306 }
6307 case 15: {
6308 // streamwise traveling wave on airfoil
6310 if(zeroPos) {
6311 t_offset = F0;
6312 }
6315 const MFloat pressureTransitionLength = m_wavePressureEndTransition - m_wavePressureBeginTransition;
6316 const MFloat pressureTransitionOutLength = m_wavePressureOutEndTransition - m_wavePressureOutBeginTransition;
6318
6319 MFloat fadeInFactor = 0;
6321
6322 if(t_offset < timeRelaxation) {
6323 fadeInFactor = (1.0 - cos(t_offset / timeRelaxation * pi)) * F1B2;
6324 } else {
6325 fadeInFactor = 1.0;
6326 }
6327
6328 if(zeroPos) {
6329 fadeInFactor = F1;
6330 }
6331
6333
6334 if(myBlockId == 0) {
6336 MFloat transitionFactor = F0;
6341
6345
6346 // upper surface
6349 transitionFactor = F0;
6353 transitionFactor = F1;
6355 transitionFactor = (1 + cos((xWall - m_waveOutBeginTransition) / transitionOutLength * pi)) * F1B2;
6356 } else {
6357 transitionFactor = F0;
6358 }
6359
6360 // lower surface
6361 if(wallNormalVec[1] < 0.0) {
6362 amp0 = m_waveAmplitudePressure;
6364 transitionFactor = F0;
6366 transitionFactor =
6369 transitionFactor = F1;
6370 } else if(xWall > m_wavePressureOutBeginTransition && xWall < m_wavePressureOutEndTransition) {
6371 transitionFactor =
6372 (1 + cos((xWall - m_wavePressureOutBeginTransition) / pressureTransitionOutLength * pi)) * F1B2;
6373 } else {
6374 transitionFactor = F0;
6375 }
6376 }
6377
6383
6384 MFloat yTransitionFactor = F1;
6387 yTransitionFactor = F1;
6389 yTransitionFactor = (1 + cos((normalDist - m_waveYBeginTransition) / yTransitionLength * pi)) * F1B2;
6390 } else {
6391 yTransitionFactor = F0;
6392 }
6393
6395 fadeInFactor
6396 * (amp0 * yTransitionFactor * transitionFactor
6399 fadeInFactor
6402
6403 m_grid->m_coordinates[0][pIJK] = xInit + normalDisplacement * wallNormalVec[0];
6404 m_grid->m_coordinates[1][pIJK] = yInit + normalDisplacement * wallNormalVec[1];
6405
6406 m_grid->m_velocity[0][pIJK] = normalVel * wallNormalVec[0];
6407 m_grid->m_velocity[1][pIJK] = normalVel * wallNormalVec[1];
6408 }
6409 }
6410 }
6411 }
6412 break;
6413 }
6414 default: {
6416 }
6417 }
6418
6420
6423 m_grid->extrapolateGhostPointCoordinates();
6424 }
6426 }
6427
6431 }
6432 }
6434
6436 m_grid->computeCellCenterCoordinates();
6438
6440 m_grid->computeMetrics();
6442
6444 m_grid->computeJacobian();
6445
6447}
MFloat m_wavePressureOutEndTransition
Definition: fvstructuredsolver.h:373
MFloat m_waveGradientSuction
Definition: fvstructuredsolver.h:362
MFloat m_waveYBeginTransition
Definition: fvstructuredsolver.h:374
MFloat m_wavePressureBeginTransition
Definition: fvstructuredsolver.h:370
MBool m_mgExchangeCoordinates
Definition: fvstructuredsolver.h:344
MFloat m_waveAmplitudeSuction
Definition: fvstructuredsolver.h:360
MFloat m_wavePressureOutBeginTransition
Definition: fvstructuredsolver.h:372
MFloat m_wavePressureEndTransition
Definition: fvstructuredsolver.h:371
MInt m_gridMovingMethod
Definition: fvstructuredsolver.h:345
MFloat * m_airfoilCoords
Definition: fvstructuredsolver.h:865
MFloat m_waveYEndTransition
Definition: fvstructuredsolver.h:375
MFloat m_waveGradientPressure
Definition: fvstructuredsolver.h:363
MFloat m_waveAmplitudePressure
Definition: fvstructuredsolver.h:361
@ MGCellCenterCoordinates
Definition: fvstructuredtimers.h:40
◆ Muscl()
|
overridevirtual |
Implements FvStructuredSolver< 3 >.
Definition at line 3646 of file fvstructuredsolver3d.cpp.
3646 {
3647 TRACE();
3648
3649
3654 m_grid->saveGrid();
3655 m_grid->saveCellJacobian();
3657 }
3658
3662
3663 // compute the volume fluxes
3668 }
3669
3672 }
3673
3675 (this->*reconstructSurfaceData)();
3677
3681 applySandpaperTripAirfoil();
3682 } else {
3683 applySandpaperTrip();
3684 }
3686 }
3687}
void applySandpaperTripAirfoil()
Definition: fvstructuredsolver3d.cpp:4135
void(FvStructuredSolver3D::* reconstructSurfaceData)()
Definition: fvstructuredsolver3d.h:103
void applyBodyForce(const MBool isRestart, const MBool zeroPos) override
Definition: fvstructuredsolver3d.cpp:7863
void applySandpaperTrip()
Definition: fvstructuredsolver3d.cpp:4077
void moveGrid(const MBool isRestart, const MBool zeroPos) override
Definition: fvstructuredsolver3d.cpp:5798
◆ Muscl_AusmDV()
template<MInt noVars>
| template void FvStructuredSolver3D::Muscl_AusmDV< 7 > | ( | ) |
Definition at line 3330 of file fvstructuredsolver3d.cpp.
3330 {
3331 TRACE();
3332
3333 // stencil identifier
3335
3343
3348
3352
3356
3359#ifdef _OPENMP
3360#pragma omp parallel for
3361#endif
3368 }
3369 }
3370 }
3371 }
3372
3374 }
3375
3376
3378#ifdef _OPENMP
3379#pragma omp parallel for
3380#endif
3383#if defined(MAIA_INTEL_COMPILER)
3384#pragma ivdep
3385#pragma vector always
3386#endif
3392
3396
3399
3400 // unrolled the loop so the compiler
3401 // can optimize better
3407
3413
3419
3425
3431
3436
3439
3444
3445
3446 // compute lenght of metric vector for normalization
3449
3450 // scale by metric length to get velocity in the new basis (get normalized basis vectors)
3452
3453
3455
3456 MFloat AL = FRHOL * PL;
3457 MFloat AR = FRHOR * PR;
3458
3462
3463 AL = sqrt(gamma * AL);
3464 AR = sqrt(gamma * AR);
3465 AL = mMax(AL, AR);
3466 AR = AL;
3467
3470
3471 AL = AL * DGRAD;
3472 AR = AR * DGRAD;
3473
3476
3481 MFloat SV = 2.0 * DGRAD / (m_cells->cellJac[I] + m_cells->cellJac[IP1]) * (FDV + (F1 - FDV) * getPSI(I, dim));
3485
3488
3489 MFloat FXMA = F1B2 * (XMAL1 + fabs(XMAL1));
3491 FXMA = F1B2 * (XMAR1 - fabs(XMAR1));
3492 const MFloat XMARM = ALPHAR * (-F1B4 * POW2(XMAR1 - F1) - FXMA) + FXMA + (mMin(-F1, XMAR) + F1);
3493
3497
3503
3506 UL = UL - UUL2;
3507 UR = UR - UUR2;
3510 VL = VL - UUL3;
3511 VR = VR - UUR3;
3514 WL = WL - UUL4;
3515 WR = WR - UUR4;
3516
3517 flux[CV->RHO_U][I] = RHOU2 * (UL + UR) + ARHOU2 * (UL - UR) + PLR * surf0 + RHOUL * UUL2 + RHOUR * UUR2;
3518 flux[CV->RHO_V][I] = RHOU2 * (VL + VR) + ARHOU2 * (VL - VR) + PLR * surf1 + RHOUL * UUL3 + RHOUR * UUR3;
3519 flux[CV->RHO_W][I] = RHOU2 * (WL + WR) + ARHOU2 * (WL - WR) + PLR * surf2 + RHOUL * UUL4 + RHOUR * UUR4;
3520 flux[CV->RHO_E][I] = RHOU2 * (e0 + e1) + ARHOU2 * (e0 - e1) + PLR * dxdtau;
3521 flux[CV->RHO][I] = RHOU;
3522 }
3523 }
3524 }
3525
3526 // FLUX BALANCE
3528#ifdef _OPENMP
3529#pragma omp parallel for
3530#endif
3533#if defined(MAIA_INTEL_COMPILER)
3534#pragma ivdep
3535#pragma vector always
3536#endif
3542 rhs[I] += flux[v][IM1] - flux[v][I];
3543 }
3544 }
3545 }
3546 }
3547 }
3548}
◆ Muscl_AusmLES()
template<MInt noVars>
| template void FvStructuredSolver3D::Muscl_AusmLES< 7 > | ( | ) |
Definition at line 2980 of file fvstructuredsolver3d.cpp.
2980 {
2981 TRACE();
2982
2985
2993
2997
2998
2999
3002 maia::parallelFor<true, nDim>(beginP0(), endM1(), [=](const MInt& i, const MInt& j, const MInt& k) {
3006 });
3007 }
3008
3010 }
3011
3012
3014 maia::parallelFor<true, nDim>(beginP1(), endM2(), [=](const MInt& i, const MInt& j, const MInt& k) {
3019
3023
3026
3027 // unrolled the loop so the compiler
3028 // can optimize better
3034
3040
3046
3052
3058
3063
3064 // compute length of metric vector for normalization
3067
3068 // scale by metric length to get velocity in the new basis (get normalized basis vectors)
3070
3071
3073
3074
3075 // speed of sound
3078
3081
3084
3087
3090
3093
3094 const MFloat RHOU = F1B2 * (MALR * (RHO_AL + RHO_AR) + fabs(MALR) * (RHO_AL - RHO_AR)) * metricLength;
3096 // multiply by metric length to take surface area into account
3098
3099 flux[CV->RHO_U][I] = RHOU2 * (UL + UR) + AbsRHO_U2 * (UL - UR) + PLR * surf0;
3100 flux[CV->RHO_V][I] = RHOU2 * (VL + VR) + AbsRHO_U2 * (VL - VR) + PLR * surf1;
3101 flux[CV->RHO_W][I] = RHOU2 * (WL + WR) + AbsRHO_U2 * (WL - WR) + PLR * surf2;
3102 flux[CV->RHO_E][I] = RHOU2 * (e0 + e1) + AbsRHO_U2 * (e0 - e1) + PLR * dxdtau;
3103 flux[CV->RHO][I] = RHOU;
3104 });
3105
3106 // FLUX BALANCE
3108 maia::parallelFor<true, nDim>(beginP2(), endM2(), [=](const MInt& i, const MInt& j, const MInt& k) {
3113 rhs[I] += flux[v][IM1] - flux[v][I];
3114 });
3115 }
3116 }
3117}
std::array< MInt, nDim > beginP0()
Definition: fvstructuredsolver.h:280
std::array< MInt, nDim > beginP2()
Definition: fvstructuredsolver.h:294
std::array< MInt, nDim > beginP1()
Definition: fvstructuredsolver.h:287
◆ Muscl_AusmLES_PTHRC()
template<MInt noVars>
| template void FvStructuredSolver3D::Muscl_AusmLES_PTHRC< 7 > | ( | ) |
Definition at line 3124 of file fvstructuredsolver3d.cpp.
3124 {
3125 TRACE();
3126
3127 // stencil identifier
3129
3138
3143
3147
3151
3154#ifdef _OPENMP
3155#pragma omp parallel for
3156#endif
3163 }
3164 }
3165 }
3166 }
3167
3169 }
3170
3172#ifdef _OPENMP
3173#pragma omp parallel for
3174#endif
3177#if defined(MAIA_INTEL_COMPILER)
3178#pragma ivdep
3179#pragma vector always
3180#endif
3186
3190
3193
3194 // unrolled the loop so the compiler
3195 // can optimize better
3201
3207
3213
3219
3225
3230
3231 // compute lenght of metric vector for normalization
3234
3235 // scale by metric length to get velocity in the new basis (get normalized basis vectors)
3237
3238
3240
3241
3242 // speed of sound
3245
3248
3250
3251 // 4th order pressure damping
3256
3261
3266
3270
3274
3276
3279
3282
3285
3286 const MFloat RHOU = F1B2 * (MALR * (RHO_AL + RHO_AR) + fabs(MALR) * (RHO_AL - RHO_AR)) * metricLength;
3288 // multiply by metric length to take surface area into account
3290
3291 flux[CV->RHO_U][I] = RHOU2 * (UL + UR) + AbsRHO_U2 * (UL - UR) + PLR * surf0;
3292 flux[CV->RHO_V][I] = RHOU2 * (VL + VR) + AbsRHO_U2 * (VL - VR) + PLR * surf1;
3293 flux[CV->RHO_W][I] = RHOU2 * (WL + WR) + AbsRHO_U2 * (WL - WR) + PLR * surf2;
3294 flux[CV->RHO_E][I] = RHOU2 * (e0 + e1) + AbsRHO_U2 * (e0 - e1) + PLR * dxdtau;
3295 flux[CV->RHO][I] = RHOU;
3296 }
3297 }
3298 }
3299
3300 // FLUX BALANCE
3302#ifdef _OPENMP
3303#pragma omp parallel for
3304#endif
3307#if defined(MAIA_INTEL_COMPILER)
3308#pragma ivdep
3309#pragma vector always
3310#endif
3316 rhs[I] += flux[v][IM1] - flux[v][I];
3317 }
3318 }
3319 }
3320 }
3321 }
3322}
◆ MusclAlbada()
| void FvStructuredSolver3D::MusclAlbada | ( | ) |
Definition at line 2360 of file fvstructuredsolver3d.cpp.
2360 {
2361 TRACE();
2362 // stencil identifier
2364
2365 // reduce to onedimensional arrays
2371
2373 // MFloat epsi=F1;
2374 // MFloat kappa=F1B3;
2380 // cell ids
2385
2386 // distances q_i+1 - q_i
2388 // distances q_i - q_i-1
2390 const MFloat DSP1 = sqrt(POW2(x[IP2] - x[IP1]) + POW2(y[IP2] - y[IP1]) + POW2(z[IP2] - z[IP1]));
2393
2398
2404
2410 F1B2
2411 - (F1B2
2412 - mMax(F0, (DQP1 * DQM1 * smps + F1B2 * epsLim) / (POW2(DQP1 * DS) + POW2(DQM1 * DSP1) + epsLim)))
2413 * psi;
2414
2415 m_QLeft[var] = pvars[var][I] + DSM * (DSM1 * DQ + DS * DQM1) * phi;
2416 m_QRight[var] = pvars[var][IP1] - DSP * (DS * DQP1 + DSP1 * DQ) * phi;
2417 }
2418
2420 }
2421 }
2422 }
2423
2424 // FLUX BALANCE
2432 }
2433 }
2434 }
2435 }
2436 }
2437}
void AusmLES(MFloat *QLeft, MFloat *QRight, const MInt dim, const MInt cellId)
AUSM central.
Definition: fvstructuredsolver3d.cpp:2692
◆ MusclMinModLimiter()
| void FvStructuredSolver3D::MusclMinModLimiter | ( | ) |
Definition at line 2284 of file fvstructuredsolver3d.cpp.
2284 {
2285 TRACE();
2286 // stencil identifier
2288 // switch for order
2290
2291 // reduce to onedimensional arrays
2296
2299
2304 // cell ids
2309
2310 // distances q_i+1 - q_i
2312 // distances q_i - q_i-1
2314 const MFloat DSP1 = sqrt(POW2(x[IP2] - x[IP1]) + POW2(y[IP2] - y[IP1]) + POW2(z[IP2] - z[IP1]));
2317
2318 if(sword == 2) {
2322
2326
2329
2332
2333 m_QLeft[var] = vars[I] + (DQ * DSM1 + DQM1 * DS) * DSM * phii;
2334 m_QRight[var] = vars[IP1] - (DQP1 * DS + DQ * DSP1) * DSP * phiip;
2335 }
2336 }
2337
2339 }
2340 }
2341 }
2342
2343 // FLUX BALANCE
2349
2352 }
2353 }
2354 }
2355 }
2356 }
2357}
◆ MusclRANS()
| void FvStructuredSolver3D::MusclRANS | ( | ) |
Definition at line 2281 of file fvstructuredsolver3d.cpp.
2281{ m_ransSolver->Muscl(); }
◆ MusclStretched_()
template<FvStructuredSolver3D::fluxmethod ausm, MInt noVars>
| void FvStructuredSolver3D::MusclStretched_ |
Definition at line 3556 of file fvstructuredsolver3d.cpp.
3556 {
3557 TRACE();
3558
3559 // stencil identifier
3569 // MFloat epsi=F1;
3576
3584
3590
3596
3600 // left variables
3603 qleft[v] = vars[I] + f * (f1 * DQ + f2 * DQM1);
3604
3605 // right variables
3608 qright[v] = vars[IP1] - fa * (fb * DQP1 + fc * DQ1);
3609 }
3610
3612 }
3613 }
3614 }
3615
3616 // FLUX BALANCE
3625 rhs[I] += flux[v][IM1] - flux[v][I];
3626 }
3627 }
3628 }
3629 }
3630 }
3631}
◆ MusclVenkatakrishnan3D()
| void FvStructuredSolver3D::MusclVenkatakrishnan3D | ( | ) |
Here, MUSCL and AUSM are run through separately The values of QLeft and QRight (of every cell) from the MUSCL are stored in the scratchspace before passing it to the AUSM scheme. Pros: the limiter can include the values of the neighbouring cells, better results for waves Cons: higher computational effort
- Date
- 2015
Definition at line 2452 of file fvstructuredsolver3d.cpp.
2452 {
2453 TRACE();
2456 // reduce to onedimensional arrays
2462
2466
2467 QLeft.fill(F0);
2468 QRight.fill(F0);
2469
2476
2477 MFloat minNghbrDelta = F0;
2478 MFloat maxNghbrDelta = F0;
2479 MFloat effNghbrDelta = F0;
2480
2481 // 1. get the abs min value of the max and min value from the reconstruction neighbours
2487
2488 MFloat rcnstrctnNghbrValue = F0;
2489 if(rcnstructnNghbr == 0) {
2490 rcnstrctnNghbrValue = pvars[IP1]; //(i/j/k)+1
2491 } else {
2492 rcnstrctnNghbrValue = pvars[IM1]; //(i/j/k)-1
2493 }
2494
2496 maxNghbrDelta = mMax(maxNghbrDelta, tmpDelta);
2497 minNghbrDelta = mMin(minNghbrDelta, tmpDelta);
2498 }
2499 }
2500
2501 effNghbrDelta = mMin(maxNghbrDelta, abs(minNghbrDelta));
2502
2503 MFloat srfcDelta = F0;
2504 MFloat dxEpsSqr = F1;
2509
2510 const MFloat DS = sqrt(POW2(x[IP1] - x[cellId]) + POW2(y[IP1] - y[cellId]) + POW2(z[IP1] - z[cellId]));
2511 // distances q_i - q_i-1
2512 const MFloat DSM1 = sqrt(POW2(x[cellId] - x[IM1]) + POW2(y[cellId] - y[IM1]) + POW2(z[cellId] - z[IM1]));
2514
2515 // 2. get srfcDelta and compute the minimum phi
2517 DSM
2521 (pvars[IP1] - pvars[IM1]) / sqrt(POW2(x[IP1] - x[IM1]) + POW2(y[IP1] - y[IM1]) + POW2(z[IP1] - z[IM1]));
2522 srfcDelta += abs(DQ * dx1);
2523 dxEpsSqr *= dx1;
2524 }
2525
2526 MInt cellPos = 0;
2527
2528 // calling limiter function
2529 (this->*Venkatakrishnan_function)(effNghbrDelta, srfcDelta, dxEpsSqr, cellPos, var, minPhi);
2530
2531 minNghbrDelta = F0;
2532 maxNghbrDelta = F0;
2533 effNghbrDelta = F0;
2534
2535 // 1. get the abs min value of the max and min value from the reconstruction neighbours
2541
2542 MFloat rcnstrctnNghbrValue = F0;
2543 if(rcnstructnNghbr == 0) {
2544 rcnstrctnNghbrValue = pvars[IP2]; //(i/j/k)+2
2545 } else {
2547 }
2548
2550
2551 maxNghbrDelta = mMax(maxNghbrDelta, tmpDelta);
2552 minNghbrDelta = mMin(minNghbrDelta, tmpDelta);
2553 }
2554 }
2555
2556 effNghbrDelta = mMin(maxNghbrDelta, abs(minNghbrDelta));
2557
2558 srfcDelta = F0;
2559 dxEpsSqr = F1;
2564
2565 const MFloat DS = sqrt(POW2(x[IP1] - x[cellId]) + POW2(y[IP1] - y[cellId]) + POW2(z[IP1] - z[cellId]));
2566 // distances q_i - q_i-1
2567 const MFloat DSP1 = sqrt(POW2(x[IP2] - x[IP1]) + POW2(y[IP2] - y[IP1]) + POW2(z[IP2] - z[IP1]));
2569
2570 // 2. get srfcDelta and compute the minimum phi
2571
2573 DSP
2575 + DSP1 * sqrt(POW2(x[IP1] - x[cellId]) + POW2(y[IP1] - y[cellId]) + POW2(z[IP1] - z[cellId])));
2578
2579 srfcDelta += abs(DQ * dx2);
2580 dxEpsSqr *= dx2;
2581 }
2582
2583 cellPos = 1;
2584 (this->*Venkatakrishnan_function)(effNghbrDelta, srfcDelta, dxEpsSqr, cellPos, var, minPhi); // calling Venk
2585
2591
2592 const MFloat DS = sqrt(POW2(x[IP1] - x[cellId]) + POW2(y[IP1] - y[cellId]) + POW2(z[IP1] - z[cellId]));
2593 // distances q_i - q_i-1
2594 const MFloat DSM1 = sqrt(POW2(x[cellId] - x[IM1]) + POW2(y[cellId] - y[IM1]) + POW2(z[cellId] - z[IM1]));
2595 const MFloat DSP1 = sqrt(POW2(x[IP2] - x[IP1]) + POW2(y[IP2] - y[IP1]) + POW2(z[IP2] - z[IP1]));
2598
2599 QLeft(cellId, var, dim1) =
2600 pvars[cellId]
2601 + DSM * (DSM1 * (pvars[IP1] - pvars[cellId]) + DS * (pvars[cellId] - pvars[IM1])) * minPhi(var, 0);
2602 QRight(cellId, var, dim1) =
2603 pvars[IP1]
2604 - DSP * (DS * (pvars[IP2] - pvars[IP1]) + DSP1 * (pvars[IP1] - pvars[cellId])) * minPhi(var, 1);
2605 }
2606 }
2607 }
2608 }
2609 }
2610
2616 // cell ids
2618
2620 m_QLeft[v] = QLeft(cellId, v, dim);
2621 m_QRight[v] = QRight(cellId, v, dim);
2622 }
2623
2625 }
2626 }
2627 }
2628
2629
2630 // FLUX BALANCE
2638 }
2639 }
2640 }
2641 }
2642 }
2643}
void(FvStructuredSolver3D::* Venkatakrishnan_function)(MFloat, MFloat, MFloat, MInt, MInt, MFloatScratchSpace &)
Definition: fvstructuredsolver3d.h:119
◆ nonReflectingBC()
| void FvStructuredSolver3D::nonReflectingBC | ( | ) |
Definition at line 2269 of file fvstructuredsolver3d.cpp.
◆ pointIndex()
Definition at line 5007 of file fvstructuredsolver3d.cpp.
◆ pressure()
Definition at line 5779 of file fvstructuredsolver3d.cpp.
◆ randnormal()
Definition at line 2233 of file fvstructuredsolver3d.cpp.
2233 {
2234 TRACE();
2236 static float storedDeviate; // deviate from previous calculation
2237 MFloat polar, rsquared, var1, var2;
2238
2239 // If no deviate has been stored, the polar Box-Muller transformation is
2240 // performed, producing two independent normally-distributed random
2241 // deviates. One is stored for the next round, and one is returned.
2242 if(!deviateAvailable) {
2243 // choose pairs of uniformly distributed deviates, discarding those
2244 // that don't fall within the unit circle
2245 do {
2248 rsquared = var1 * var1 + var2 * var2;
2250
2251 // calculate polar tranformation for each deviate
2252 polar = sqrt(-2.0 * log(rsquared) / rsquared);
2253
2254 // store first deviate and set flag
2255 storedDeviate = var1 * polar;
2256 deviateAvailable = true;
2257
2258 // return second deviate
2259 return var2 * polar * sigma + mu;
2260
2261 // If a deviate is available from a previous call to this function, it is
2262 // eturned, and the flag is set to false.
2263 } else {
2264 deviateAvailable = false;
2265 return storedDeviate * sigma + mu;
2266 }
2267}
◆ rungeKuttaStep()
|
virtual |
Implements FvStructuredSolver< 3 >.
Definition at line 4418 of file fvstructuredsolver3d.cpp.
4418 {
4419 TRACE();
4424
4430
4431 // set old variables
4437 maia::parallelFor<true>(0, m_noCells, [=](MInt cellId) { oldCellVars[cellId] = cellVars[cellId]; });
4438 }
4439 }
4440
4442 case 2: {
4443 // for moving grids we take the old Jacobian into account
4450
4451 maia::parallelFor<true, nDim>(beginP2(), endM2(), [=](const MInt& i, const MInt& j, const MInt& k) {
4456 });
4457 }
4464
4465 maia::parallelFor<true, nDim>(beginP2(), endM2(), [=](const MInt& i, const MInt& j, const MInt& k) {
4467 cellVars[cellId] =
4469 });
4470 }
4471 } else {
4477
4478 maia::parallelFor<true, nDim>(beginP2(), endM2(), [=](const MInt& i, const MInt& j, const MInt& k) {
4482 });
4483 }
4484 }
4485 break;
4486 }
4487 case 3: {
4493
4494 maia::parallelFor<true, nDim>(beginP2(), endM2(), [=](const MInt& i, const MInt& j, const MInt& k) {
4497 cellVars[cellId] =
4498 rkAlpha * cellVars[cellId] + (F1 - rkAlpha) * oldCellVars[cellId] - factor * cellRhs[cellId];
4499 });
4500 }
4501 break;
4502 }
4503 default: {
4504 stringstream errorMessage;
4505 errorMessage << "Given RungeKutta Order " << m_rungeKuttaOrder << " not implemented! " << endl;
4506 mTerm(1, AT_, errorMessage.str());
4507 }
4508 }
4509
4510 ++m_RKStep;
4511
4515
4516 m_RKStep = 0;
4517
4518 return true;
4519 } else {
4520 return false;
4521 }
4522}
MInt m_rungeKuttaOrder
Definition: fvstructuredsolver.h:560
◆ saveInterpolatedPoints()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 9253 of file fvstructuredsolver3d.cpp.
9253 {
9254 TRACE();
9258
9259 // if it a moving grid the interpolation
9260 // coefficients need to be computed every time
9263 m_intpPointsHasPartnerLocal[pointId] = 0;
9264 m_intpPointsHasPartnerGlobal[pointId] = 0;
9266 m_intpPointsVarsLocal[var][pointId] = F0;
9267 m_intpPointsVarsGlobal[var][pointId] = F0;
9268 }
9269 }
9270
9271 m_pointInterpolation =
9272 make_unique<StructuredInterpolation<3>>(m_nCells, m_cells->coordinates, m_cells->pvariables, m_StructuredComm);
9273 // allocate domains points of the lines
9275 m_intpPointsHasPartnerLocal);
9276 MPI_Allreduce(m_intpPointsHasPartnerLocal, m_intpPointsHasPartnerGlobal, m_intpPointsNoPointsTotal, MPI_INT,
9277 MPI_SUM, m_StructuredComm, AT_, "m_intpPointsHasPartnerLocal", "m_intpPointsHasPartnerGlobal");
9278 }
9279
9280 // calculation of interpolated variables only for the points in domain
9284 m_intpPointsVarsLocal[var][pointId] = m_pointInterpolation->getInterpolatedVariable(pointId, var);
9285 }
9286 }
9287 }
9288
9291 "m_intpPointsVarsLocal[0][0]", "m_intpPointsVarsGlobal[0][0]");
9292
9293 // calculation of right value, if a point is assigned to more than one domain
9297 m_intpPointsVarsGlobal[var][pointId] =
9299 }
9300 }
9301 }
9302
9304 m_pointInterpolation.reset();
9305 }
9306
9310 stringstream fileName;
9312
9314
9317
9319
9321 ParallelIo::size_type dataOffset[2] = {0, 0};
9322 ParallelIo::size_type dataSize[2] = {m_intpPointsNoPoints2D[lineId], m_intpPointsNoPoints[lineId]};
9324
9325 ParallelIo::size_type noDims = 1;
9326 dataSize[0] = dataSize[1];
9328 noDims = 2;
9329 }
9330
9331 stringstream path;
9332 path << lineId;
9334 solutionpath += path.str();
9335 const char* dsetname = solutionpath.c_str();
9336
9338 pio.defineArray(maia::parallel_io::PIO_FLOAT, dsetname, m_pvariableNames[v], noDims, dataSize);
9339 }
9340
9344
9347 pio.writeArray(&m_intpPointsVarsGlobal[v][fieldOffset], dsetname, m_pvariableNames[v], noDims, dataOffset,
9348 dataSize);
9349 }
9350
9351 pio.writeArray(&m_intpPointsCoordinates[0][fieldOffset], dsetname, "x", noDims, dataOffset, dataSize);
9352 pio.writeArray(&m_intpPointsCoordinates[1][fieldOffset], dsetname, "y", noDims, dataOffset, dataSize);
9353 pio.writeArray(&m_intpPointsCoordinates[2][fieldOffset], dsetname, "z", noDims, dataOffset, dataSize);
9354 } else {
9355 dataSize[0] = 0;
9356 dataSize[1] = 0;
9357 MFloat empty = 0;
9359 pio.writeArray(&empty, dsetname, m_pvariableNames[v], 1, dataOffset, dataSize);
9360 }
9361
9362 pio.writeArray(&empty, dsetname, "x", noDims, dataOffset, dataSize);
9363 pio.writeArray(&empty, dsetname, "y", noDims, dataOffset, dataSize);
9364 pio.writeArray(&empty, dsetname, "z", noDims, dataOffset, dataSize);
9365 }
9366 }
9367}
MInt * m_intpPointsOffsets
Definition: fvstructuredsolver.h:499
MInt * m_intpPointsNoPoints2D
Definition: fvstructuredsolver.h:498
MFloat ** m_intpPointsVarsLocal
Definition: fvstructuredsolver.h:507
MString m_intpPointsOutputDir
Definition: fvstructuredsolver.h:667
MInt * m_intpPointsNoPoints
Definition: fvstructuredsolver.h:497
virtual void writePropertiesAsAttributes(ParallelIoHdf5 *pio, MString path)
Overloaded version of writePropertiesAsAttributes that receives ParallelIoHdf5 object pointer instead...
Definition: fvstructuredsolver.cpp:5033
MInt m_intpPointsNoLines
Definition: fvstructuredsolver.h:494
MFloat ** m_intpPointsVarsGlobal
Definition: fvstructuredsolver.h:506
virtual void writeHeaderAttributes(ParallelIoHdf5 *pio, MString fileType)
Overloaded version of writeHeaderAttributes that receives ParallelIoHdf5 object pointer instead of 'f...
Definition: fvstructuredsolver.cpp:5007
◆ saveNodalBoxes()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 8424 of file fvstructuredsolver3d.cpp.
8424 {
8425 stringstream filename;
8426
8427 filename << m_nodalBoxOutputDir << "nodalBoxOutput" << m_outputIterationNumber << m_outputFormat;
8428
8430
8434
8436 m_nodalBoxInterpolation = make_unique<StructuredInterpolation<nDim>>(m_nCells, m_cells->coordinates,
8438
8441 mAlloc(m_nodalBoxLocalDomainOffset, m_nodalBoxNoBoxes, nDim, "m_nodalBoxLocalDomainOffset", 0, AT_);
8444
8445 m_nodalBoxTotalLocalSize = 0;
8459 && m_nOffsetPoints[0]
8461
8464 && m_nodalBoxOffset[b][dim] + m_nodalBoxPoints[b][dim] < m_nOffsetPoints[dim] + m_nActivePoints[dim]) {
8469 m_nodalBoxLocalPoints[b][dim] = (m_nOffsetPoints[dim] + m_nActivePoints[dim]) - m_nodalBoxOffset[b][dim];
8478 } else {
8483 }
8484 }
8485
8489 }
8490
8492 }
8493
8496 }
8497 }
8498
8500 mAlloc(m_nodalBoxCoordinates, nDim, m_nodalBoxTotalLocalSize, "m_nodalBoxCoordinates", 0.0, AT_);
8501 mAlloc(m_nodalBoxVariables, m_maxNoVariables, m_nodalBoxTotalLocalSize, "m_nodalBoxVariables", 0.0, AT_);
8503 }
8504
8508 ++k) {
8511 ++j) {
8514 ++i) {
8521 + (nodalBoxI + (nodalBoxJ + nodalBoxK * m_nodalBoxLocalPoints[b][1]) * m_nodalBoxLocalPoints[b][2]);
8522
8525 }
8526 }
8527 }
8528 }
8529 }
8530
8532 m_nodalBoxPartnerLocal);
8534 }
8535
8536
8538 stringstream pathName;
8540
8544
8548
8550
8551 MInt hasCoordinates = 0;
8552
8553 ParallelIo::size_type size[3] = {m_nodalBoxPoints[b][0], m_nodalBoxPoints[b][1], m_nodalBoxPoints[b][2]};
8556 }
8557
8558 // create datasets for the variables
8560 hasCoordinates = 1;
8564 }
8565 // write output to check if coordinates are contained within the variable list
8566 pio.setAttribute(hasCoordinates, "hasCoordinates", pathName.str());
8567 }
8568
8569
8571 // interpolate primitive variables to nodal grid points
8573 }
8574
8575
8577 // check if the box is contained your inputsolverId
8578 ParallelIo::size_type localOffset[3] = {m_nodalBoxLocalOffset[b][0], m_nodalBoxLocalOffset[b][1],
8580 ParallelIo::size_type localSize[3] = {m_nodalBoxLocalPoints[b][0], m_nodalBoxLocalPoints[b][1],
8582 stringstream pathName;
8586 pio.writeArray(&m_nodalBoxVariables[v][m_nodalBoxTotalLocalOffset[b]], pathName.str(), m_pvariableNames[v],
8587 nDim, localOffset, localSize);
8588 }
8589
8591 pio.writeArray(&m_nodalBoxCoordinates[0][m_nodalBoxTotalLocalOffset[b]], pathName.str(), "x", nDim, localOffset,
8592 localSize);
8593 pio.writeArray(&m_nodalBoxCoordinates[1][m_nodalBoxTotalLocalOffset[b]], pathName.str(), "y", nDim, localOffset,
8594 localSize);
8595 pio.writeArray(&m_nodalBoxCoordinates[2][m_nodalBoxTotalLocalOffset[b]], pathName.str(), "z", nDim, localOffset,
8596 localSize);
8597 }
8598 } else { // write out nothing as box is not contained
8599 ParallelIo::size_type localBoxPoints[3] = {0, 0, 0};
8600 ParallelIo::size_type localBoxOffset[3] = {0, 0, 0};
8601 MFloat empty = 0;
8602
8604 pio.writeArray(&empty, pathName.str(), m_pvariableNames[v], nDim, localBoxOffset, localBoxPoints);
8605 }
8606
8611 }
8612 }
8613 }
8614}
MFloat ** m_nodalBoxVariables
Definition: fvstructuredsolver.h:530
MInt ** m_nodalBoxLocalOffset
Definition: fvstructuredsolver.h:524
MString m_nodalBoxOutputDir
Definition: fvstructuredsolver.h:519
MInt m_nodalBoxTotalLocalSize
Definition: fvstructuredsolver.h:532
MInt ** m_nodalBoxPoints
Definition: fvstructuredsolver.h:523
MInt * m_nodalBoxLocalSize
Definition: fvstructuredsolver.h:527
MInt ** m_nodalBoxLocalDomainOffset
Definition: fvstructuredsolver.h:526
MInt m_outputIterationNumber
Definition: fvstructuredsolver.h:405
std::unique_ptr< StructuredInterpolation< nDim > > m_nodalBoxInterpolation
Definition: fvstructuredsolver.h:517
MInt ** m_nodalBoxLocalPoints
Definition: fvstructuredsolver.h:525
MInt m_nodalBoxNoBoxes
Definition: fvstructuredsolver.h:518
MInt * m_nodalBoxTotalLocalOffset
Definition: fvstructuredsolver.h:533
MBool m_nodalBoxInitialized
Definition: fvstructuredsolver.h:534
MFloat ** m_nodalBoxCoordinates
Definition: fvstructuredsolver.h:529
MBool m_nodalBoxWriteCoordinates
Definition: fvstructuredsolver.h:531
MInt * m_nodalBoxPartnerLocal
Definition: fvstructuredsolver.h:528
MInt ** m_nodalBoxOffset
Definition: fvstructuredsolver.h:522
◆ savePointsToAsciiFile()
|
overridevirtual |
- Date
- 2020
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 2107 of file fvstructuredsolver3d.cpp.
2107 {
2111 m_intpPointsHasPartnerLocal[pointId] = 0;
2112 m_intpPointsHasPartnerGlobal[pointId] = 0;
2114 m_intpPointsVarsLocal[var][pointId] = F0;
2115 m_intpPointsVarsGlobal[var][pointId] = F0;
2116 }
2117 }
2118
2119 m_pointsToAsciiInterpolation = make_unique<StructuredInterpolation<3>>(m_nCells, m_cells->coordinates,
2121 // allocate domains points of the lines
2122 m_pointsToAsciiInterpolation->prepareInterpolation(m_pointsToAsciiNoPoints, m_pointsToAsciiCoordinates,
2124 MPI_Allreduce(m_pointsToAsciiHasPartnerLocal, m_pointsToAsciiHasPartnerGlobal, m_intpPointsNoPointsTotal, MPI_INT,
2126 "m_pointsToAsciiHasPartnerGlobal");
2127 }
2128
2129 // calculation of interpolated variables only for the points in domain
2135 }
2136 }
2137
2138 MPI_Allreduce(MPI_IN_PLACE, &m_pointsToAsciiVars[m_pointsToAsciiCounter][3], m_pointsToAsciiNoPoints, MPI_DOUBLE,
2140
2143 m_pointsToAsciiVars[m_pointsToAsciiCounter][3 + pointId] /= (MFloat)m_pointsToAsciiHasPartnerGlobal[pointId];
2144 }
2145 }
2146
2150
2151 m_pointsToAsciiCounter++;
2153 }
2154
2155
2159 cout << "globalTimeStep: " << globalTimeStep << " writing point data out to ascii file" << endl;
2161 FILE* f_forces;
2162 f_forces = fopen(filename.c_str(), "a+");
2163
2170 }
2171 fprintf(f_forces, "\n");
2172 }
2173 fclose(f_forces);
2174 }
2175
2176 m_pointsToAsciiCounter = 0;
2178 }
2179}
MInt m_pointsToAsciiVarId
Definition: fvstructuredsolver.h:543
MFloat ** m_pointsToAsciiVars
Definition: fvstructuredsolver.h:542
MInt m_pointsToAsciiCounter
Definition: fvstructuredsolver.h:544
MInt m_pointsToAsciiLastComputationStep
Definition: fvstructuredsolver.h:540
MInt m_pointsToAsciiLastOutputStep
Definition: fvstructuredsolver.h:539
◆ scatter()
|
override |
Definition at line 4711 of file fvstructuredsolver3d.cpp.
4712 {
4713 // the ordering of the grid points can be different from
4714 // sending instance ==> reorder it and copy it to the
4715 // right place
4716
4717 for(auto& rcv : rcvComm) {
4720
4721 std::array<MInt, nDim> begin{rcv->startInfoCells[0], rcv->startInfoCells[1], rcv->startInfoCells[2]};
4722 std::array<MInt, nDim> end{rcv->endInfoCells[0], rcv->endInfoCells[1], rcv->endInfoCells[2]};
4723 std::array<MInt, nDim> size{rcv->endInfoCells[0] - rcv->startInfoCells[0],
4724 rcv->endInfoCells[1] - rcv->startInfoCells[1],
4725 rcv->endInfoCells[2] - rcv->startInfoCells[2]};
4727
4728 std::array<MInt, nDim> stepBuffer{0};
4729 std::array<MInt, nDim> startBuffer{0};
4730 std::array<MInt, nDim> endBuffer{0};
4731 std::array<MInt, nDim> sizeBuffer{0};
4732
4734 stepBuffer[rcv->orderInfo[j]] = rcv->stepInfo[j];
4735 }
4736
4738 endBuffer[j] = size[j] - 1;
4739 sizeBuffer[rcv->orderInfo[j]] = size[j];
4740 if(stepBuffer[j] < 0) {
4741 std::swap(startBuffer[j], endBuffer[j]);
4742 }
4743 }
4744
4745 // Aliasing the unique_pointer in rcv to a raw point is needed for PSTL on NVHPC
4746 auto* orderInfo = rcv->orderInfo.data();
4747 auto* startInfoCells = rcv->startInfoCells.data();
4748 auto* cellBuffer = rcv->cellBuffer.get();
4749 auto* variables = &(rcv->variables[0]);
4752 std::array<MInt, nDim> start{};
4753 start[orderInfo[0]] = startBuffer[0] + (i - startInfoCells[0]) * stepBuffer[0];
4754 start[orderInfo[1]] = startBuffer[1] + (j - startInfoCells[1]) * stepBuffer[1];
4755 start[orderInfo[2]] = startBuffer[2] + (k - startInfoCells[2]) * stepBuffer[2];
4756
4757 const MInt bufferId = var * totalSize + start[0] + (start[1] + start[2] * sizeBuffer[1]) * sizeBuffer[0];
4759 variables[var][cellId] = cellBuffer[bufferId];
4760 });
4761 }
4762 }
4763}
◆ shiftAverageCellValues()
|
protected |
- Date
- 2016
Definition at line 9636 of file fvstructuredsolver3d.cpp.
9636 {
9637 TRACE();
9638 MInt cellId_org = 0;
9640 MInt i_new, j_new, k_new;
9641
9642 // accounting for the ghost layers and shift the values to the right place
9647 i_new = i + m_noGhostLayers;
9648 j_new = j + m_noGhostLayers;
9649 k_new = k + m_noGhostLayers;
9651
9654 m_summedVars[var][cellId_org] = F0;
9655 }
9656
9660 m_favre[var][cellId_org] = F0;
9661 }
9662 }
9663 }
9664 }
9665 }
9666}
◆ shiftAverageCellValuesRestart()
|
protected |
- Date
- 12.01.2015
Definition at line 9580 of file fvstructuredsolver3d.cpp.
9580 {
9581 TRACE();
9582 MInt cellId_org = 0;
9584 MInt i_new, j_new, k_new;
9585
9586 // accounting for the ghost layers and shift the values to the right place
9591 i_new = i + m_noGhostLayers;
9592 j_new = j + m_noGhostLayers;
9593 k_new = k + m_noGhostLayers;
9595
9598 m_summedVars[var][cellId_org] = F0;
9599 }
9600
9604 m_favre[var][cellId_org] = F0;
9605 }
9606 }
9607
9610 m_square[var][cellId_org] = F0;
9611 }
9612
9616 m_cube[var][cellId_org] = F0;
9617 }
9618 }
9619
9623 m_fourth[var][cellId_org] = F0;
9624 }
9625 }
9626 }
9627 }
9628 }
9629}
MInt getNoPPSquareVars()
Returns number of pp Square variables.
Definition: fvstructuredpostprocessing.cpp:2412
MFloat ** m_square
Definition: fvstructuredpostprocessing.h:94
MFloat ** m_fourth
Definition: fvstructuredpostprocessing.h:96
MFloat ** m_cube
Definition: fvstructuredpostprocessing.h:95
◆ spanwiseAvgZonal()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 4303 of file fvstructuredsolver3d.cpp.
4303 {
4304 if(!variables.empty()) {
4309
4314 MInt localCellId3D =
4315 i + m_noGhostLayers + ((j + m_noGhostLayers) + (k + m_noGhostLayers) * m_nCells[1]) * m_nCells[2];
4316 MInt globalCellId2D = (i + m_nOffsetCells[2]) + (j + m_nOffsetCells[1]) * m_grid->getMyBlockNoCells(2);
4317 localSpannwiseVars(globalCellId2D, varPos) += variables[varPos][localCellId3D];
4318 }
4319 }
4320 }
4321 }
4322
4323 MPI_Allreduce(&localSpannwiseVars(0, 0), &globalSpannwiseVars(0, 0), totalNoCellsIJ * noVars, MPI_DOUBLE, MPI_SUM,
4325
4328 globalSpannwiseVars(cellId, varPos) /= m_grid->getMyBlockNoCells(0);
4329 }
4330 }
4331
4336 MInt localCellId3D =
4337 i + m_noGhostLayers + ((j + m_noGhostLayers) + (k + m_noGhostLayers) * m_nCells[1]) * m_nCells[2];
4338 MInt globalCellId2D = (i + m_nOffsetCells[2]) + (j + m_nOffsetCells[1]) * m_grid->getMyBlockNoCells(2);
4339 variables[varPos][localCellId3D] = globalSpannwiseVars(globalCellId2D, varPos);
4340 }
4341 }
4342 }
4343 }
4344 }
4345}
◆ spanwiseWaveReorder()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 4872 of file fvstructuredsolver3d.cpp.
4872 {
4874
4876 MInt waveZeroPos = ((MInt)round((globalTimeStep - m_movingGridStepOffset) / m_waveNoStepsPerCell)) % allCellsK;
4877
4879 waveZeroPos =
4880 allCellsK - ((MInt)round((globalTimeStep - m_movingGridStepOffset) / m_waveNoStepsPerCell)) % allCellsK;
4881 }
4882
4883 m_windowInfo->createWaveWindowMapping(waveZeroPos);
4884
4887 noVars += (2 * nDim - 3);
4888 computeVorticity();
4889 }
4890
4892 waveExchange();
4893
4895}
void waveExchange() override
Definition: fvstructuredsolver3d.cpp:4765
void computeVorticity() override
Definition: fvstructuredsolver3d.cpp:9093
MInt m_movingGridStepOffset
Definition: fvstructuredsolver.h:346
std::vector< std::unique_ptr< StructuredComm< nDim > > > m_waveRcvComm
Definition: fvstructuredsolver.h:734
std::vector< std::unique_ptr< StructuredComm< nDim > > > m_waveSndComm
Definition: fvstructuredsolver.h:733
MInt m_waveNoStepsPerCell
Definition: fvstructuredsolver.h:348
@ WaveSpanwiseReordering
Definition: fvstructuredtimers.h:46
◆ surfId()
◆ tripForceCoefficients()
| void FvStructuredSolver3D::tripForceCoefficients | ( | MFloat * | modes, |
| MFloat * | forceCoef, | ||
| MFloat * | coords, | ||
| MInt | noCells, | ||
| MInt | noModes | ||
| ) |
Definition at line 4214 of file fvstructuredsolver3d.cpp.
4215 {
4216 MFloat maxLocalValue = -999999.0;
4217 MFloat minLocalValue = 999999.0;
4218 MFloat* ak = &modes[0];
4219 MFloat* phik = &modes[noModes];
4220
4224 forceCoef[k] += sin(z * ak[n] + phik[n]);
4225 }
4226
4227 maxLocalValue = mMax(maxLocalValue, forceCoef[k]);
4228 minLocalValue = mMin(minLocalValue, forceCoef[k]);
4229 }
4230
4231 // find out min and max to normalize coefficients to interval [-1,1]
4232 MFloat maxGlobalValue = 0.0;
4233 MFloat minGlobalValue = 0.0;
4234 MPI_Allreduce(&maxLocalValue, &maxGlobalValue, 1, MPI_DOUBLE, MPI_MAX, m_StructuredComm, AT_, "maxLocalValue",
4235 "maxGlobalValue");
4236 MPI_Allreduce(&minLocalValue, &minGlobalValue, 1, MPI_DOUBLE, MPI_MIN, m_StructuredComm, AT_, "minLocalValue",
4237 "minGlobalValue");
4238
4239 // normalize the series
4241 forceCoef[k] = 2 * (forceCoef[k] - minGlobalValue) / (maxGlobalValue - minGlobalValue) - 1.0;
4242 }
4243}
◆ tripFourierCoefficients()
| void FvStructuredSolver3D::tripFourierCoefficients | ( | MFloat * | modes, |
| MInt | noModes, | ||
| MFloat | maxWaveLength, | ||
| MFloat | minWaveLength | ||
| ) |
Definition at line 4246 of file fvstructuredsolver3d.cpp.
◆ updateSpongeLayer()
|
virtual |
Implements FvStructuredSolver< 3 >.
Definition at line 4412 of file fvstructuredsolver3d.cpp.
◆ VENKATAKRISHNAN_FCT()
| void FvStructuredSolver3D::VENKATAKRISHNAN_FCT | ( | MFloat | effNghbrDelta, |
| MFloat | srfcDelta, | ||
| MFloat | dxEpsSqr, | ||
| MInt | cellPos, | ||
| MInt | var, | ||
| MFloatScratchSpace & | minPhi | ||
| ) |
- Date
- 2015
Definition at line 2664 of file fvstructuredsolver3d.cpp.
◆ VENKATAKRISHNAN_MOD_FCT()
| void FvStructuredSolver3D::VENKATAKRISHNAN_MOD_FCT | ( | MFloat | effNghbrDelta, |
| MFloat | srfcDelta, | ||
| MFloat | dxEpsSqr, | ||
| MInt | cellPos, | ||
| MInt | var, | ||
| MFloatScratchSpace & | minPhi | ||
| ) |
- Date
- 2015
Definition at line 2652 of file fvstructuredsolver3d.cpp.
2653 {
2655 minPhi(var, cellPos) = (pow(effNghbrDelta, F2) + epsSqr + F2 * effNghbrDelta * srfcDelta)
2656 / (pow(effNghbrDelta, F2) + F2 * pow(srfcDelta, F2) + effNghbrDelta * srfcDelta + epsSqr);
2657}
◆ viscousFlux()
|
virtual |
Implements FvStructuredSolver< 3 >.
Definition at line 5012 of file fvstructuredsolver3d.cpp.
5012{ (this->*viscFluxMethod)(); }
void(FvStructuredSolver3D::* viscFluxMethod)()
Definition: fvstructuredsolver3d.h:129
◆ viscousFluxCorrection()
| void FvStructuredSolver3D::viscousFluxCorrection | ( | ) |
Definition at line 5308 of file fvstructuredsolver3d.cpp.
5308 {
5313
5319
5323
5324 MInt dim = 0;
5325 MInt start[3], end[3], nghbr[20];
5326 MInt len1[3];
5327 MInt totalCells;
5328
5330 // only correct for bc 6000 not for bc 4000-5000
5332 totalCells = 1;
5335 if(len1[j] != 0) totalCells *= len1[j];
5336 }
5337
5340 dim = n;
5341 // start[n]=m_singularity[i].start[n]+1;
5344 } else {
5347 }
5348 }
5349
5350 MFloat u[20], v[20], w[20], T[20];
5351 MFloat U, V, W, t, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, dTdx, dTdy, dTdz;
5352 MFloat qx, qy, qz;
5353 MFloat tau1, tau2, tau3, tau4, tau5, tau6;
5354 MFloat mueOverRe, mue, mueH;
5358 MInt count = 0;
5359 MInt temp[3] = {0, 0, 0};
5361 kk + m_singularity[i].Viscous[2]);
5362
5364 m_cells->cornerMetrics[0][IJK], m_cells->cornerMetrics[1][IJK], m_cells->cornerMetrics[2][IJK],
5365 m_cells->cornerMetrics[3][IJK], m_cells->cornerMetrics[4][IJK], m_cells->cornerMetrics[5][IJK],
5366 m_cells->cornerMetrics[6][IJK], m_cells->cornerMetrics[7][IJK], m_cells->cornerMetrics[8][IJK]};
5367
5368 temp[dim] = 1;
5369 nghbr[count++] = cellIndex(ii, jj, kk);
5370 nghbr[count++] = cellIndex(ii + temp[0], jj + temp[1], kk + temp[2]);
5371
5374 nghbr[count++] = cellIndex(ii + change[0], jj + change[1], kk + change[2]);
5375 nghbr[count++] = cellIndex(ii + temp[0] + change[0], jj + temp[1] + change[1], kk + temp[2] + change[2]);
5376 }
5377
5379 cout << "what the hell! it is wrong!!!" << endl;
5380 }
5381
5383 u[m] = rhou[nghbr[m]] / rho[nghbr[m]];
5384 v[m] = rhov[nghbr[m]] / rho[nghbr[m]];
5385 w[m] = rhow[nghbr[m]] / rho[nghbr[m]];
5386 T[m] = (m_gamma * gammaMinusOne
5388 / rho[nghbr[m]];
5389 }
5390
5391 U = F0;
5392 V = F0;
5393 W = F0;
5394 t = F0;
5395 dudx = F0;
5396 dudy = F0;
5397 dudz = F0;
5398 dvdx = F0;
5399 dvdy = F0;
5400 dvdz = F0;
5401 dwdx = F0;
5402 dwdy = F0;
5403 dwdz = F0;
5404 dTdx = F0;
5405 dTdy = F0;
5406 dTdz = F0;
5407
5408 MInt id2 = ii - start[0] + ((jj - start[1]) + (kk - start[2]) * len1[1]) * len1[0];
5409
5411 MInt ID = id2 * count + n;
5416
5421
5426
5431 }
5432
5433 tau1 = 2 * dudx - 2 / 3 * (dudx + dvdy + dwdz);
5434 tau2 = dudy + dvdx;
5435 tau3 = dudz + dwdx;
5436 tau4 = 2 * dvdy - 2 / 3 * (dudx + dvdy + dwdz);
5437 tau5 = dvdz + dwdy;
5438 tau6 = 2 * dwdz - 2 / 3 * (dudx + dvdy + dwdz);
5439
5440 mue = SUTHERLANDLAW(t);
5441 mueOverRe = mue * rRe;
5442 tau1 *= mueOverRe;
5443 tau2 *= mueOverRe;
5444 tau3 *= mueOverRe;
5445 tau4 *= mueOverRe;
5446 tau5 *= mueOverRe;
5447 tau6 *= mueOverRe;
5448 mueH = FgammaMinusOne * mueOverRe * rPr;
5449
5450 qx = mueH * dTdx + U * tau1 + V * tau2 + W * tau3;
5451 qy = mueH * dTdy + U * tau2 + V * tau4 + W * tau5;
5452 qz = mueH * dTdz + U * tau3 + V * tau5 + W * tau6;
5453
5454
5455 // efluxes
5458
5461
5464
5467
5468 // ffluxes
5471
5474
5477
5480
5481 // gfluxes
5484
5487
5490
5493 }
5494 }
5495 }
5496 }
5497 }
5498}
◆ viscousFluxLES()
template<MBool twoEqRans>
| void FvStructuredSolver3D::viscousFluxLES |
Definition at line 5020 of file fvstructuredsolver3d.cpp.
5020 {
5021 TRACE();
5027
5036
5041
5042 // only relevant for 2-eq Rans model (Waiting for if constexpr)
5043 const MFloat* const RESTRICT TKE = (twoEqRans) ? ALIGNED_F(m_cells->pvariables[PV->RANS_VAR[0]]) : nullptr;
5044
5045 maia::parallelFor<true, nDim>(beginP1(), endM1(), [=](const MInt& i, const MInt& j, const MInt& k) {
5048 muLam[I] = SUTHERLANDLAW(T[I]);
5049 });
5050
5051 maia::parallelFor<true, nDim>(beginP1(), endM1(), [=](const MInt& i, const MInt& j, const MInt& k) {
5052 // get the adjacent cells;
5061
5063 m_cells->cornerMetrics[0][IJK], m_cells->cornerMetrics[1][IJK], m_cells->cornerMetrics[2][IJK],
5064 m_cells->cornerMetrics[3][IJK], m_cells->cornerMetrics[4][IJK], m_cells->cornerMetrics[5][IJK],
5065 m_cells->cornerMetrics[6][IJK], m_cells->cornerMetrics[7][IJK], m_cells->cornerMetrics[8][IJK]};
5066
5067
5068 const MFloat dudxi = F1B4 * (u[IPJPKP] + u[IPJPK] + u[IPJKP] + u[IPJK] - u[IJPKP] - u[IJPK] - u[IJKP] - u[IJK]);
5069 const MFloat dudet = F1B4 * (u[IPJPKP] + u[IJPKP] + u[IPJPK] + u[IJPK] - u[IPJKP] - u[IJKP] - u[IPJK] - u[IJK]);
5070 const MFloat dudze = F1B4 * (u[IPJPKP] + u[IJPKP] + u[IPJKP] + u[IJKP] - u[IPJPK] - u[IJPK] - u[IPJK] - u[IJK]);
5071
5072 const MFloat dvdxi = F1B4 * (v[IPJPKP] + v[IPJPK] + v[IPJKP] + v[IPJK] - v[IJPKP] - v[IJPK] - v[IJKP] - v[IJK]);
5073 const MFloat dvdet = F1B4 * (v[IPJPKP] + v[IJPKP] + v[IPJPK] + v[IJPK] - v[IPJKP] - v[IJKP] - v[IPJK] - v[IJK]);
5074 const MFloat dvdze = F1B4 * (v[IPJPKP] + v[IJPKP] + v[IPJKP] + v[IJKP] - v[IPJPK] - v[IJPK] - v[IPJK] - v[IJK]);
5075
5076 const MFloat dwdxi = F1B4 * (w[IPJPKP] + w[IPJPK] + w[IPJKP] + w[IPJK] - w[IJPKP] - w[IJPK] - w[IJKP] - w[IJK]);
5077 const MFloat dwdet = F1B4 * (w[IPJPKP] + w[IJPKP] + w[IPJPK] + w[IJPK] - w[IPJKP] - w[IJKP] - w[IPJK] - w[IJK]);
5078 const MFloat dwdze = F1B4 * (w[IPJPKP] + w[IJPKP] + w[IPJKP] + w[IJKP] - w[IPJPK] - w[IJPK] - w[IPJK] - w[IJK]);
5079
5080 const MFloat dTdxi = F1B4 * (T[IPJPKP] + T[IPJPK] + T[IPJKP] + T[IPJK] - T[IJPKP] - T[IJPK] - T[IJKP] - T[IJK]);
5081 const MFloat dTdet = F1B4 * (T[IPJPKP] + T[IJPKP] + T[IPJPK] + T[IJPK] - T[IPJKP] - T[IJKP] - T[IPJK] - T[IJK]);
5082 const MFloat dTdze = F1B4 * (T[IPJPKP] + T[IJPKP] + T[IPJKP] + T[IJKP] - T[IPJPK] - T[IJPK] - T[IPJK] - T[IJK]);
5083
5084 const MFloat uAvg = F1B8 * (u[IPJPKP] + u[IJPKP] + u[IJPK] + u[IPJPK] + u[IPJKP] + u[IJKP] + u[IJK] + u[IPJK]);
5085 const MFloat vAvg = F1B8 * (v[IPJPKP] + v[IJPKP] + v[IJPK] + v[IPJPK] + v[IPJKP] + v[IJKP] + v[IJK] + v[IPJK]);
5086 const MFloat wAvg = F1B8 * (w[IPJPKP] + w[IJPKP] + w[IJPK] + w[IPJPK] + w[IPJKP] + w[IJKP] + w[IJK] + w[IPJK]);
5087
5089 * (muLam[IPJPKP] + muLam[IJPKP] + muLam[IJPK] + muLam[IPJPK] + muLam[IPJKP] + muLam[IJKP]
5090 + muLam[IJK] + muLam[IPJK]);
5091
5092 // turbulent viscosity is set to zero for LES
5093 // and is only non-zero for RANS
5095 * (muTurb[IPJPKP] + muTurb[IJPKP] + muTurb[IJPK] + muTurb[IPJPK] + muTurb[IPJKP]
5096 + muTurb[IJKP] + muTurb[IJK] + muTurb[IPJK]);
5097
5098 MFloat TKEcorner = F0;
5099 if(twoEqRans)
5100 TKEcorner =
5101 -2 / 3 * F1B8
5102 * (rho[IPJPKP] * TKE[IPJPKP] + rho[IJPKP] * TKE[IJPKP] + rho[IJPK] * TKE[IJPK] + rho[IPJPK] * TKE[IPJPK]
5103 + rho[IPJKP] * TKE[IPJKP] + rho[IJKP] * TKE[IJKP] + rho[IJK] * TKE[IJK] + rho[IPJK] * TKE[IPJK]);
5104
5105 // compute tau1 = 2 du/dx - 2/3 ( du/dx + dv/dy + dw/dz )
5106
5107 // tau_xx = 4/3*( du/dxi * dxi/dx + du/deta * deta/dx + du/dzeta * dzeta/dx )
5108 // - 2/3*( dv/dxi * dxi/dy + dv/deta * deta/dy + dv/dzeta * dzeta/dy)
5109 // - 2/3*( dw/dxi * dxi/dz + dw/deta * deta/dz + dw/dzeta * dzeta/dz )
5110 MFloat tau1 = F4B3
5113 -
5114
5115 F2B3
5118 -
5119
5120 F2B3
5123
5124 // compute tau2 = du/dy + dv/dx
5125
5126 // tau_xy = du/dxi * dxi/dy + du/deta * deta/dy + du/dzeta * dzeta/dy
5127 // + dv/dxi * dxi/dx + dv/deta * deta/dx + dv/dzeta * dzeta/dx
5130
5133
5134 // compute tau3 = du/dz + dw/dx
5135
5136 // tau_xz = du/dxi * dxi/dz + du/deta * deta/dz + du/dzeta * dzeta/dz
5137 // + dw/dxi * dxi/dx + dw/deta * deta/dx + dw/dzeta * dzeta/dx
5140
5143
5144 // compute tau4 = 2 dv/dy - 2/3 ( du/dx + dv/dy + dw/dz )
5145
5146 // tau_yy = 4/3*( dv/dxi * dxi/dy + dv/deta * deta/dy + dv/dzeta * dzeta/dy )
5147 // - 2/3*( du/dxi * dxi/dx + du/deta * deta/dx + du/dzeta * dzeta/dx)
5148 // - 2/3*( dw/dxi * dxi/dz + dw/deta * deta/dz + dw/dzeta * dzeta/dz )
5149 MFloat tau4 = F4B3
5152 -
5153
5154 F2B3
5157 -
5158
5159 F2B3
5162
5163 // compute tau5 = dv/dz + dw/dy
5164
5165 // tau_yz = dv/dxi * dxi/dz + dv/deta * deta/dz + dv/dzeta * dzeta/dz
5166 // + dw/dxi * dxi/dy + dw/deta * deta/dy + dw/dzeta * dzeta/dy
5169
5172
5173 // compute tau6 = 2 dw/dz - 2/3 ( du/dx + dv/dy + dw/dz )
5174
5175 // tau_zz = 4/3*( dw/dxi * dxi/dz + dw/deta * deta/dz + dw/dzeta * dzeta/dz )
5176 // - 2/3*( du/dxi * dxi/dx + du/deta * deta/dx + du/dzeta * dzeta/dx)
5177 // - 2/3*( dv/dxi * dxi/dy + dv/deta * deta/dy + dv/dzeta * dzeta/dy)
5178 MFloat tau6 = F4B3
5181 -
5182
5183 F2B3
5186 -
5187
5188 F2B3
5191
5192
5193 const MFloat dTdx = dTdxi * cornerMetrics[xsd * 3 + xsd] + dTdet * cornerMetrics[ysd * 3 + xsd]
5195
5196 const MFloat dTdy = dTdxi * cornerMetrics[xsd * 3 + ysd] + dTdet * cornerMetrics[ysd * 3 + ysd]
5198
5199 const MFloat dTdz = dTdxi * cornerMetrics[xsd * 3 + zsd] + dTdet * cornerMetrics[ysd * 3 + zsd]
5201
5204 tau1 = mueOverRe * tau1 + TKEcorner;
5205 tau2 *= mueOverRe;
5206 tau3 *= mueOverRe;
5207 tau4 = mueOverRe * tau4 + TKEcorner;
5208 tau5 *= mueOverRe;
5209 tau6 = mueOverRe * tau6 + TKEcorner;
5210
5211 const MFloat muCombined = FgammaMinusOne * rRe * fJac * (rPrLam * muLamAvg + rPrTurb * muTurbAvg);
5212
5216
5217 // efluxes
5218 eflux[0][IJK] =
5219 tau1 * cornerMetrics[xsd * 3 + xsd] + tau2 * cornerMetrics[xsd * 3 + ysd] + tau3 * cornerMetrics[xsd * 3 + zsd];
5220
5221 eflux[1][IJK] =
5222 tau2 * cornerMetrics[xsd * 3 + xsd] + tau4 * cornerMetrics[xsd * 3 + ysd] + tau5 * cornerMetrics[xsd * 3 + zsd];
5223
5224 eflux[2][IJK] =
5225 tau3 * cornerMetrics[xsd * 3 + xsd] + tau5 * cornerMetrics[xsd * 3 + ysd] + tau6 * cornerMetrics[xsd * 3 + zsd];
5226
5227 eflux[3][IJK] =
5228 qx * cornerMetrics[xsd * 3 + xsd] + qy * cornerMetrics[xsd * 3 + ysd] + qz * cornerMetrics[xsd * 3 + zsd];
5229
5230 // ffluxes
5231 fflux[0][IJK] =
5232 tau1 * cornerMetrics[ysd * 3 + xsd] + tau2 * cornerMetrics[ysd * 3 + ysd] + tau3 * cornerMetrics[ysd * 3 + zsd];
5233
5234 fflux[1][IJK] =
5235 tau2 * cornerMetrics[ysd * 3 + xsd] + tau4 * cornerMetrics[ysd * 3 + ysd] + tau5 * cornerMetrics[ysd * 3 + zsd];
5236
5237 fflux[2][IJK] =
5238 tau3 * cornerMetrics[ysd * 3 + xsd] + tau5 * cornerMetrics[ysd * 3 + ysd] + tau6 * cornerMetrics[ysd * 3 + zsd];
5239
5240 fflux[3][IJK] =
5241 qx * cornerMetrics[ysd * 3 + xsd] + qy * cornerMetrics[ysd * 3 + ysd] + qz * cornerMetrics[ysd * 3 + zsd];
5242
5243 // gfluxes
5244 gflux[0][IJK] =
5245 tau1 * cornerMetrics[zsd * 3 + xsd] + tau2 * cornerMetrics[zsd * 3 + ysd] + tau3 * cornerMetrics[zsd * 3 + zsd];
5246
5247 gflux[1][IJK] =
5248 tau2 * cornerMetrics[zsd * 3 + xsd] + tau4 * cornerMetrics[zsd * 3 + ysd] + tau5 * cornerMetrics[zsd * 3 + zsd];
5249
5250 gflux[2][IJK] =
5251 tau3 * cornerMetrics[zsd * 3 + xsd] + tau5 * cornerMetrics[zsd * 3 + ysd] + tau6 * cornerMetrics[zsd * 3 + zsd];
5252
5253 gflux[3][IJK] =
5254 qx * cornerMetrics[zsd * 3 + xsd] + qy * cornerMetrics[zsd * 3 + ysd] + qz * cornerMetrics[zsd * 3 + zsd];
5255 });
5256
5257
5258 // viscous flux correction for the singular points
5259 // m_hasSingularity=0 means no singular points in this block, otherwise do flux correction
5261 viscousFluxCorrection();
5262 }
5263
5266 maia::parallelFor<true, nDim>(beginIM1, endM2(), [=](const MInt& i, const MInt& j, const MInt& k) {
5271
5272 vflux[0][IJK] = F1B4 * (eflux[var][IJK] + eflux[var][IJKM] + eflux[var][IJMK] + eflux[var][IJMKM]);
5273 });
5274
5275
5277 maia::parallelFor<true, nDim>(beginJM1, endM2(), [=](const MInt& i, const MInt& j, const MInt& k) {
5282
5283 vflux[1][IJK] = F1B4 * (fflux[var][IJK] + fflux[var][IJKM] + fflux[var][IMJK] + fflux[var][IMJKM]);
5284 });
5285
5287 maia::parallelFor<true, nDim>(beginKM1, endM2(), [=](const MInt& i, const MInt& j, const MInt& k) {
5292
5293 vflux[2][IJK] = F1B4 * (gflux[var][IJK] + gflux[var][IMJK] + gflux[var][IJMK] + gflux[var][IMJMK]);
5294 });
5295
5296 maia::parallelFor<true, nDim>(beginP2(), endM2(), [=](const MInt& i, const MInt& j, const MInt& k) {
5302 vflux[0][IJK] - vflux[0][IMJK] + vflux[1][IJK] - vflux[1][IJMK] + vflux[2][IJK] - vflux[2][IJKM];
5303 });
5304 }
5305}
void viscousFluxCorrection()
Definition: fvstructuredsolver3d.cpp:5308
◆ viscousFluxRANS()
| void FvStructuredSolver3D::viscousFluxRANS | ( | ) |
Definition at line 5014 of file fvstructuredsolver3d.cpp.
5014{ m_ransSolver->viscousFluxRANS(); }
◆ waveExchange()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 4765 of file fvstructuredsolver3d.cpp.
4765 {
4766 std::vector<MPI_Request> sndRequests;
4767 std::vector<MPI_Request> rcvRequests;
4768 std::vector<MPI_Status> sndStatus;
4769 std::vector<MPI_Status> rcvStatus;
4770 sndRequests.reserve(m_waveSndComm.size());
4771 rcvRequests.reserve(m_waveRcvComm.size());
4772
4773 waveGather();
4774 waveSend(sndRequests);
4775 waveReceive(rcvRequests);
4776
4777 sndStatus.resize(sndRequests.size());
4778 MPI_Waitall(sndRequests.size(), &sndRequests[0], &sndStatus[0], AT_);
4779 rcvStatus.resize(rcvRequests.size());
4780 MPI_Waitall(rcvRequests.size(), &rcvRequests[0], &rcvStatus[0], AT_);
4781
4782 waveScatter();
4783}
void waveSend(std::vector< MPI_Request > &)
Definition: fvstructuredsolver3d.cpp:4818
void waveReceive(std::vector< MPI_Request > &)
Definition: fvstructuredsolver3d.cpp:4829
◆ waveGather()
| void FvStructuredSolver3D::waveGather | ( | ) |
Definition at line 4786 of file fvstructuredsolver3d.cpp.
4786 {
4788 MInt pos = 0;
4789
4796 pos++;
4797 }
4798 }
4799 }
4800 }
4801
4809 pos++;
4810 }
4811 }
4812 }
4813 }
4814 }
4815 }
4816}
◆ waveReceive()
| void FvStructuredSolver3D::waveReceive | ( | std::vector< MPI_Request > & | rcvRequests | ) |
Definition at line 4829 of file fvstructuredsolver3d.cpp.
◆ waveScatter()
| void FvStructuredSolver3D::waveScatter | ( | ) |
Definition at line 4840 of file fvstructuredsolver3d.cpp.
4840 {
4842 MInt pos = 0;
4843
4850 pos++;
4851 }
4852 }
4853 }
4854 }
4855
4863 pos++;
4864 }
4865 }
4866 }
4867 }
4868 }
4869 }
4870}
MFloat ** m_tempWaveSample
Definition: fvstructuredpostprocessing.h:97
◆ waveSend()
| void FvStructuredSolver3D::waveSend | ( | std::vector< MPI_Request > & | sndRequests | ) |
Definition at line 4818 of file fvstructuredsolver3d.cpp.
◆ zonalAllreduce()
| void FvStructuredSolver3D::zonalAllreduce | ( | ) |
◆ zonalBufferAverage()
| void FvStructuredSolver3D::zonalBufferAverage | ( | ) |
◆ zonalCheck()
| void FvStructuredSolver3D::zonalCheck | ( | ) |
◆ zonalExchange()
|
overridevirtual |
Reimplemented from FvStructuredSolver< 3 >.
Definition at line 546 of file fvstructuredsolver3d.cpp.
◆ zonalGather()
| void FvStructuredSolver3D::zonalGather | ( | ) |
Definition at line 554 of file fvstructuredsolver3d.cpp.
554 {
555 // gather sending variables in buffer
559 MInt pos = 0;
561 if(zBC->localCommReceiverIds[cellId] == zBC->sndComm[noSnd]->localId) {
564 zBC->sndComm[noSnd]->buffer[pos + 0 * noCells] =
566 zBC->sndComm[noSnd]->buffer[pos + 1 * noCells] =
568 zBC->sndComm[noSnd]->buffer[pos + 2 * noCells] =
570 zBC->sndComm[noSnd]->buffer[pos + 3 * noCells] =
572 zBC->sndComm[noSnd]->buffer[pos + 4 * noCells] =
574 zBC->sndComm[noSnd]->buffer[pos + 5 * noCells] =
576 } else {
577 zBC->sndComm[noSnd]->buffer[pos + 0 * noCells] =
579 zBC->sndComm[noSnd]->buffer[pos + 1 * noCells] =
581 zBC->sndComm[noSnd]->buffer[pos + 2 * noCells] =
583 zBC->sndComm[noSnd]->buffer[pos + 3 * noCells] =
585 zBC->sndComm[noSnd]->buffer[pos + 4 * noCells] =
587 zBC->sndComm[noSnd]->buffer[pos + 5 * noCells] =
589 }
590 pos++;
591 }
592 }
593 }
594 }
595}
◆ zonalRealInterpolation()
| void FvStructuredSolver3D::zonalRealInterpolation | ( | ) |
◆ zonalReceive()
| void FvStructuredSolver3D::zonalReceive | ( | ) |
Definition at line 613 of file fvstructuredsolver3d.cpp.
613 {
616 // access only rcvDomains
619 MInt err = MPI_Irecv(&zBC->rcvComm[noRcv]->buffer[0], zBC->rcvComm[noRcv]->bufferSize * zBC->noZonalVariables,
620 MPI_DOUBLE, zBC->rcvComm[noRcv]->localId, 0, m_StructuredComm, &zBC->rcvRequest[noRcv],
621 AT_, "&zBC->bufferRcvZonal[noRcv][0]");
623 }
624 }
625 }
626 MPI_Waitall(zBC->noSndNghbrDomains, zBC->sndRequest, zBC->sndStatus, AT_);
627 MPI_Waitall(zBC->noRcvNghbrDomains, zBC->rcvRequest, zBC->rcvStatus, AT_);
628 }
629}
◆ zonalScatter()
| void FvStructuredSolver3D::zonalScatter | ( | ) |
Definition at line 632 of file fvstructuredsolver3d.cpp.
632 {
638 if(zBC->hasSTG) {
645 } else {
652 }
653 }
654 }
655 }
656}
◆ zonalSend()
| void FvStructuredSolver3D::zonalSend | ( | ) |
Definition at line 597 of file fvstructuredsolver3d.cpp.
597 {
599 // access only sndDomains
603 MInt err = MPI_Isend(&zBC->sndComm[noSnd]->buffer[0], zBC->sndComm[noSnd]->bufferSize * zBC->noZonalVariables,
604 MPI_DOUBLE, zBC->sndComm[noSnd]->localId, 0, m_StructuredComm, &zBC->sndRequest[noSnd],
605 AT_, "zBC->buffer[noSnd][0]");
607 }
608 }
609 }
610 }
611}
Friends And Related Function Documentation
◆ FvStructuredSolver3DRans
|
friend |
Definition at line 35 of file fvstructuredsolver3d.h.
◆ StructuredBndryCnd3D
template<MBool isRans>
|
friend |
Definition at line 34 of file fvstructuredsolver3d.h.
Member Data Documentation
◆ m_kineticEOld
|
protected |
Definition at line 190 of file fvstructuredsolver3d.h.
◆ m_pointInterpolation
|
protected |
Definition at line 208 of file fvstructuredsolver3d.h.
◆ m_ransSolver
|
protected |
Definition at line 203 of file fvstructuredsolver3d.h.
◆ m_structuredBndryCnd
|
protected |
Definition at line 183 of file fvstructuredsolver3d.h.
◆ m_structuredInterpolation
|
protected |
Definition at line 207 of file fvstructuredsolver3d.h.
◆ m_zonalBC
| std::vector<std::unique_ptr<StructuredZonalBC> > FvStructuredSolver3D::m_zonalBC |
Definition at line 69 of file fvstructuredsolver3d.h.
◆ nDim
|
staticconstexprprotected |
Definition at line 227 of file fvstructuredsolver3d.h.
◆ reconstructSurfaceData
| void(FvStructuredSolver3D::* FvStructuredSolver3D::reconstructSurfaceData) () |
Definition at line 103 of file fvstructuredsolver3d.h.
◆ Venkatakrishnan_function
| void(FvStructuredSolver3D::* FvStructuredSolver3D::Venkatakrishnan_function) (MFloat, MFloat, MFloat, MInt, MInt, MFloatScratchSpace &) |
Definition at line 119 of file fvstructuredsolver3d.h.
◆ viscFluxMethod
| void(FvStructuredSolver3D::* FvStructuredSolver3D::viscFluxMethod) () |
Definition at line 129 of file fvstructuredsolver3d.h.
◆ xsd
|
staticconstexprprotected |
Definition at line 186 of file fvstructuredsolver3d.h.
◆ ysd
|
staticconstexprprotected |
Definition at line 187 of file fvstructuredsolver3d.h.
◆ zsd
|
staticconstexprprotected |
Definition at line 188 of file fvstructuredsolver3d.h.
The documentation for this class was generated from the following files:
- /home/gitlab-runner/scratch/builds/oxpnswJ6/1/aia/m-AIA/m-AIA/src/FV/fvstructuredsolver3d.h
- /home/gitlab-runner/scratch/builds/oxpnswJ6/1/aia/m-AIA/m-AIA/src/FV/fvstructuredsolver3d.cpp
