Loading...
Searching...
No Matches
Classes |
Public Member Functions |
Private Types |
Private Member Functions |
Private Attributes |
Static Private Attributes |
List of all members
AcaSolver< nDim > Class Template Reference
Acoustic extrapolation of near/near-far field data to the true far-field using an integral method (currently only FWH, but should allow for other methods to be implemented as well) More...
#include <acasolver.h>
Inheritance diagram for AcaSolver< nDim >:
Collaboration diagram for AcaSolver< nDim >:
Classes | |
| struct | AcaSymBc |
| Hold data for symmetry boundary condition. More... | |
| struct | ConversionFactors |
| Conversion-Factors for the input conversion. More... | |
| struct | FOURIER |
| struct | FWH |
| struct | FWH_APE |
| struct | FWH_TIME |
| struct | PP |
| struct | WINDOW |
Public Member Functions | |
| AcaSolver (const MInt solverId, const MPI_Comm comm) | |
| ~AcaSolver () override | |
| MBool | solutionStep () override |
| Main compute method to perform acoustic far-field prediction. More... | |
| MBool | solverConverged () |
| MBool | isActive () const override |
| void | initSolver () override |
| Initialize solver. More... | |
| void | finalizeInitSolver () override |
| MFloat | time () const override |
| Return the time. More... | |
| MInt | noVariables () const override |
| Return the number of variables. More... | |
| MInt | noInternalCells () const override |
| Return the number of internal cells within this solver. More... | |
| void | preTimeStep () override |
| void | postTimeStep () |
| virtual void | saveSolverSolution (const MBool, const MBool) |
| virtual void | cleanUp () |
 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 () |
Private Types | |
| enum | NONDIMBASIS { NONDIMACA = 0 , NONDIMSTAGNATION = 1 , NONDIMLB = 2 } |
| using | SurfaceDataCollector = maia::acoustic_analogy::collector::SurfaceDataCollector< nDim > |
| using | ObserverDataCollector = maia::acoustic_analogy::observer_collector::ObserverDataCollector< nDim > |
| using | SourceParameters = maia::acoustic::SourceParameters |
| using | SourceVars = maia::acoustic::SourceVars |
| using | Timers = maia::acoustic::Timers_ |
Private Member Functions | |
| void | initTimers () |
| Initialize solver timers. More... | |
| void | averageTimer () |
| Average all FWH timer over participating ranks. More... | |
| void | setInputOutputProperties () |
| Read/set all I/O related properties. More... | |
| void | setNumericalProperties () |
| Read/set all numerical properties. More... | |
| MInt | noSurfaces () |
| Return number of surfaces. More... | |
| MBool | hasSurface (const MInt surfaceId) |
| Return if this rank has data of this surface. More... | |
| MInt | noSurfaceElements (const MInt surfaceId) |
| Return number of (local) surfaces elements for a surface. More... | |
| MInt | totalNoSurfaceElements () |
| Return total number of (local) surfaces elements. More... | |
| MInt | localSurfaceOffset (const MInt sId) |
| Return local offset of a surface in m_surfaceData. More... | |
| MInt | noSamples () |
| Return number of samples (i.e. number of time steps) More... | |
| MInt | noObservers () |
| Return local number of observer points. More... | |
| MInt | observerOffset () |
| Return local number of observer points. More... | |
| MInt | globalNoObservers () |
| Return number of observer points. More... | |
| MInt | a_localObserverId (const MInt globalObserverId) |
| Return local observer id. More... | |
| MInt | a_globalObserverId (const MInt localObserverId) |
| Return global observer id. More... | |
| MFloat | a_globalObserverCoordinate (const MInt globalObserverId, const MInt dir) |
| Accessor for global observer coordinates. More... | |
| std::array< MFloat, nDim > | getGlobalObserverCoordinates (const MInt globalObserverId) |
| Accessor for global observer coordinates. More... | |
| MInt | getObserverDomainId (const MInt globalObserverId) |
| Get domain id of (global) observer. More... | |
| void | initObserverPoints () |
| Initialize observer points. More... | |
| void | initParallelization () |
| Initialize parallelization (distribution of surface segments on all domains) More... | |
| void | distributeObservers () |
| Distribute observers on all ranks. More... | |
| MString | getSurfaceDataFileName (const MInt surfaceId, const MInt fileNo) |
| IO methods. More... | |
| void | readSurfaceData () |
| Read the surface data (or generate for an analytical case) More... | |
| void | readSurfaceDataFromFile (const MInt surfaceId) |
| Read the data for one surface. More... | |
| void | storeSurfaceData (const MInt surfaceId) |
| Store the surface data. Note: used for validating the data input and for the output of generated data. More... | |
| void | readNoSegmentsAndSamples (const MInt surfaceId, MInt &noSegments, MInt &noSamples, MString &inputFileName) |
| Read the number of segments and samples for the given surface id. More... | |
| void | printMessage (const MString msg) |
| Info functions. More... | |
| void | calcSourceTerms () |
| Main compute functions. More... | |
| void | calcWindowFactors (MFloat *const p_window) |
| Compute the window function factors and the normalization factor. More... | |
| void | windowAndTransformSources () |
| Apply window function to source terms and transform into frequency domain. More... | |
| void | windowAndTransformObservers () |
| Apply window function to the observer time series data and transform it into frequency domain. More... | |
| void | calcFrequencies () |
| Calculate the frequencies. More... | |
| void | calcSurfaceIntegralsForObserver (const std::array< MFloat, nDim > coord, MFloat *const p_complexVariables) |
| Calculate the FWH surface integrals for a certain observer location. More... | |
| void | calcSurfaceIntegralsAtSourceTime (const MInt globalObserverId) |
| Calculate the FWH surface integrals for a certain observer using the time-domain formulation. More... | |
| MFloat | calcTimeDerivative (const MInt segmentId, const MInt varId, const MInt tau) |
| void | interpolateFromSourceToObserverTime (const MFloat distMinObservers, MFloat *const p_perturbedFWH) |
| Interpolate the observer pressure signal from source time to observer time using the advanced time approach by Casalino. Reference: Casalino, D., 2003, An advanced time approach ... doi.org/10.1016/s0022-460x(02)00986-0. More... | |
| void | combineSurfaceIntegrals (const MInt globalObserverId, MFloat *const p_complexVariables) |
| Combine the FWH surface integrals of all ranks. More... | |
| void | combineSurfaceIntegralsInTime (const MInt globalObserverId, MFloat *const p_perturbedFWH) |
| Combine the FWH surface integrals of the time domain method of all ranks. More... | |
| void | calcMinDistanceForObservers (MFloat *const distMinObservers) |
| Find the minimum distance from the surface elements to the observers. More... | |
| void | calculateObserverInFrequency () |
| Calculate the observer signals in frequency domain. More... | |
| void | calculateObserverInTime () |
| Calculate the observer signals in time domain. More... | |
| void | backtransformObservers () |
| Backtransform the observer frequency signals into the time domain. More... | |
| void | saveObserverSignals () |
| Save the observer signals. Save observer data: 1. pressure data in time domain and 2. pressure data in frequency domain. More... | |
| void | loadObserverSignals () |
| Load the observer signals from a file for further postprocessing. More... | |
| void | postprocessObserverSignals () |
| Postprocessing of observer signals. More... | |
| void | FastFourierTransform (MFloat *dataArray, const MInt size, const MInt dir, MFloat *result, const MBool inPlace) |
| Functions for Fourier Transformation. More... | |
| void | DiscreteFourierTransform (MFloat *dataArray, const MInt size, const MInt dir, MFloat *result, const MBool inPlace) |
| Discrete Fourier transform (computational expensive O(n^2) instead of O(n*logn) compared to FFT) More... | |
| void | checkNoSamplesPotencyOfTwo () |
| void | calcSourceTermsFwhFreq (const MInt segmentId) |
| Source term compute kernels. More... | |
| void | calcSourceTermsFwhTime (const MInt segmentId) |
| Calculate FWH source terms for the time domain FWH formulation Reference: [1] Brentner, K. S. and F. Farassat, 1998, Analytical Comparison ... doi.org/10.2514/2.558 [2] Casalino, D., 2003, An advanced time approach ... doi.org/10.1016/s0022-460x(02)00986-0 [3] Brès, G., 2010, A Ffowcs Williams - Hawkings ... doi.org/10.2514/6.2010-3711. More... | |
| void | transformCoordinatesToFlowDirection2D (MFloat *const fX, MFloat *const fY) |
| Coordinate transformation. More... | |
| void | transformCoordinatesToFlowDirection3D (MFloat *const fX, MFloat *const fY, MFloat *const fZ) |
| void | transformBackToOriginalCoordinate2D (MFloat *const fX, MFloat *const fY) |
| void | transformBackToOriginalCoordinate3D (MFloat *const fX, MFloat *const fY, MFloat *const fZ) |
| void | generateSurfaces () |
| Generation of surfaces. More... | |
| void | generateSurfacePlane (const MInt sId) |
| Generate a surface plane (Cartesian). More... | |
| void | generateSurfaceCircle (const MInt sId) |
| Generate a circular surface in 2D. More... | |
| MInt | readNoSegmentsSurfaceGeneration (const MInt sId, std::vector< MInt > &noSegments, const MInt lengthCheck=-1) |
| Determine the number of segments to be generated for a surface, returns the total number of segments and a vector with the number of segments in each direction for the specific surface type. More... | |
| void | generateSurfaceData () |
| Generation of analytical input data. More... | |
| void | genMonopoleAnalytic2D (const MFloat *coord, const SourceParameters ¶m, SourceVars &vars) |
| Monopole. More... | |
| void | genMonopoleAnalytic3D (const MFloat *coord, const SourceParameters ¶m, SourceVars &vars) |
| 3D analytic monopole generation More... | |
| void | genDipoleAnalytic2D (const MFloat *coord, const SourceParameters ¶m, SourceVars &vars) |
| Dipole. More... | |
| void | genDipoleAnalytic3D (const MFloat *coord, const SourceParameters ¶m, SourceVars &vars) |
| 3D dipole analytic More... | |
| void | genQuadrupoleAnalytic2D (const MFloat *coord, const SourceParameters ¶m, SourceVars &vars) |
| Quadrupole. More... | |
| void | genQuadrupoleAnalytic3D (const MFloat *coord, const SourceParameters ¶m, SourceVars &vars) |
| 3D quadrupole analytic More... | |
| void | genVortexConvectionAnalytic2D (const MFloat *obsCoord, const SourceParameters &m_sourceParameters, SourceVars &sourceVars) |
| Inviscid vortex convection. More... | |
| void | applySymmetryBc (const std::array< MFloat, nDim > coord, MFloat *const p_complexVariables) |
| Symmetry boundary condition. More... | |
| void | calcObsPressureAnalytic () |
| void | initConversionFactors () |
| Initialize the conversion factors required to convert between. More... | |
| void | changeDimensionsSurfaceData () |
| Change of the non-dimensionalization of input data from the stagnation state (0) non-dim. to the undisturbed freestream state (inf) non-dim. More... | |
| void | changeDimensionsObserverData () |
| void | computeDt () |
| void | computeAccuracy () |
Private Attributes | |
| std::vector< std::unique_ptr< AcaPostProcessing > > | m_post |
| MBool | m_isFinished = false |
| MInt | m_acaResultsNondimMode = 0 |
| struct AcaSolver::ConversionFactors | m_input2Aca |
| struct AcaSolver::ConversionFactors | m_aca2Output |
| MInt | m_noSurfaces = -1 |
| Number of surfaces. More... | |
| std::vector< MInt > | m_surfaceIds {} |
| Ids of the surface input file. More... | |
| MBool | m_useMergedInputFile |
| MInt | m_inputFileIndexStart |
| MInt | m_inputFileIndexEnd |
| std::vector< MInt > | m_noSurfaceElements {} |
| Number of (local) surface elements for each surface. More... | |
| MBool | m_allowMultipleSurfacesPerRank = false |
| Allow multiple surfaces per rank. More... | |
| std::vector< MInt > | m_noSurfaceElementsGlobal {} |
| Total number of surface elements for each surface. More... | |
| std::vector< MInt > | m_surfaceElementOffsets {} |
| Local surface element offsets for each surface. More... | |
| std::vector< MPI_Comm > | m_mpiCommSurface {} |
| Surface local MPI communicators. More... | |
| std::vector< MInt > | m_noDomainsSurface {} |
| Surface local number of domains. More... | |
| std::vector< MInt > | m_domainIdSurface {} |
| Surface local domain id. More... | |
| std::vector< MFloat > | m_surfaceWeightingFactor {} |
| Weighting factor for each surface (for surface averaging) More... | |
| std::vector< MString > | m_surfaceInputFileName {} |
| Surface input file names that were used for sampling. More... | |
| MInt | m_noSamples = -1 |
| Number of samples. More... | |
| MInt | m_totalNoSamples = -1 |
| Total number of samples in input data files. More... | |
| MInt | m_sampleStride = 1 |
| Stride of used samples in input data files. More... | |
| MInt | m_sampleOffset = 0 |
| Offset in input data files from which to start reading samples. More... | |
| SurfaceDataCollector | m_surfaceData |
| Collector for surface data. More... | |
| MInt | m_acousticMethod = -1 |
| Acoustic extrapolation method to use, see enum ExtrapolationMethods in enums.h. More... | |
| MInt | m_fwhTimeShiftType = 0 |
| MInt | m_fwhMeanStateType = 0 |
| MInt | m_noVariables = -1 |
| MInt | m_noComplexVariables = -1 |
| std::map< MString, MInt > | m_inputVars {} |
| Set of input variables. More... | |
| MInt | m_windowType = 0 |
| Window type. More... | |
| MFloat | m_windowNormFactor = -1.0 |
| Normalization factor of window. More... | |
| MInt | m_transformationType = 0 |
| Fourier type. More... | |
| MInt | m_noGlobalObservers = -1 |
| Number of samples. More... | |
| std::vector< MFloat > | m_globalObserverCoordinates |
| ObserverDataCollector | m_observerData |
| Coordinates of (global) observer points. More... | |
| MInt | m_offsetObserver = 0 |
| Collector for (local) observer data. More... | |
| MInt | m_noObservers = 0 |
| Offset of local observer. More... | |
| std::array< MFloat, nDim *nDim > | m_transformationMatrix |
| Number of local observer. More... | |
| MBool | m_zeroFlowAngle = false |
| MString | m_inputDir |
| I/O. More... | |
| MString | m_outputFilePrefix = "" |
| Input directory. More... | |
| MBool | m_storeSurfaceData = false |
| Prefix for all output files. More... | |
| MString | m_observerFileName |
| File containing observer points. More... | |
| MBool | m_generateObservers = false |
| Enable observer point generation (instead of reading from file) More... | |
| MBool | m_generateSurfaceData = false |
| Enable generation of analytical surface data. More... | |
| MInt | m_sourceType = -1 |
| Source Type. More... | |
| maia::acoustic::SourceParameters | m_sourceParameters |
| Parameter vector. More... | |
| MFloat | m_omega_factor = -1.0 |
| Omega factor. More... | |
| MFloat | m_lambdaZero = -1.0 |
| Wave length lambda. More... | |
| MFloat | m_lambdaMach = -1.0 |
| MFloat | m_cInfty = -1.0 |
| Non-dimensionalization variables. More... | |
| MFloat | m_rhoInfty = -1.0 |
| std::vector< MFloat > | m_Analyticfreq {} |
| Analytic Fourier Transformation, also used for Validation. More... | |
| std::vector< MFloat > | m_rmsP_Analytic {} |
| Analytic rms pressure for accuracy function. More... | |
| MFloat | m_dt = -1.0 |
| Constant time step size. More... | |
| MFloat | m_MaVec [3] |
| Mach vector. More... | |
| MFloat | m_Ma = -1.0 |
| Mach number. More... | |
| MFloat | m_MaDim = -1.0 |
| Mach number used for change of nondimensionalization. More... | |
| MBool | m_FastFourier = false |
| Can FFT be used? More... | |
| std::vector< MFloat > | m_times {} |
| Storage for times and frequencies. More... | |
| std::vector< MFloat > | m_frequencies {} |
| MBool | m_hasTimes = false |
| MInt | m_noPostprocessingOps = 0 |
| Post Processing of observer signals. More... | |
| std::vector< MInt > | m_postprocessingOps {} |
| MBool | m_acaPostprocessingMode = false |
| MString | m_acaPostprocessingFile = "" |
| MInt | m_timerGroup = -1 |
| std::array< MInt, Timers::_count > | m_timers {} |
| std::vector< AcaSymBc > | m_symBc |
Static Private Attributes | |
| static constexpr MFloat | m_gamma = 1.4 |
Additional Inherited Members | |
| Public Attributes inherited from Solver | |
| std::set< MInt > | m_freeIndices |
| MBool | m_singleAdaptation = false |
| MBool | m_splitAdaptation = true |
| MBool | m_saveSensorData = false |
| 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 Attributes inherited from Solver | |
| MFloat | m_Re {} |
| the Reynolds number More... | |
| MFloat | m_Ma {} |
| the Mach number More... | |
| MInt | m_solutionInterval |
| The number of timesteps before writing the next solution file. More... | |
| MInt | m_solutionOffset {} |
| std::set< MInt > | m_solutionTimeSteps |
| MInt | m_restartInterval |
| The number of timesteps before writing the next restart file. More... | |
| MInt | m_restartTimeStep |
| MInt | m_restartOffset |
| MString | m_solutionOutput |
| MBool | m_useNonSpecifiedRestartFile = false |
| MBool | m_initFromRestartFile |
| MInt | m_residualInterval |
| The number of timesteps before writing the next residual. More... | |
| const MInt | m_solverId |
| a unique solver identifier More... | |
| MFloat * | m_outerBandWidth = nullptr |
| MFloat * | m_innerBandWidth = nullptr |
| MInt * | m_bandWidth = nullptr |
| MBool | m_restart = false |
| MBool | m_restartFile = false |
Detailed Description
template<MInt nDim>
class AcaSolver< nDim >
class AcaSolver< nDim >
Reference for implementation/validation of FWH solver: 'Predition of far-field noise using the Ffowcs Williams-Hawkings method', Maximilian Kerstan, Master thesis, 2020.
Definition at line 78 of file acasolver.h.
Member Typedef Documentation
◆ ObserverDataCollector
template<MInt nDim>
|
private |
Definition at line 82 of file acasolver.h.
◆ SourceParameters
template<MInt nDim>
|
private |
Definition at line 83 of file acasolver.h.
◆ SourceVars
template<MInt nDim>
|
private |
Definition at line 84 of file acasolver.h.
◆ SurfaceDataCollector
template<MInt nDim>
|
private |
Definition at line 81 of file acasolver.h.
◆ Timers
template<MInt nDim>
|
private |
Definition at line 85 of file acasolver.h.
Member Enumeration Documentation
◆ NONDIMBASIS
template<MInt nDim>
|
private |
| Enumerator | |
|---|---|
| NONDIMACA | |
| NONDIMSTAGNATION | |
| NONDIMLB | |
Definition at line 87 of file acasolver.h.
Constructor & Destructor Documentation
◆ AcaSolver()
template<MInt nDim>
| AcaSolver< nDim >::AcaSolver | ( | const MInt | solverId, |
| const MPI_Comm | comm | ||
| ) |
Definition at line 24 of file acasolver.cpp.
Parent class of all solvers This class is the base for all solvers. I.e. all solver class (e....
Definition: solver.h:29
◆ ~AcaSolver()
Definition at line 96 of file acasolver.h.
96 {
97 averageTimer();
99 };
Member Function Documentation
◆ a_globalObserverCoordinate()
template<MInt nDim>
|
inlineprivate |
Definition at line 165 of file acasolver.h.
◆ a_globalObserverId()
template<MInt nDim>
|
inlineprivate |
Definition at line 163 of file acasolver.h.
◆ a_localObserverId()
template<MInt nDim>
|
inlineprivate |
Definition at line 155 of file acasolver.h.
155 {
156 if(globalObserverId < m_offsetObserver || (m_offsetObserver + m_noObservers) < globalObserverId) {
157 return -1;
158 } else {
160 }
161 };
◆ applySymmetryBc()
template<MInt nDim>
|
private |
Apply symmetry boundary condition.
- Date
- 15.06.2023
- Parameters
-
[in] coord Location of the observer [out] p_complexVariables Output of complex p', length of 2*noSamples()
Calculates surface integral for mirrored observer points. The resulting complex pressure is then added to the image observer points.
Definition at line 4547 of file acasolver.cpp.
4547 {
4548 TRACE();
4550 if(noSym == 0) return;
4551 TERMM_IF_COND(noSym > 1, "WARNING: only implemented for one symmetry BC so far");
4552 auto transform = [](std::array<MFloat, nDim>& coord_, const AcaSymBc& bc) {
4553 const MFloat coordTimesNormal = std::inner_product(&coord_[0], &coord_[nDim], &bc.normal[0], .0);
4555 coord_[d] += -2 * coordTimesNormal * bc.normal[d] + 2 * bc.origin[d];
4556 }
4557 };
4559 // transform observer position
4560 std::array<MFloat, nDim> symCoord = coord;
4561 transform(symCoord, m_symBc[0]);
4562 // Calculate pressure for transformed point before adding to its 'image' observer
4563 calcSurfaceIntegralsForObserver(symCoord, complexVariables_.data());
4565 p_complexVariables[2 * i + 0] += complexVariables_[2 * i + 0];
4566 p_complexVariables[2 * i + 1] += complexVariables_[2 * i + 1];
4567 }
4568}
void calcSurfaceIntegralsForObserver(const std::array< MFloat, nDim > coord, MFloat *const p_complexVariables)
Calculate the FWH surface integrals for a certain observer location.
Definition: acasolver.cpp:1573
◆ averageTimer()
- Date
- 21.02.2024
Definition at line 74 of file acasolver.cpp.
74 {
75 // 0) map timer ids for safety
76 std::vector<MInt> timerIds_;
77 timerIds_.reserve(11);
90 // 1) fill buffer with local timer values
91 std::vector<MFloat> timerValues_;
92 timerValues_.reserve(noTimers);
94 timerValues_.emplace_back(RETURN_TIMER_TIME(timerIds_[i]));
95 }
96 // 2) collect values from all ranks and use max value
97 MPI_Allreduce(MPI_IN_PLACE, timerValues_.data(), noTimers, maia::type_traits<MFloat>::mpiType(), MPI_MAX, mpiComm(),
98 AT_, "MPI_IN_PLACE", "timerValues_");
100 SET_RECORD(timerIds_[i], timerValues_[i]);
101 }
102}
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
@ WindowAndTransformSources
Definition: acasolver.h:36
Definition: typetraits.h:22
◆ backtransformObservers()
Definition at line 2348 of file acasolver.cpp.
2348 {
2349 TRACE();
2352
2353 // Data buffer for backtransformed pressure data
2355
2356 // BACKTRANSFORM
2358 // Copy the first calculated half to the other half starting from behind going to the middle
2359 // (since it is a two-sided frequency plot)
2361 m_observerData.complexVariables(obs, 0, noSamples() - i, 0) = m_observerData.complexVariables(obs, 0, i, 0);
2364 }
2365
2366 std::fill_n(&backtransform[0], 2 * noSamples(), 0.0);
2367 // Backtransformation
2369 FastFourierTransform(&m_observerData.complexVariables(obs, 0), noSamples(), -1, &backtransform[0], false);
2371 DiscreteFourierTransform(&m_observerData.complexVariables(obs, 0), noSamples(), -1, &backtransform[0], false);
2372 }
2373
2374 // Storing the results in observer.variables since it is back in time domain (only real part)
2377 }
2378 }
2379
2381}
ObserverDataCollector m_observerData
Coordinates of (global) observer points.
Definition: acasolver.h:369
void DiscreteFourierTransform(MFloat *dataArray, const MInt size, const MInt dir, MFloat *result, const MBool inPlace)
Discrete Fourier transform (computational expensive O(n^2) instead of O(n*logn) compared to FFT)
Definition: acasolver.cpp:1500
void FastFourierTransform(MFloat *dataArray, const MInt size, const MInt dir, MFloat *result, const MBool inPlace)
Functions for Fourier Transformation.
Definition: acasolver.cpp:1411
MFloat & complexVariables(const MInt id, const MInt var)
Accessor for complex variables.
Definition: acaobserverdatacollector.h:302
MFloat & variables(const MInt id, const MInt var)
Accessor for variables (return a reference to the start of the given variable in memory).
Definition: acaobserverdatacollector.h:261
◆ calcFrequencies()
Definition at line 1549 of file acasolver.cpp.
1549 {
1550 TRACE();
1552 // Calculate frequencies, m_frequencies is already resized with noSamples
1553 // Numerical Recipes in C by Cambridge page 503
1554 // Structure: 0 +1 +2 +3 ... +noSamples/2 ... -3 -2 -1
1558 }
1560 }
1561 // Frequencies from Hz --> rad/s: now omega
1563 m_frequencies[i] *= 2.0 * PI;
1564 }
1565}
◆ calcMinDistanceForObservers()
template<MInt nDim>
|
private |
Definition at line 2069 of file acasolver.cpp.
2069 {
2071
2072 // Calculate Mach number and the Prandtl-Glauert-Factor
2075
2076 printMessage("- Loop over observers to find minimum distance from the source to the observers");
2078 MFloat distMinLocal = std::numeric_limits<MFloat>::max();
2079
2080 // Get observer coordinates
2082
2083 // SURFACE LOOP
2086
2087 if constexpr(nDim == 3) {
2088 // Geometric Distance & Coordinate Transformation
2089 MFloat distX = coord[0] - surfCoordinates[0];
2090 MFloat distY = coord[1] - surfCoordinates[1];
2091 MFloat distZ = coord[2] - surfCoordinates[2];
2092 transformCoordinatesToFlowDirection3D(&distX, &distY, &distZ);
2093
2094 // Effective Acoustic Distance
2097
2098 // Find the minimum distance of all segments
2100 } else if constexpr(nDim == 2) {
2101 TERMM(1, "Time Domain FW-H not implimented in 2D: Requires free space Green function in 2D");
2102 } // End of the dimension switch
2103 } // End of surface loop
2104
2105 distMinObservers[k] = distMinLocal;
2106 } // End of observer loop
2107
2108 MPI_Allreduce(MPI_IN_PLACE, &distMinObservers[0], globalNoObservers(), type_traits<MFloat>::mpiType(), MPI_MIN,
2110
2112 const MFloat disMinGlobal = *std::min_element(distMinObservers, distMinObservers + globalNoObservers());
2113 std::fill(distMinObservers, distMinObservers + globalNoObservers(), disMinGlobal);
2115 std::cout << "Apply universal time shift for all observers = " << disMinGlobal << std::endl;
2116 }
2117 }
2118}
void transformCoordinatesToFlowDirection3D(MFloat *const fX, MFloat *const fY, MFloat *const fZ)
Definition: acasolver.cpp:2006
std::array< MFloat, nDim > getGlobalObserverCoordinates(const MInt globalObserverId)
Accessor for global observer coordinates.
Definition: acasolver.h:169
MInt totalNoSurfaceElements()
Return total number of (local) surfaces elements.
Definition: acasolver.h:138
MFloat & surfaceCoordinates(const MInt id)
Accessor for surface element coordinates.
Definition: acasurfacedatacollector.h:259
◆ calcObsPressureAnalytic()
Definition at line 2583 of file acasolver.cpp.
2583 {
2584 TRACE();
2585
2588 vars.fill(0.0);
2589
2591 // Pointers to pass to functions
2593 SourceVars sourceVars;
2596 sourceVars.u = &vars(FWH::U, obsId, 0);
2597 sourceVars.v = &vars(FWH::V, obsId, 0);
2598 sourceVars.w = &vars(FWH::W, obsId, 0);
2599 sourceVars.rho = &vars(FWH::RHO, obsId, 0);
2600 } else {
2601 sourceVars.p = &vars(FWH_APE::PP, obsId, 0);
2602 sourceVars.u = &vars(FWH_APE::UP, obsId, 0);
2603 sourceVars.v = &vars(FWH_APE::VP, obsId, 0);
2604 sourceVars.w = &vars(FWH_APE::WP, obsId, 0);
2605 }
2606
2607 if constexpr(nDim == 2) {
2609 case 0: {
2611 break;
2612 }
2613 case 1: {
2615 break;
2616 }
2617 case 2: {
2619 break;
2620 }
2621 case 3: {
2623 break;
2624 }
2625 default: {
2626 TERMM(1, "source type not implemented");
2627 break;
2628 }
2629 }
2630 } else {
2632 case 0: {
2634 break;
2635 }
2636 case 1: {
2638 break;
2639 }
2640 case 2: {
2642 break;
2643 }
2644 default: {
2645 TERMM(1, "source type not implemented");
2646 break;
2647 }
2648 }
2649 }
2650 }
2651
2652 // OUTPUT
2654 const MString fileName = outputDir() + m_outputFilePrefix + "obsPressAnalytic" + ParallelIo::fileExt();
2655
2657
2659 file.defineArray(PIO_FLOAT, "obsPressAnalytic", 2, &dimSizes[0]);
2660
2664 } else {
2666 }
2667}
void genQuadrupoleAnalytic2D(const MFloat *coord, const SourceParameters ¶m, SourceVars &vars)
Quadrupole.
Definition: acasolver.cpp:4307
void genDipoleAnalytic2D(const MFloat *coord, const SourceParameters ¶m, SourceVars &vars)
Dipole.
Definition: acasolver.cpp:4160
MInt m_acousticMethod
Acoustic extrapolation method to use, see enum ExtrapolationMethods in enums.h.
Definition: acasolver.h:347
void genQuadrupoleAnalytic3D(const MFloat *coord, const SourceParameters ¶m, SourceVars &vars)
3D quadrupole analytic
Definition: acasolver.cpp:4409
void genDipoleAnalytic3D(const MFloat *coord, const SourceParameters ¶m, SourceVars &vars)
3D dipole analytic
Definition: acasolver.cpp:4233
void genMonopoleAnalytic3D(const MFloat *coord, const SourceParameters ¶m, SourceVars &vars)
3D analytic monopole generation
Definition: acasolver.cpp:4097
void genMonopoleAnalytic2D(const MFloat *coord, const SourceParameters ¶m, SourceVars &vars)
Monopole.
Definition: acasolver.cpp:4035
void genVortexConvectionAnalytic2D(const MFloat *obsCoord, const SourceParameters &m_sourceParameters, SourceVars &sourceVars)
Inviscid vortex convection.
Definition: acasolver.cpp:4506
MFloat & observerCoordinates(const MInt id)
Accessor for surface element coordinates.
Definition: acaobserverdatacollector.h:229
Definition: parallelio.h:26
◆ calcSourceTerms()
Compute source terms for all surface elements depending on the selected method.
Definition at line 879 of file acasolver.cpp.
879 {
880 TRACE();
883
884 // Note: totalNoSurfaceElements() returns the local, total number of elements for a surface
885 // (may differ from domain to domain)
886
890 calcSourceTermsFwhFreq(i);
891 }
895 calcSourceTermsFwhTime(i);
896 }
897 } else {
898 TERMM(1, "Unknown solverMethod");
899 }
900
902}
void calcSourceTermsFwhTime(const MInt segmentId)
Calculate FWH source terms for the time domain FWH formulation Reference: [1] Brentner,...
Definition: acasolver.cpp:1062
void calcSourceTermsFwhFreq(const MInt segmentId)
Source term compute kernels.
Definition: acasolver.cpp:909
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
Definition: enums.cpp:20
◆ calcSourceTermsFwhFreq()
Calculate FWH source terms for the frequency FWH formulation Reference: [1] Lockard, D.P., 2000. An efficient, two-dimensional ... doi.org/10.1006/jsvi.1999.2522 [2] Lockard, D.P., 2002, A comparison of Ffowcs ... doi.org/10.2514/6.2002-2580.
Definition at line 909 of file acasolver.cpp.
909 {
910 TRACE();
911
912 // Get pointer to normal vector
914 // First two normal vector components (2D)
917 // Mach numbers (2D)
920
925
926 if constexpr(nDim == 3) {
927 // Third normal vector component for 3D
929 // Third Mach number for 3D
931
932 // Compute mean velocities
933 MFloat uMean = 0.0;
934 MFloat vMean = 0.0;
935 MFloat wMean = 0.0;
938 case 0: { // Default: use infnity properties to compute u'
939 uMean = m_MaVec[0];
940 vMean = m_MaVec[1];
941 wMean = m_MaVec[2];
942 break;
943 }
944 case 1: { // Use mean velocities to compute u'
949 }
950 uMean *= noSamplesInv;
951 vMean *= noSamplesInv;
952 wMean *= noSamplesInv;
953 break;
954 }
955 default: {
956 TERMM(1, "fwh mean state type not implemented");
957 break;
958 }
959 }
960 }
961
962 // Time loop over all samples for 3D
964 // Get required variables (u', v', w', p, rho)
965 std::vector<MInt> VarDim;
966 MFloat up = 0.0;
967 MFloat vp = 0.0;
968 MFloat wp = 0.0;
970 MFloat rho = 0.0;
985 }
986
989
990 // Calculation of source terms (3D) Fx, Fy, Fz and Q for each sample.
995 });
996 } else {
997 // Compute mean velocities
998 MFloat uMean = 0.0;
999 MFloat vMean = 0.0;
1002 case 0: { // Default: use infnity properties to compute u'
1003 uMean = m_MaVec[0];
1004 vMean = m_MaVec[1];
1005 break;
1006 }
1007 case 1: { // Use mean velocities to compute u'
1011 }
1012 uMean = uMean * noSamplesInv;
1013 vMean = vMean * noSamplesInv;
1014 break;
1015 }
1016 default: {
1017 TERMM(1, "fwh mean state type not implemented");
1018 break;
1019 }
1020 }
1021 }
1022
1023 // Time loop over all samples for 2D
1025 // Get required variables (u', v', p, rho)
1026 std::vector<MInt> VarDim;
1027 MFloat up = 0.0;
1028 MFloat vp = 0.0;
1030 MFloat rho = 0.0;
1043 }
1044
1047
1048 // Calculation of source terms (2D) Fx, Fy and Q for each sample.
1049 m_surfaceData.variables(segmentId, VarDim[0], t) = p * dfdx + rho * (up - Max) * f1 + Max * f2;
1050 m_surfaceData.variables(segmentId, VarDim[1], t) = p * dfdy + rho * (vp - May) * f1 + May * f2;
1052 });
1053 }
1054}
MFloat & variables(const MInt id, const MInt var)
Accessor for variables.
Definition: acasurfacedatacollector.h:325
MFloat & surfaceNormal(const MInt id)
Accessor for surface normal.
Definition: acasurfacedatacollector.h:303
constexpr std::underlying_type< FcCell >::type p(const FcCell property)
Converts property name to underlying integer value.
Definition: fccellproperties.h:32
◆ calcSourceTermsFwhTime()
Definition at line 1062 of file acasolver.cpp.
1062 {
1063 TRACE();
1064
1065 // Get pointer to normal vector
1071
1072 if constexpr(nDim == 3) {
1073 // Compute mean variables
1074 MFloat uMean = 0.0;
1075 MFloat vMean = 0.0;
1076 MFloat wMean = 0.0;
1077 MFloat pMean = 0.0;
1080 case 0: { // Default: use infnity properties to compute u'
1081 uMean = m_MaVec[0];
1082 vMean = m_MaVec[1];
1083 wMean = m_MaVec[2];
1084 pMean = p_inf;
1085 break;
1086 }
1087 case 1: { // Use mean velocities to compute u'
1093 }
1094 uMean *= noSamplesInv;
1095 vMean *= noSamplesInv;
1096 wMean *= noSamplesInv;
1097 pMean *= noSamplesInv;
1098 break;
1099 }
1100 default: {
1101 TERMM(1, "fwh mean state type not implemented");
1102 break;
1103 }
1104 }
1105 }
1106
1107 // Set surface geometry properties
1108 // TODO: This assumes a stationary sampling surface! Impliment the version for moving surface...
1109 MFloat nx = normal[0];
1110 MFloat ny = normal[1];
1111 MFloat nz = normal[2];
1115
1117 // Get surface (f=0) velocity --> based on the ABSOLUTE frame of reference
1121
1122 // Get primative variables --> based on the ABSOLUTE frame of reference
1123 MFloat up = 0.0;
1124 MFloat vp = 0.0;
1125 MFloat wp = 0.0;
1126 MFloat pp = 0.0;
1127 MFloat rho = 0.0;
1139 rho = rho_inf + pp; // = rho_inf + p'/(c_inf * c_inf)
1140 }
1141
1142 // Compute sources
1146
1151
1152 // Transform the source terms to the flow direction coordinate
1153 transformCoordinatesToFlowDirection3D(&m_surfaceData.variables(segmentId, FWH_TIME::surfMX, tau),
1156 transformCoordinatesToFlowDirection3D(&m_surfaceData.variables(segmentId, FWH_TIME::srcLX, tau),
1159
1160 /*// Debugging
1161 MFloat* surfCoordinates = &m_surfaceData.surfaceCoordinates(segmentId);
1162 if(segmentId == 0 && tau == 0) {
1163 m_log << "tau = " << tau
1164 << ", surfCoordinates[0] = " << surfCoordinates[0]
1165 << ", surfCoordinates[1] = " << surfCoordinates[1]
1166 << ", surfCoordinates[2] = " << surfCoordinates[2]
1167 << ", nx = " << normal[0]
1168 << ", ny = " << normal[1]
1169 << ", nz = " << normal[2]
1170 << ", surfVX = " << surfVX
1171 << ", surfVY = " << surfVY
1172 << ", surfVZ = " << surfVZ
1173 << ", up = " << up
1174 << ", vp = " << vp
1175 << ", wp = " << wp
1176 << ", rho = " << rho
1177 << ", pp = " << pp
1178 << ", un = " << un
1179 << ", vn = " << vn
1180 << ", deltaUn = " << deltaUn
1181 << ", srcUn = " << std::setprecision(16) << m_surfaceData.variables(segmentId, FWH_TIME::srcUn, tau)
1182 << ", srcLX = " << std::setprecision(16) << m_surfaceData.variables(segmentId, FWH_TIME::srcLX, tau)
1183 << ", srcLY = " << std::setprecision(16) << m_surfaceData.variables(segmentId, FWH_TIME::srcLY, tau)
1184 << ", srcLZ = " << std::setprecision(16) << m_surfaceData.variables(segmentId, FWH_TIME::srcLZ, tau)
1185 << std::endl;
1186 }*/
1187 });
1188 } else if constexpr(nDim == 2) {
1189 TERMM(1, "Time Domain FW-H not implimented in 2D: Requires free space Green function in 2D");
1190 } // End of dimension switch
1191}
◆ calcSurfaceIntegralsAtSourceTime()
template<MInt nDim>
|
private |
- Date
- 15.05.2024 \Reference: Brentner, K. S. and F. Farassat, 1998, Analytical Comparison ... doi.org/10.2514/2.558
- Parameters
-
[in] globalObserverId [out] m_surfaceData.variables(segmentId,FWH_TIME::timeObs,tau) observer time [out] m_surfaceData.variables(segmentId,FWH_TIME::fwhPP,tau) observer pressure signal
Definition at line 1866 of file acasolver.cpp.
1866 {
1867 // Get observer coordinates
1869
1870 // Calculate Mach number and the Prandtl-Glauert-Factor
1873
1874 // Loop over all the source timeSteps
1876 if constexpr(nDim == 3) {
1877 // SURFACE LOOP
1880
1881 // Geometric Distance & Coordinate Transformation
1882 MFloat distX = coord[0] - surfCoordinates[0];
1883 MFloat distY = coord[1] - surfCoordinates[1];
1884 MFloat distZ = coord[2] - surfCoordinates[2];
1885 transformCoordinatesToFlowDirection3D(&distX, &distY, &distZ);
1886
1887 // Effective Acoustic Distance
1888 // Reference: Brès, G., 2010, A Ffowcs Williams - Hawkings ... doi.org/10.2514/6.2010-3711
1891
1897
1898 // compute observer time
1901
1902 // Get surface Mach number and source terms
1910
1911 // Project variables on the radiation direction
1917
1918 // Compute time derivative in second order accuracy
1926
1927 // Project time derivative on the radiation direction
1930
1933 (dist * dotMr + Mr - mach * mach) * invsersDistSq * inverseDopplerFactor * inverseDopplerFactorSq;
1934
1935 // Compute Thickness Source
1938
1939 // Compute Loading Source
1943
1944 // Save the computed source value
1947 surfIntegralPPT_0 + surfIntegralPPT_1 + surfIntegralPPL_0 + surfIntegralPPL_1 + surfIntegralPPL_2;
1948 m_surfaceData.variables(segmentId, FWH_TIME::fwhPP, tau) *= m_surfaceData.surfaceArea(segmentId);
1949
1950 /*// Debugging
1951 if(segmentId == 0 && tau == 0) {
1952 m_log << "tau = " << tau
1953 << ", timeSource = " << timeSource
1954 << ", timeObserver = " << timeObserver
1955 << ", coord[0] = " << coord[0]
1956 << ", coord[1] = " << coord[1]
1957 << ", coord[2] = " << coord[2]
1958 << ", surfCoordinates[0] = " << surfCoordinates[0]
1959 << ", surfCoordinates[1] = " << surfCoordinates[1]
1960 << ", surfCoordinates[2] = " << surfCoordinates[2]
1961 << ", distX = " << distX
1962 << ", distY = " << distY
1963 << ", distZ = " << distZ
1964 << ", dist = " << dist
1965 << ", distUnitX = " << distUnitX
1966 << ", distUnitY = " << distUnitY
1967 << ", distUnitZ = " << distUnitZ
1968 << ", dL = " << m_surfaceData.surfaceArea(segmentId)
1969 << ", srcUn = " << srcUn
1970 << ", dotUn = " << dotUn
1971 << ", Lr = " << Lr
1972 << ", dotLr = " << dotLr
1973 << ", LM = " << LM
1974 << ", Mr = " << Mr
1975 << ", dotMr = " << dotMr
1976 << ", MSq = " << MSq
1977 << ", MrFactor0 = " << MrFactor0
1978 << ", MrFactor1 = " << MrFactor1
1979 << ", surfIntegralPPT_0= " << surfIntegralPPT_0
1980 << ", surfIntegralPPT_1= " << surfIntegralPPT_1
1981 << ", surfIntegralPPL_0= " << surfIntegralPPL_0
1982 << ", surfIntegralPPL_1= " << surfIntegralPPL_1
1983 << ", surfIntegralPPL_2= " << surfIntegralPPL_2
1984 << std::endl;
1985 }*/
1986 } // End of surface loop
1987 } else if constexpr(nDim == 2) {
1988 TERMM(1, "Time Domain FW-H not implimented in 2D: Requires free space Green function in 2D");
1989 } // End of the dimension switch
1990 }); // End of the time loop
1991}
MFloat calcTimeDerivative(const MInt segmentId, const MInt varId, const MInt tau)
Definition: acasolver.cpp:2048
MFloat & surfaceArea(const MInt id)
Accessor for surface area (or segment length in 2D).
Definition: acasurfacedatacollector.h:281
◆ calcSurfaceIntegralsForObserver()
template<MInt nDim>
|
private |
- Date
- 15.06.2023
- Parameters
-
[in] coord Location of the observer [out] p_complexVariables Output of complex p', length of 2*noSamples()
Definition at line 1573 of file acasolver.cpp.
1574 {
1576 std::vector<MInt> VarDimC;
1578 if constexpr(nDim == 3) {
1580 } else {
1582 }
1584 if constexpr(nDim == 3) {
1586 } else {
1588 }
1589 } else {
1590 TERMM(1, "Unknown acoustic method.");
1591 }
1592
1593 // Calculate Mach number
1595 // Prandtl-Glauert-Factor
1599
1600 if constexpr(nDim == 2) {
1603
1604 MFloat theta = 0;
1606 theta = atan2(May, Max);
1607 }
1608
1609 // Pre-computing sine and cosine terms
1612
1613 // INNER LOOP OVER ALL FREQUENCIES (with only positive frequencies because bessel is not valid
1614 // for negativ nor zero). Since nw = 0 is always a 0, it is left out in the loop --> start at
1615 // nw = 1
1617 // non-dimensionalized omega is already in m_frequencies
1619
1620 // Set everything to 0
1621 MFloat integralFxRe = 0.0;
1622 MFloat integralFxIm = 0.0;
1623 MFloat integralFyRe = 0.0;
1624 MFloat integralFyIm = 0.0;
1625 MFloat integralQRe = 0.0;
1626 MFloat integralQIm = 0.0;
1627
1628 // LOOP OVER ALL LOCAL SURFACE ELEMENTS
1630 // Calculation of Prandtl-Glauert-Coordinates
1631 // Division by reference length is leaved out bc it would cancel out later
1637
1640
1641 // Argument in Hankelfunction argH and its derivatives
1643 const MFloat dargHdxi = waveNumber * fbeta2 * fd * (-1.0 * X * cosTheta + betaSq * Y * sinTheta);
1644 const MFloat dargHdeta = waveNumber * fbeta2 * fd * (-1.0 * X * sinTheta - betaSq * Y * cosTheta);
1645
1646 // Argument for exp function and its derivaties
1651
1656
1659
1660 // Calculation of the Green's function. y0, j0, y1 and j1 are used for Hankel function
1663
1664 // Calculation of the derivation of the Greenfunction resp. xi and eta
1666 C1 * (dargEdxi * (-1.0 * cosargE * J0 - sinargE * Y0) + dargHdxi * (-1.0 * cosargE * Y1 + sinargE * J1));
1668 C1 * (dargEdxi * (cosargE * Y0 - sinargE * J0) + dargHdxi * (-1.0 * cosargE * J1 - sinargE * Y1));
1670 C1 * (dargEdeta * (-1.0 * cosargE * J0 - sinargE * Y0) + dargHdeta * (-1.0 * cosargE * Y1 + sinargE * J1));
1672 C1 * (dargEdeta * (cosargE * Y0 - sinargE * J0) + dargHdeta * (-1.0 * cosargE * J1 - sinargE * Y1));
1673
1674 // Get Grid size dL
1676 // Get local values
1679
1682
1685
1686 // Calculate integral terms
1687 integralFxRe += dL * (FxRe * dGreendXi_Re - FxIm * dGreendXi_Im);
1688 integralFxIm += dL * (FxRe * dGreendXi_Im + FxIm * dGreendXi_Re);
1689
1690 integralFyRe += dL * (FyRe * dGreendEta_Re - FyIm * dGreendEta_Im);
1691 integralFyIm += dL * (FyRe * dGreendEta_Im + FyIm * dGreendEta_Re);
1692
1693 integralQRe += -1.0 * waveNumber * dL * (QRe * greenValue_Im + QIm * greenValue_Re);
1694 integralQIm += waveNumber * dL * (QRe * greenValue_Re - QIm * greenValue_Im);
1695 } // end segmentId Loop
1696
1697 // sum up all integral terms
1700
1701 // Store pressure values in observerData
1702 p_complexVariables[2 * nw + 0] = integralRe;
1703 p_complexVariables[2 * nw + 1] = integralIm;
1704 });
1705 } else if constexpr(nDim == 3) {
1709
1710 MFloat theta = 0.0;
1711 MFloat alpha = 0.0;
1712 // Determination of the flow direction
1713 // Use approx(..., 0.0, MFloatEps) instead of Maz == 0.0 || May == 0.0. Compares Ma with a
1714 // really low number epsilon
1716 theta = atan2(Maz, Max);
1717 }
1719 alpha = asin(May / mach);
1720 }
1721
1722 // Precompute trigonometric values
1727
1728 // FREQUENCIE LOOP: only for positive frequencies without nw=0
1730 // Calculate wave number: Non-dimensionalized omega is already saved in m_frequencies
1732
1733 // Set values to zero
1734 MFloat integralFxRe = 0.0;
1735 MFloat integralFxIm = 0.0;
1736 MFloat integralFyRe = 0.0;
1737 MFloat integralFyIm = 0.0;
1738 MFloat integralFzRe = 0.0;
1739 MFloat integralFzIm = 0.0;
1740 MFloat integralQRe = 0.0;
1741 MFloat integralQIm = 0.0;
1742
1743 // SURFACE LOOP
1746
1753
1754 /*
1755 // Version without Prandtl-Glauert
1756 const MFloat X = x - xi;
1757 const MFloat Y = y - eta;
1758 const MFloat Z = z - zeta;
1759 const MFloat R = sqrt(POW2(X) + POW2(Y) + POW2(Z));
1760
1761 const MFloat arg = waveNumber * R;
1762 const MFloat sinarg = sin(arg);
1763 const MFloat cosarg = cos(arg);
1764 const MFloat factor1 = waveNumber / (4.0 * PI * POW2(R));
1765 const MFloat factor2 = 1.0 / (4.0 * PI * POW3(R));
1766 const MFloat factorRe = -factor1 * sinarg - factor2 * cosarg;
1767 const MFloat factorIm = -factor1 * cosarg + factor2 * sinarg;
1768
1769 const MFloat greenValue_Re = -cosarg / (4.0 * PI * R);
1770 const MFloat greenValue_Im = sinarg / (4.0 * PI * R);
1771 const MFloat dGreendXi_Re = X * factorRe;
1772 const MFloat dGreendXi_Im = X * factorIm;
1773 const MFloat dGreendEta_Re = Y * factorRe;
1774 const MFloat dGreendEta_Im = Y * factorIm;
1775 const MFloat dGreendZeta_Re = Z * factorRe;
1776 const MFloat dGreendZeta_Im = Z * factorIm;
1777 */
1778
1779 // Version with Prandtl-Glauert, cf. Lockard2002
1780 const MFloat X = (x - xi) * cosAlpha * cosTheta + (y - eta) * sinAlpha + (z - zeta) * cosAlpha * sinTheta;
1781 // The Y coordinate differs from the paper, as there is a typo of the sign in front of the (z-zeta) term
1782 const MFloat Y = -(x - xi) * sinAlpha * cosTheta + (y - eta) * cosAlpha - (z - zeta) * sinAlpha * sinTheta;
1786
1787 // partial d / partial xi
1789 const MFloat d_xi = (-X * cosAlpha * cosTheta + betaSq * (Y * sinAlpha * cosTheta + Z * sinTheta)) / d;
1791 const MFloat d_zeta = (-X * cosAlpha * sinTheta + betaSq * (-Y * sinAlpha * sinTheta - Z * cosTheta)) / d;
1792
1796
1797 // Pre-terms in the derivatives
1801
1802 // clang-format off
1805 const MFloat dGreendXi_Re =-P0 * sinArg * cosAlpha * cosTheta + P1 * d_xi * cosArg - P2 * d_xi * sinArg;
1806 const MFloat dGreendXi_Im = P0 * cosArg * cosAlpha * cosTheta + P1 * d_xi * sinArg + P2 * d_xi * cosArg;
1807 const MFloat dGreendEta_Re =-P0 * sinArg * sinAlpha + P1 * d_eta * cosArg - P2 * d_eta * sinArg;
1808 const MFloat dGreendEta_Im = P0 * cosArg * sinAlpha + P1 * d_eta * sinArg + P2 * d_eta * cosArg;
1809 const MFloat dGreendZeta_Re =-P0 * sinArg * cosAlpha * sinTheta + P1 * d_zeta * cosArg - P2 * d_zeta * sinArg;
1810 const MFloat dGreendZeta_Im = P0 * cosArg * cosAlpha * sinTheta + P1 * d_zeta * sinArg + P2 * d_zeta * cosArg;
1811 // clang-format on
1812
1814
1817
1820
1823
1826
1827 // Calculation of integrals
1828 integralFxRe += dL * (FxRe * dGreendXi_Re - FxIm * dGreendXi_Im);
1829 integralFxIm += dL * (FxRe * dGreendXi_Im + FxIm * dGreendXi_Re);
1830
1831 integralFyRe += dL * (FyRe * dGreendEta_Re - FyIm * dGreendEta_Im);
1832 integralFyIm += dL * (FyRe * dGreendEta_Im + FyIm * dGreendEta_Re);
1833
1834 integralFzRe += dL * (FzRe * dGreendZeta_Re - FzIm * dGreendZeta_Im);
1835 integralFzIm += dL * (FzRe * dGreendZeta_Im + FzIm * dGreendZeta_Re);
1836
1837 integralQRe += -1.0 * waveNumber * dL * (QRe * greenValue_Im + QIm * greenValue_Re);
1838 integralQIm += waveNumber * dL * (QRe * greenValue_Re - QIm * greenValue_Im);
1839 } // End of segmentId loop
1840
1841 // sum up all integral terms
1844
1845 // Saving Pressure in Observerdata variable p for each positive frequencie nw between 1 and
1846 // size/2+1. 0 = RE; 1 = IM
1847 p_complexVariables[2 * nw + 0] = integralRe;
1848 p_complexVariables[2 * nw + 1] = integralIm;
1849 });
1850 }
1851 // nw = 0 is not defined and therefore set to zero to avoid "NaN-errors"
1853 p_complexVariables[2 * nw + 0] = 0.0;
1854 p_complexVariables[2 * nw + 1] = 0.0;
1856}
MFloat & complexVariables(const MInt id, const MInt var)
Accessor for complex variables.
Definition: acasurfacedatacollector.h:366
◆ calcTimeDerivative()
template<MInt nDim>
|
private |
Definition at line 2048 of file acasolver.cpp.
2048 {
2049 // Compute time derivative in second order accuarate
2050 MFloat dotVar = 0.0;
2051
2052 if(tau == 0) {
2060 } else {
2063 }
2065}
◆ calculateObserverInFrequency()
Definition at line 2294 of file acasolver.cpp.
2294 {
2297 complexVariablesBuffer.fill(0.0);
2300 std::array<MFloat, nDim> coord = getGlobalObserverCoordinates(globalObserverId);
2301 calcSurfaceIntegralsForObserver(coord, p_complexVariables);
2302 applySymmetryBc(coord, p_complexVariables);
2303 combineSurfaceIntegrals(globalObserverId, p_complexVariables);
2304 }
2305}
void applySymmetryBc(const std::array< MFloat, nDim > coord, MFloat *const p_complexVariables)
Symmetry boundary condition.
Definition: acasolver.cpp:4547
void combineSurfaceIntegrals(const MInt globalObserverId, MFloat *const p_complexVariables)
Combine the FWH surface integrals of all ranks.
Definition: acasolver.cpp:2203
◆ calculateObserverInTime()
Definition at line 2309 of file acasolver.cpp.
2309 {
2310 // Set memory for min observer distance
2312 distMinBuffer.fill(0.0);
2314
2315 // Set memory for observer pressure signal
2317 perturbedPressureBuffer.fill(0.0);
2319
2321 calcMinDistanceForObservers(distMinObservers);
2322
2324 // m_log << "Compute for observer = " << globalObserverId << "/" << globalNoObservers() << std::endl;
2325 calcSurfaceIntegralsAtSourceTime(globalObserverId);
2326 interpolateFromSourceToObserverTime(distMinObservers[globalObserverId], p_perturbedFWH);
2327 combineSurfaceIntegralsInTime(globalObserverId, p_perturbedFWH);
2328 }
2329}
void interpolateFromSourceToObserverTime(const MFloat distMinObservers, MFloat *const p_perturbedFWH)
Interpolate the observer pressure signal from source time to observer time using the advanced time ap...
Definition: acasolver.cpp:2123
void combineSurfaceIntegralsInTime(const MInt globalObserverId, MFloat *const p_perturbedFWH)
Combine the FWH surface integrals of the time domain method of all ranks.
Definition: acasolver.cpp:2253
void calcSurfaceIntegralsAtSourceTime(const MInt globalObserverId)
Calculate the FWH surface integrals for a certain observer using the time-domain formulation.
Definition: acasolver.cpp:1866
void calcMinDistanceForObservers(MFloat *const distMinObservers)
Find the minimum distance from the surface elements to the observers.
Definition: acasolver.cpp:2069
◆ calcWindowFactors()
Definition at line 1195 of file acasolver.cpp.
1195 {
1197 MFloat normFactor = 1.0;
1198
1202 break;
1203 }
1206 MFloat sum_w_sq = 0.0;
1209 sum_w_sq += w * w;
1210 p_window[t] = w;
1211 }
1212 normFactor = 1.0 / sqrt(sum_w_sq / N);
1213 break;
1214 }
1217 MFloat sum_w_sq = 0.0;
1220 sum_w_sq += w * w;
1221 p_window[t] = w;
1222 }
1223 normFactor = 1.0 / sqrt(sum_w_sq / N);
1224 break;
1225 }
1228 MFloat sum_w_sq = 0.0;
1230 MFloat w = 0.0;
1231 if(t < N / 8.0 || t > 7.0 * N / 8.0) {
1232 w = 0.5 * (1.0 - cos(8.0 * PI * t / N));
1233 } else {
1234 w = 1.0;
1235 }
1236 sum_w_sq += w * w;
1237 p_window[t] = w;
1238 }
1239 normFactor = 1.0 / sqrt(sum_w_sq / N);
1240 break;
1241 }
1242 default: {
1244 break;
1245 }
1246 }
1248 m_windowNormFactor = normFactor;
1249}
◆ changeDimensionsObserverData()
Definition at line 4816 of file acasolver.cpp.
4816 {
4821 }
4822 }
4827 }
4828 }
4833 }
4834 }
4835 }
4840 }
4841 }
4842}
struct AcaSolver::ConversionFactors m_aca2Output
◆ changeDimensionsSurfaceData()
Definition at line 4784 of file acasolver.cpp.
4784 {
4785 // All variables which are read in need to be non-dimensionalized with the freestream c and rho.
4786 // Index 0 in this (and only in this) function descripes the stagnation state.
4788
4789 // APE or not APE? That is the question
4794
4795 // Apply norm factor to the variables
4804 }
4805 } // end of sample
4806 } // end of seg
4807
4808 // Norm factor for time:
4811 }
4812}
struct AcaSolver::ConversionFactors m_input2Aca
◆ checkNoSamplesPotencyOfTwo()
Definition at line 4864 of file acasolver.cpp.
4864 {
4865 // Check if noSamples is a potency of 2 to see if FFT can be used
4867 MInt prove = 1;
4869
4870 // if potency of 2: numSampl == 0; if not: numSampl == 1
4871 while(numSampl > 1) {
4872 numSampl /= 2;
4873 prove *= 2;
4874 }
4878 "be used, if not defined otherwise (transformationType = 0 for DFT)."
4879 << std::endl;
4880 } else { // prove != noSamples()
4882 m_transformationType = 0;
4884 "cannot be used, no matter what."
4885 << std::endl;
4886 }
4887}
◆ cleanUp()
◆ combineSurfaceIntegrals()
template<MInt nDim>
|
private |
- Parameters
-
[in] globalObserverId Global observer id [in] p_complexVariables Output of complex p', length of 2*noSamples()
Definition at line 2203 of file acasolver.cpp.
2203 {
2204 TRACE();
2206
2207 // Data is stored in unconvenient way (globalObserverId, 0) so that Allreduce won't work. Put Data in linear
2208 // sequence --> Allreduce --> Back in the convenient data structure
2209 // (globalObserverId, 0)
2210
2211 // First: Copy all Observerdata to a databuffer
2212 // obsDataSize is noSamples, because we have imaginary and real part but only for half of samples
2213 // (all positive frequencies including zero)
2217 std::vector<MFloat> obsDataBuff(dataBuffSize, 0.0);
2218
2221
2223 "Error in observer partitioning!");
2224
2225 // The actual copying: Whole observer Data from observer "globalObserverId" to the position globalObserverId *
2226 // obsDataSize in databuffer
2227 std::copy_n(p_complexVariables, halfNoSamplesRealAndIm, obsDataBuff.data());
2228 // obsDataBuff now contains all complex pressure values of all observers in linear sequence
2229
2231 // Sums up all parts of the surface integral and distributes the results back to all processes.
2232 MPI_Reduce(MPI_IN_PLACE, obsDataBuff.data(), dataBuffSize, type_traits<MFloat>::mpiType(), MPI_SUM, trgDomain,
2234 // obsDataBuff now contains the FWH-Integral value of each observer and each frequency in linear
2235 // sequence.
2236
2237 // Split the data buffer into the convenient data structure
2238 std::copy_n(obsDataBuff.data(), halfNoSamplesRealAndIm, &m_observerData.complexVariables(observerId, 0));
2239 } else {
2240 MPI_Reduce(obsDataBuff.data(), obsDataBuff.data(), dataBuffSize, type_traits<MFloat>::mpiType(), MPI_SUM, trgDomain,
2242 }
2243
2244 // At first the segments were split up on processers. Each processes computed
2245 // the integral for all observers for his segments. Then all integral were
2246 // summed up. Now, all processes have the correct (local) observer data.
2247
2249}
MInt getObserverDomainId(const MInt globalObserverId)
Get domain id of (global) observer.
Definition: acasolver.h:177
MInt a_localObserverId(const MInt globalObserverId)
Return local observer id.
Definition: acasolver.h:155
int MPI_Reduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, int root, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Reduce
Definition: mpioverride.cpp:912
◆ combineSurfaceIntegralsInTime()
template<MInt nDim>
|
private |
Definition at line 2253 of file acasolver.cpp.
2253 {
2254 TRACE();
2255
2256 // Data is stored in unconvenient way (globalObserverId, 0) so that Allreduce won't work. Put Data in linear
2257 // sequence --> Allreduce --> Back in the convenient data structure (globalObserverId, 0)
2258
2259 // First: Copy all Observerdata to a databuffer
2261 std::vector<MFloat> obsDataBuff(dataBuffSize, 0.0);
2262
2265
2267 "Error in observer partitioning!");
2268
2269 // The actual copying: Whole observer Data from observer "globalObserverId" to
2270 // the position globalObserverId * obsDataSize in databuffer
2271 std::copy_n(p_perturbedFWH, dataBuffSize, obsDataBuff.data());
2272 // obsDataBuff now contains all complex pressure values of all observers in linear sequence
2273
2275 // Sums up all parts of the surface integral and distributes the results back to all processes.
2276 MPI_Reduce(MPI_IN_PLACE, obsDataBuff.data(), dataBuffSize, type_traits<MFloat>::mpiType(), MPI_SUM, trgDomain,
2278 // obsDataBuff now contains the FWH-Integral value of each observer and each frequency in linear sequence.
2279
2280 // Split the data buffer into the convenient data structure
2282 } else {
2283 MPI_Reduce(obsDataBuff.data(), obsDataBuff.data(), dataBuffSize, type_traits<MFloat>::mpiType(), MPI_SUM, trgDomain,
2285 }
2286
2287 // At first the segments were split up on processers. Each processes computed
2288 // the integral for all observers for his segments. Then all integral were
2289 // summed up. Now, all processes have the correct (local) observer data.
2290}
◆ computeAccuracy()
Definition at line 4892 of file acasolver.cpp.
4892 {
4893 TERMM(1, "TODO check this function");
4895 std::cout << "8. Starting accuarcy calculation." << std::endl;
4896
4897 // These parameters have to be set
4898 // -------------------
4900 MInt noSegs = 100;
4901 noSegs = Context::getSolverProperty<MFloat>("noSegments_" + std::to_string(0), m_solverId, AT_, 0);
4902 if(noSegs == 0) {
4903 noSegs = Context::getSolverProperty<MFloat>("noSegments_" + std::to_string(0), m_solverId, AT_, 1);
4904 }
4905
4906 MFloat box_size = 2.0 * 2;
4907 box_size = Context::getSolverProperty<MFloat>("surfaceCoordinates_" + std::to_string(0), m_solverId, AT_, 0);
4908 box_size = 2.0 * std::fabs(box_size);
4909 // --------------------
4910
4911 // Alternative expression: SamplePoints per Period even with flow
4912 MFloat noPeriods = 1.0;
4914
4915 MFloat omega_factor = 4.0;
4916 omega_factor = Context::getSolverProperty<MFloat>("omega_factor", m_solverId, AT_, &omega_factor);
4917
4918 MFloat Delta_Seg = box_size / noSegs;
4919 if constexpr(nDim == 3) {
4920 Delta_Seg *= Delta_Seg;
4921 }
4922
4923 // Wave length lambda
4926
4927 // minimal resolution of one period according to the thesis.
4931
4932 // Read in the analytic solution for even periods to be able to compare with uneven periods
4934 std::vector<MFloat> Analytic_EvenPeriods_Amplitude(globalNoObservers() * noSamples() / 2, 0.0);
4935 std::vector<MFloat> Analytic_EvenPeriods_Frequencies(globalNoObservers() * noSamples() / 2, 0.0);
4936 std::vector<MFloat> RMS_EvenPeriods(globalNoObservers(), 0.0);
4937
4940 std::ifstream file;
4941 MString line;
4943 file.open(Name, std::ios::in);
4944 MInt count = 0;
4945 MFloat i1, i2, i3, i4, i5;
4946 if(file.is_open()) {
4948 std::stringstream linebuffer(line);
4949 linebuffer >> i1 >> i2 >> i3 >> i4 >> i5;
4950 Analytic_EvenPeriods[5 * count] = i1;
4951 Analytic_EvenPeriods[5 * count + 1] = i2;
4952 Analytic_EvenPeriods[5 * count + 2] = i3;
4953 Analytic_EvenPeriods[5 * count + 3] = i4;
4954 Analytic_EvenPeriods[5 * count + 4] = i5;
4955
4956 Analytic_EvenPeriods_Amplitude[count] = i4;
4957 Analytic_EvenPeriods_Frequencies[count] = i1;
4958 count += 1;
4959 }
4960 } else {
4961 std::cerr << "Fehler beim Öffnen der Datei: " << Name << "." << std::endl;
4962 }
4963
4964 std::ifstream file1;
4965 MString line1;
4967 file1.open(Name1, std::ios::in);
4968 MInt count1 = 0;
4969 MFloat i_1, i_2;
4970 if(file1.is_open()) {
4972 std::stringstream linebuffer(line1);
4973 linebuffer >> i_1 >> i_2;
4974 RMS_EvenPeriods[count1] = i_2;
4975 count1 += 1;
4976 }
4977 } else {
4978 std::cerr << "Fehler beim Öffnen der Datei: " << Name1 << "." << std::endl;
4979 }
4980 }
4981 // RMS PRESSURE
4982 std::vector<MFloat> p_rms(globalNoObservers(), 0.0);
4983
4987 }
4988 p_rms[obs] = sqrt(p_rms[obs] / noSamples());
4989 }
4990
4991 // Quality (Variance) check
4992 MFloat max_diff = 0.0;
4993 MFloat max_diff_2 = 0.0;
4994 MFloat diff_normalized = 0.0;
4995 MFloat diff_2_normalized = 0.0;
4996 MFloat sum_diff_sq = 0.0;
4997 MFloat sum_diff_2_sq = 0.0;
4998 MFloat sum_analytic_sq = 0.0;
4999 MFloat sum_a_ep_sq = 0.0;
5000 std::vector<MFloat> diff(observer, 0.0);
5001 std::vector<MFloat> diff_2(observer, 0.0);
5003 diff[i] = std::fabs(m_rmsP_Analytic[i] - p_rms[i]);
5004 sum_diff_sq += diff[i] * diff[i];
5005 sum_analytic_sq += m_rmsP_Analytic[i];
5006 diff_normalized = diff[i] / m_rmsP_Analytic[i] * 100;
5007
5008 // In case of uneven periods
5010 diff_2[i] = std::fabs(RMS_EvenPeriods[i] - p_rms[i]);
5011 sum_diff_2_sq += diff_2[i] * diff_2[i];
5012 sum_a_ep_sq += RMS_EvenPeriods[i];
5013 diff_2_normalized = diff_2[i] / RMS_EvenPeriods[i] * 100;
5014 }
5015 if(diff_normalized >= max_diff) {
5017 max_diff = diff_normalized; // [%]
5018 }
5019 if(EP_Samples == 64) {
5020 max_diff_2 = diff_2_normalized; //[%]
5021 }
5022 }
5023 }
5027
5028 // In case of uneven periods
5029 MFloat Q_C_UnevenPeriods = 0.0;
5033 Q_C_UnevenPeriods = 100 * value_top_UP / value_bot_UP;
5034 }
5035
5036 // Frequency, Amplitude comparison for all observers
5038 std::vector<MFloat> max_Ampl(globalNoObservers(), 0.0);
5039 std::vector<MFloat> max_Ampl_2(globalNoObservers(), 0.0);
5040 std::vector<MFloat> max_Ampl_a(globalNoObservers(), 0.0);
5041 std::vector<MFloat> max_Ampl_A_EvenPeriods(globalNoObservers(), 0.0);
5042 std::vector<MFloat> max_Ampl_Diff(globalNoObservers(), 0.0);
5043 std::vector<MFloat> max_Ampl_Diff_2(globalNoObservers(), 0.0);
5044 std::vector<MFloat> max_Ampl_Diff_normalized(globalNoObservers(), 0.0);
5045 std::vector<MFloat> max_Ampl_Diff_normalized_2(globalNoObservers(), 0.0);
5046 std::vector<MFloat> freq_At_Max_Ampl(globalNoObservers(), 0.0);
5047 std::vector<MFloat> freq_At_Max_Ampl_2(globalNoObservers(), 0.0);
5048 std::vector<MFloat> freq_At_Max_Ampl_a(globalNoObservers(), 0.0);
5049 std::vector<MFloat> freq_At_Max_Ampl_A_EvenPeriods(globalNoObservers(), 0.0);
5050 std::vector<MFloat> freq_Diff(globalNoObservers(), 0.0);
5051 std::vector<MFloat> freq_Diff_2(globalNoObservers(), 0.0);
5052 std::vector<MFloat> freq_Diff_normalized(globalNoObservers(), 0.0);
5053 std::vector<MFloat> freq_Diff_normalized_2(globalNoObservers(), 0.0);
5054 std::vector<MFloat> mean_Ampl_analytic(globalNoObservers(), 0.0);
5055 std::vector<MFloat> sum_Diff_Amplitudes_UnevenPeriods(globalNoObservers(), 0.0);
5056 std::vector<MFloat> SideLobes(globalNoObservers(), 0.0);
5057 MFloat sum_max_Ampl_Diff_normalized = 0.0;
5058 MFloat sum_freq_Diff_normalized = 0.0;
5059 MFloat sum_max_Ampl_Diff_normalized_2 = 0.0;
5060 MFloat sum_freq_Diff_normalized_2 = 0.0;
5061 MFloat sum_SideLobes = 0.0;
5062 MFloat Amplitude_EvenPeriods = 0.0;
5071 // Comparison with even periods in case uneven periods are used
5073 Amplitude_EvenPeriods = Analytic_EvenPeriods_Amplitude[obs * noSamples() / 2 + (i - 1)];
5074
5075 // Get Frequency at which the analytic source for even periods oscillates
5076 if(Amplitude_EvenPeriods > max_Ampl_A_EvenPeriods[obs]) {
5077 freq_At_Max_Ampl_A_EvenPeriods[obs] = Analytic_EvenPeriods_Frequencies[obs * noSamples() / 2 + (i - 1)];
5078 max_Ampl_A_EvenPeriods[obs] = Amplitude_EvenPeriods;
5079 }
5080 }
5081
5082 // Get Frequency at which analytic source oscillates
5083 if(Amplitude_a > max_Ampl_a[obs]) {
5084 freq_At_Max_Ampl_a[obs] = m_frequencies[i];
5085 max_Ampl_a[obs] = Amplitude_a;
5086 }
5087
5088 // Get max Amplitude computed and its frequency
5089 if(Amplitude_c > max_Ampl[obs]) {
5090 freq_At_Max_Ampl[obs] = m_frequencies[i];
5091 max_Ampl[obs] = Amplitude_c;
5092 }
5093
5094 // Side lobes
5095 SideLobes[obs] += std::fabs(Amplitude_EvenPeriods - Amplitude_c);
5096 }
5097 // Normalize the maximum difference in amplitude and frequency with the maximum analytical
5098 // ampl. and freq. Index 2 is for uneven Periods comparison
5099 max_Ampl_Diff[obs] = std::fabs(max_Ampl_a[obs] - max_Ampl[obs]);
5100 freq_Diff[obs] = std::fabs(freq_At_Max_Ampl[obs] - freq_At_Max_Ampl_a[obs]);
5101
5102 max_Ampl_Diff_normalized[obs] = max_Ampl_Diff[obs] / max_Ampl_a[obs];
5103 freq_Diff_normalized[obs] = freq_Diff[obs] / freq_At_Max_Ampl_a[obs];
5104
5105 sum_max_Ampl_Diff_normalized += max_Ampl_Diff_normalized[obs];
5106 sum_freq_Diff_normalized += freq_Diff_normalized[obs];
5107
5109 max_Ampl_Diff_2[obs] = std::fabs(max_Ampl_A_EvenPeriods[obs] - max_Ampl[obs]);
5110 freq_Diff_2[obs] = std::fabs(freq_At_Max_Ampl[obs] - freq_At_Max_Ampl_A_EvenPeriods[obs]);
5111 sum_SideLobes += (SideLobes[obs] / (noSamples() / 2));
5112
5113 max_Ampl_Diff_normalized_2[obs] = max_Ampl_Diff_2[obs] / max_Ampl_A_EvenPeriods[obs];
5114 freq_Diff_normalized_2[obs] = freq_Diff_2[obs] / freq_At_Max_Ampl_A_EvenPeriods[obs];
5115
5116 // In case of uneven and even periods comparison
5117 sum_max_Ampl_Diff_normalized_2 += max_Ampl_Diff_normalized_2[obs];
5118 sum_freq_Diff_normalized_2 += freq_Diff_normalized_2[obs];
5119 }
5120 }
5123
5124 // In case of uneven and even periods comparison
5125 MFloat Side_Lobes = 0.0;
5126 MFloat Freq_aver_Diff_UnevenPeriods = 0.0;
5127 MFloat Ampl_aver_Diff_UnevenPeriods = 0.0;
5129 Ampl_aver_Diff_UnevenPeriods = 100 * sum_max_Ampl_Diff_normalized_2 / globalNoObservers(); //[%]
5132 // Those minor deviations are due to the writing out in a file and re-reading it again
5133 // (Analytic Even Periods). Thus, they match perfectly but a deviation is falsly shown
5134 if(Freq_aver_Diff_UnevenPeriods < 1e-10) {
5135 Freq_aver_Diff_UnevenPeriods = 0.0;
5136 }
5137 }
5138
5139 // OUTPUT OF max Amplitude at given frequencie
5141 fileFreq += std::to_string(noSegs);
5142 fileFreq += ".dat";
5143 std::ofstream outfile;
5144 outfile.open(fileFreq);
5145 outfile << Y << " " << Freq_aver_Diff << " " << Ampl_aver_Diff << " " << Side_Lobes << std::endl;
5146 outfile.close();
5147
5148 // Freq for samples
5150 fileF += std::to_string(noSamples());
5151 fileF += ".dat";
5152 outfile.open(fileF);
5153 outfile << R_min << " " << Freq_aver_Diff << " " << Ampl_aver_Diff << " " << Side_Lobes << std::endl;
5154
5155 outfile.close();
5156
5158 // Freq for periods
5160 fileFr += std::to_string(noPeriods);
5161 fileFr += ".dat";
5162 outfile.open(fileFr);
5163 outfile << noPeriods << " " << Freq_aver_Diff_UnevenPeriods << " " << Ampl_aver_Diff_UnevenPeriods << " "
5164 << Side_Lobes << " " << Ampl_aver_Diff << std::endl;
5165 outfile.close();
5166 }
5167
5168 // OUTPUT IF ONE VARIES NUMBER OF SEGMENTS
5170 ofile += std::to_string(noSegs);
5171 ofile += ".dat";
5172 outfile.open(ofile);
5173 outfile << Y << " " << Q_C << " " << max_diff << " " << Z << " " << Side_Lobes << std::endl;
5174 outfile.close();
5175
5176 // OUTPUT IF ONE VARIES NUMBER OF SAMPLES
5178 fileName += std::to_string(noSamples());
5179 fileName += ".dat";
5180 outfile.open(fileName);
5181 outfile << R_min << " " << Q_C << " " << max_diff << " " << Side_Lobes << std::endl;
5182 outfile.close();
5183
5185 // OUTPUT IF ONE VARIES NUMBER OF SAMPLES
5187 fileN += std::to_string(noPeriods);
5188 fileN += ".dat";
5189 outfile.open(fileN);
5190 outfile << noPeriods << " " << Q_C_UnevenPeriods << " " << max_diff_2 << " " << Side_Lobes << std::endl;
5191 outfile.close();
5192 }
5193
5194 // OUTPUT IF ONE VARIES NUMBER OF SAMPLES
5196 fileNa += std::to_string(omega_factor);
5197 fileNa += ".dat";
5198 outfile.open(fileNa);
5199 outfile << omega_factor << " " << Q_C << " " << max_diff << " " << Side_Lobes << std::endl;
5200 outfile.close();
5201
5202 std::cout << "8. Accuarcy calculation finished." << std::endl;
5203 }
5204}
std::vector< MFloat > m_rmsP_Analytic
Analytic rms pressure for accuracy function.
Definition: acasolver.h:415
std::vector< MFloat > m_Analyticfreq
Analytic Fourier Transformation, also used for Validation.
Definition: acasolver.h:412
std::istream & getline(std::istream &input, std::string &line)
Definition: cpptoml.h:1519
◆ computeDt()
Definition at line 4845 of file acasolver.cpp.
4845 {
4846 // Compute time step size
4848
4849 // Check for equidistant timesteps - necessary for the computation!
4853 if(relErr > 1e-10) {
4854 std::ostringstream err;
4855 err << std::setprecision(12) << std::scientific << "Error: time step size mismatch; m_dt = " << m_dt
4856 << "; delta = " << delta << "; relErr = " << relErr << "; i = " << i << "; times[i] = " << m_times[i]
4858 TERMM(1, err.str());
4859 }
4860 }
4861}
◆ DiscreteFourierTransform()
template<MInt nDim>
|
private |
- Parameters
-
dataArray Pointer to complex data array. size Size of the complex data array (number of samples). dir Direction of transformation: 1=forward; -1=backward. result Optional storage for transformed data, used if inPlace==false. inPlace Determines if result overwrites the input data or stored in result array.
Definition at line 1500 of file acasolver.cpp.
1501 {
1502 TRACE();
1506 std::vector<MFloat> resultArray;
1507
1508 if(inPlace) {
1509 resultArray.resize(size2); // Temporary solution storage
1510 resultPtr = &resultArray[0];
1511 } else {
1512 ASSERT(result != nullptr, "Pointer to result storage is nullptr.");
1513 resultPtr = result; // Store solution in provided storage
1514 }
1515
1516 if(dir == 1) { // FOURIER FORWARD
1520 // cos(-a) = cos(a) and sin(-a) = -sin(a)
1521 // Actual Discrete Fourier Transform
1522 resultPtr[2 * k] += dataArray[2 * i] * cos(k * con * i) + dataArray[2 * i + 1] * sin(k * con * i);
1523 resultPtr[2 * k + 1] += -1.0 * dataArray[2 * i] * sin(k * con * i) + dataArray[2 * i + 1] * cos(k * con * i);
1524 }
1525 resultPtr[2 * k] /= sizeF;
1526 resultPtr[2 * k + 1] /= sizeF;
1527 }
1528 } else if(dir == -1) { // FOURIER BACKWARD
1532 resultPtr[2 * t] += dataArray[2 * k] * cos(arg * k) - dataArray[2 * k + 1] * sin(arg * k);
1533 resultPtr[2 * t + 1] += dataArray[2 * k] * sin(arg * k) + dataArray[2 * k + 1] * cos(arg * k);
1534 }
1535 }
1536 } else {
1537 TERMM(1, "Please insert 1 for forward or -1 for backward");
1538 }
1539
1540 if(inPlace) {
1541 // Copy results back into dataArray overwriting the input data
1542 std::copy_n(&resultPtr[0], size2, &dataArray[0]);
1543 }
1544}
◆ distributeObservers()
Evenly distribute observers on all ranks.
Definition at line 2333 of file acasolver.cpp.
2333 {
2334 TRACE();
2338 m_noObservers = minNoObsPerDomain + 1;
2340 } else {
2341 m_noObservers = minNoObsPerDomain;
2343 }
2344}
◆ FastFourierTransform()
template<MInt nDim>
|
private |
Fast Fourier transformation.
Reference: 'Numerical Recipes in C', Cambridge, p. 507 ff.
- Parameters
-
dataArray Pointer to complex data array. size Size of the complex data array (number of samples). dir Direction of transformation: 1=forward; -1=backward. result Optional storage for transformed data, used if inPlace==false. inPlace Determines if result overwrites the input data or stored in result array.
Definition at line 1411 of file acasolver.cpp.
1412 {
1413 TRACE();
1414 TERMM_IF_NOT_COND(dir == 1 || dir == -1, "Invalid direction");
1415 TERMM_IF_NOT_COND(size % 2 == 0, "FFT: not a multiple of 2 samples.");
1416
1420
1421 // Set pointer to result storage
1423 // Copy input data if not performing transform in-place (overwriting input)
1424 if(!inPlace) {
1425 ASSERT(result != nullptr, "Pointer to result storage is nullptr.");
1426 std::copy_n(&dataArray[0], size2, result);
1427 }
1428
1429#define SWAP(a, b) \
1430 rtemp = (a); \
1431 (a) = (b); \
1432 (b) = rtemp
1433 // Bit reversing: Reorder by inverse bit notation, e.g. first number: 0000 --> 0000, second
1434 // number: 0001 -> 1000
1436 MInt n = nn << 1;
1437 MInt m;
1438 MFloat rtemp;
1439 MInt j = 1;
1441 if(j > i) {
1442 SWAP(data[j - 1], data[i - 1]);
1443 SWAP(data[j], data[i]);
1444 }
1445 m = nn;
1446 while(m >= 2 && j > m) {
1447 j -= m;
1448 m >>= 1;
1449 }
1450 j += m;
1451 }
1452
1453 // Fourier Transformation (Danielson-Lanczos)
1455 MFloat wtemp, wr, wpr, wpi, wi, theta, tempr, tempi;
1457 MUlong mmax = 2;
1458 while(nLong > mmax) {
1459 istep = mmax << 1;
1460 theta = isign * (6.28318530717959 / mmax);
1461 wtemp = sin(0.5 * theta);
1462 wpr = -2.0 * wtemp * wtemp;
1463 wpi = sin(theta);
1464 wr = 1.0;
1465 wi = 0.0;
1467 for(i = p; i <= nLong; i += istep) {
1468 k = i + mmax;
1469 tempr = wr * data[k - 1] - wi * data[k];
1470 tempi = wr * data[k] + wi * data[k - 1];
1471 data[k - 1] = data[i - 1] - tempr;
1472 data[k] = data[i] - tempi;
1473 data[i - 1] += tempr;
1474 data[i] += tempi;
1475 }
1476 wr = (wtemp = wr) * wpr - wi * wpi + wr;
1477 wi = wi * wpr + wtemp * wpi + wi;
1478 }
1479 mmax = istep;
1480 }
1481
1482 // Dimensionierung
1483 if(dim) {
1485 data[z] /= size;
1486 }
1487 }
1488}
◆ finalizeInitSolver()
◆ genDipoleAnalytic2D()
template<MInt nDim>
|
private |
2D analytical dipole generation
Definition at line 4160 of file acasolver.cpp.
4160 {
4161 MFloat x = coord[0];
4164
4173
4174 // Precompute time constant values/factors
4177
4178 MFloat k = omega / c0;
4179
4189
4193
4197
4201
4204
4206 C1 * omega * (-1.0 * sinarg * J0 + cosarg * Y0) + C2 * dHdx * omega * (cosarg * J1 + sinarg * Y1);
4207 const MFloat dPhidx = -1.0 * C1 * (dargdx * (sinarg * J0 - cosarg * Y0) + dHdx * (cosarg * J1 + sinarg * Y1))
4208 + C2
4209 * (dHdxdx * (sinarg * J1 - cosarg * Y1) + dHdx * dargdx * (cosarg * J1 + sinarg * Y1)
4210 + dHdx * dHdx * (sinarg * dJ1 - cosarg * dY1));
4211
4213 + C2 * (dHdxdy * (sinarg * J1 - cosarg * Y1) + dHdx * dHdy * (sinarg * dJ1 - cosarg * dY1));
4215
4216 if(param.perturbed) {
4217 vars.u[sample] = dPhidx;
4218 vars.v[sample] = dPhidy;
4219 vars.p[sample] = pp;
4220 } else {
4221 vars.u[sample] = dPhidx + u0;
4222 vars.v[sample] = dPhidy;
4223 vars.p[sample] = pp + p0;
4224 vars.rho[sample] = rho0 + pp / (c0 * c0);
4225 }
4226 transformBackToOriginalCoordinate2D(&vars.u[sample], &vars.v[sample]);
4227 }
4228}
void transformBackToOriginalCoordinate2D(MFloat *const fX, MFloat *const fY)
Definition: acasolver.cpp:2020
void transformCoordinatesToFlowDirection2D(MFloat *const fX, MFloat *const fY)
Coordinate transformation.
Definition: acasolver.cpp:1994
◆ genDipoleAnalytic3D()
template<MInt nDim>
|
private |
Definition at line 4233 of file acasolver.cpp.
4233 {
4234 MFloat x = coord[0];
4236 MFloat z = coord[2];
4238
4247
4248 MFloat k = omega / c0;
4249
4258
4262
4266
4268 // Note: adapted from acoupot3d/acou_pot.cpp::calc_p/ux/uy/uz_analytic()
4271
4274
4276 C1 * omega * fd3 * x * sinarg + C1 * omega * fd2 * x * C2 * cosarg - C1 * omega * C2 * (mach / d) * cosarg;
4277 const MFloat dPhidx = -1.0 * C1 * P1 * cosarg + C1 * x * fd3 * dargdx * sinarg + C1 * P2 * C2 * sinarg
4278 + C1 * x * fd2 * C2 * dargdx * cosarg - C1 * C2 * P3 * mach * sinarg
4279 - C1 * C2 * (mach / d) * dargdx * cosarg;
4280 const MFloat dPhidy = C1 * x * 3 * y * beta2 * fd3 * fd2 * cosarg + C1 * x * fd3 * dargdy * sinarg
4281 - C1 * x * C2 * 2 * y * beta2 * fd2 * fd2 * sinarg + C1 * x * C2 * fd2 * dargdy * cosarg
4282 + C1 * C2 * y * beta2 * mach * fd3 * sinarg - C1 * C2 * (mach / d) * dargdy * cosarg;
4283 const MFloat dPhidz = C1 * x * 3 * z * beta2 * fd3 * fd2 * cosarg + C1 * x * fd3 * dargdz * sinarg
4284 - C1 * x * C2 * 2 * z * beta2 * fd2 * fd2 * sinarg + C1 * x * C2 * fd2 * dargdz * cosarg
4285 + C1 * C2 * z * beta2 * mach * fd3 * sinarg - C1 * C2 * (mach / d) * dargdz * cosarg;
4287
4288 if(param.perturbed) {
4289 vars.u[sample] = dPhidx;
4290 vars.v[sample] = dPhidy;
4291 vars.w[sample] = dPhidz;
4292 vars.p[sample] = pp;
4293 } else {
4294 vars.u[sample] = dPhidx + u0;
4295 vars.v[sample] = dPhidy;
4296 vars.w[sample] = dPhidz;
4297 vars.p[sample] = pp + p0;
4298 vars.rho[sample] = rho0 + pp / (c0 * c0);
4299 }
4300 transformBackToOriginalCoordinate3D(&vars.u[sample], &vars.v[sample], &vars.w[sample]);
4301 }
4302}
void transformBackToOriginalCoordinate3D(MFloat *const fX, MFloat *const fY, MFloat *const fZ)
Definition: acasolver.cpp:2033
◆ generateSurfaceCircle()
Definition at line 3853 of file acasolver.cpp.
3853 {
3854 TRACE();
3855
3856 // Weighting factor for surface averaging
3858
3859 if constexpr(nDim != 2) {
3860 TERMM(1, "Only implemented in 2D.");
3861 }
3862
3863 // TODO labels:ACA Fix this
3864 TERMM(1, "The 2D implementation is not tested and seems not to work properly.");
3865
3866 // Creating Circle with radius and number of Elements
3868 totalNoSurfaceElements = Context::getSolverProperty<MInt>("noSegments_" + std::to_string(m_surfaceIds[sId]),
3870
3871 MFloat radius = 1.0;
3873
3875
3876 std::vector<MFloat> surfaceNormal(2, 0.0);
3880 // Add a new surface element
3882 // Check that the total size of m_surfaceData is correct
3884
3885 // First surface point
3888
3889 // Second surface point
3892
3893 // Segment center
3896
3897 // Determine segment length
3901
3902 // Store segment center coordinates
3904 surfCoordinates[0] = xcenter;
3905 surfCoordinates[1] = ycenter;
3906
3907 // Store segment area (weighted with surface averaging factor)
3909
3910 // Calculate surface normal vector. Simply the gradient and normed
3911 surfaceNormal[0] = 2 * xcenter / (sqrt(4 * xcenter * xcenter + 4 * ycenter * ycenter));
3912 surfaceNormal[1] = 2 * ycenter / (sqrt(4 * xcenter * xcenter + 4 * ycenter * ycenter));
3913
3914 // Store segment normal vector
3916 surfNormal[0] = surfaceNormal[0];
3917 surfNormal[1] = surfaceNormal[1];
3918 }
3919}
std::vector< MInt > m_noSurfaceElements
Number of (local) surface elements for each surface.
Definition: acasolver.h:310
std::vector< MFloat > m_surfaceWeightingFactor
Weighting factor for each surface (for surface averaging)
Definition: acasolver.h:329
MInt localSurfaceOffset(const MInt sId)
Return local offset of a surface in m_surfaceData.
Definition: acasolver.h:141
constexpr MInt size() const
Return size (i.e., currently used number of nodes)
Definition: container.h:89
◆ generateSurfaceData()
Generate surface data for an analytical test case.
Definition at line 3924 of file acasolver.cpp.
3924 {
3925 TRACE();
3926
3929
3930 // TODO labels:ACA Read in analytic or semi-analytic
3932 /* analytic = Context::getSolverProperty<MBool>("analyticSourceGeneration", m_solverId, AT_,
3933 * &analytic); */
3934
3935 // Initialize every variable with zeros
3936 std::fill_n(&m_surfaceData.variables(0, 0, 0), totalNoSurfaceElements() * noSamples() * noVars, 0.0);
3937
3939 // Initialize times (constant timestep size)
3940
3943 }
3944
3945 // TODO labels:ACA,toenhance add possibility for superposition of multiple sources with different source parameters
3947
3949 // TODO labels:ACA,toenhance change structure (performance)?
3950 // Iterate over all surface segments and all observers.
3953
3954 SourceVars sourceVars;
3961 } else {
3966 }
3967
3968 if(analytic) {
3969 // Full analytic formulation
3970 if constexpr(nDim == 2) {
3972 case 0: {
3974 break;
3975 }
3976 case 1: {
3978 break;
3979 }
3980 case 2: {
3982 break;
3983 }
3984 case 3: {
3986 break;
3987 }
3988 default: {
3989 TERMM(1, "source type not implemented");
3990 break;
3991 }
3992 }
3993 } else {
3995 case 0: {
3997 break;
3998 }
3999 case 1: {
4001 break;
4002 }
4003 case 2: {
4005 break;
4006 }
4007
4008 default: {
4009 TERMM(1, "source type not implemented");
4010 break;
4011 }
4012 }
4013 }
4014 } else {
4015 TERMM(1, "Error: semi-analytic source generation not implemented.");
4016 // TODO labels:ACA Semi-analytic formulation (numeric differentiation of acoustic potential)
4017 // 1. generate the acoustic potential (without derivatives) for the current
4018 // segment and a surrounding stencil of points
4019 // 2. use finite differences to:
4020 // 2.1 compute the derivatives to get the correct acoustic potential (eg. dipole/quadrupole)
4021 // 2.2 compute the derivatives to get the variables p', u', ...
4022 /* if constexpr(nDim == 2) { */
4023 /* genSurfaceDataNumeric2D(segmentId, m_sourceParameters); */
4024 /* } else { */
4025 /* genSurfaceDataNumeric3D(segmentId, m_sourceParameters); */
4026 /* } */
4027 }
4028 } // end of segmentId for loop
4029 } // end of sourceType loop
4030}
◆ generateSurfacePlane()
Definition at line 3593 of file acasolver.cpp.
◆ generateSurfaces()
Generation of surfaces (for analytical cases).
Definition at line 3498 of file acasolver.cpp.
◆ genMonopoleAnalytic2D()
template<MInt nDim>
|
private |
2D analytical monopole generation
Definition at line 4035 of file acasolver.cpp.
4035 {
4036 MFloat x = coord[0];
4039
4048
4049 // Precompute time constant values/factors
4051
4065
4067 // Note: adapted from acoupot3d/acou_pot.cpp::calc_p/ux/uy/uz_analytic()
4070
4073
4076 -1.0 * mach * amplkf4beta3 * (cosarg * J0 + sinarg * Y0) + amplkf4beta3 * xfd * (sinarg * J1 - cosarg * Y1);
4079
4080 if(param.perturbed) {
4081 vars.u[sample] = dPhidx;
4082 vars.v[sample] = dPhidy;
4083 vars.p[sample] = pp;
4084 } else {
4085 vars.u[sample] = dPhidx + u0;
4086 vars.v[sample] = dPhidy;
4087 vars.p[sample] = pp + p0;
4088 vars.rho[sample] = rho0 + pp / (c0 * c0);
4089 }
4090 transformBackToOriginalCoordinate2D(&vars.u[sample], &vars.v[sample]);
4091 }
4092}
◆ genMonopoleAnalytic3D()
template<MInt nDim>
|
private |
Definition at line 4097 of file acasolver.cpp.
4097 {
4098 MFloat x = coord[0];
4100 MFloat z = coord[2];
4102
4111
4112 // Precompute time constant values/factors
4113 // Wavenumber k needs to be dimensionless
4114 MFloat k = omega / c0;
4115
4126
4128 // Note: adapted from acoupot3d/acou_pot.cpp::calc_p/ux/uy/uz_analytic()
4131
4134
4136 const MFloat dPhidx = -1.0 * ampf4pid * xfd2 * cosarg + ampf4pid * k * fbeta2 * sinarg * (xfd - mach);
4140
4141 if(param.perturbed) {
4142 vars.u[sample] = dPhidx;
4143 vars.v[sample] = dPhidy;
4144 vars.w[sample] = dPhidz;
4145 vars.p[sample] = pp;
4146 } else {
4147 vars.u[sample] = dPhidx + u0;
4148 vars.v[sample] = dPhidy;
4149 vars.w[sample] = dPhidz;
4150 vars.p[sample] = pp + p0;
4151 vars.rho[sample] = rho0 + pp / (c0 * c0);
4152 }
4153 transformBackToOriginalCoordinate3D(&vars.u[sample], &vars.v[sample], &vars.w[sample]);
4154 }
4155}
◆ genQuadrupoleAnalytic2D()
template<MInt nDim>
|
private |
2D quadrupole analytic
Definition at line 4307 of file acasolver.cpp.
4307 {
4308 MFloat x = coord[0];
4311
4320
4321 MFloat k = omega / c0;
4322
4330
4333
4338
4341
4345
4346
4350
4354
4358
4361
4363 + K2 * dHdxdy * cosarg * omega * J1 + K2 * dHdxdy * sinarg * omega * Y1
4364 + K2 * dHdx * dHdy * cosarg * omega * J0 - K2 * dHdxfH * dHdy * cosarg * omega * J1
4365 + K2 * dHdx * dHdy * sinarg * omega * Y0 - K2 * dHdxfH * dHdy * sinarg * omega * Y1;
4366
4368 -1.0 * K1 * dHdydx * cosarg * J1 + K1 * dHdy * sinarg * dargdx * J1 - K1 * dHdy * cosarg * dJ1 * dHdx
4369 - K1 * dHdydx * sinarg * Y1 - K1 * dHdy * cosarg * dargdx * Y1 - K1 * dHdy * sinarg * dY1 * dHdx
4370 + K2 * dHdxdydx * sinarg * J1 + K2 * dHdxdy * cosarg * dargdx * J1 + K2 * dHdxdy * sinarg * dJ1 * dHdx
4371 - K2 * dHdxdydx * cosarg * Y1 + K2 * dHdxdy * sinarg * dargdx * Y1 - K2 * dHdxdy * cosarg * dY1 * dHdx
4372 + K2 * dHdxdx * dHdy * sinarg * J0 + K2 * dHdx * dHdydx * sinarg * J0 + K2 * dHdx * dHdy * cosarg * dargdx * J0
4373 - K2 * dHdx * dHdy * sinarg * J1 * dHdx - K2 * dHdxfHdx * dHdy * sinarg * J1
4374 - K2 * dHdxfH * dHdydx * sinarg * J1 - K2 * dHdxfH * dHdy * cosarg * dargdx * J1
4375 - K2 * dHdxfH * dHdy * sinarg * dJ1 * dHdx - K2 * dHdxdx * dHdy * cosarg * Y0 - K2 * dHdx * dHdydx * cosarg * Y0
4376 + K2 * dHdx * dHdy * sinarg * dargdx * Y0 + K2 * dHdx * dHdy * cosarg * Y1 * dHdx
4377 + K2 * dHdxfHdx * dHdy * cosarg * Y1 + K2 * dHdxfH * dHdydx * cosarg * Y1
4378 - K2 * dHdxfH * dHdy * sinarg * dargdx * Y1 + K2 * dHdxfH * dHdy * cosarg * dY1 * dHdx;
4379
4381 -1.0 * K1 * dHdydy * cosarg * J1 - K1 * dHdy * cosarg * dJ1 * dHdy - K1 * dHdydy * sinarg * Y1
4382 - K1 * dHdy * sinarg * dY1 * dHdy + K2 * dHdxdydy * sinarg * J1 + K2 * dHdxdy * sinarg * dJ1 * dHdy
4383 - K2 * dHdxdydy * cosarg * Y1 - K2 * dHdxdy * cosarg * dY1 * dHdy + K2 * dHdxdy * dHdy * sinarg * J0
4384 + K2 * dHdx * dHdydy * sinarg * J0 - K2 * dHdx * dHdy * sinarg * J1 * dHdy - K2 * dHdxfHdy * dHdy * sinarg * J1
4385 - K2 * dHdxfH * dHdydy * sinarg * J1 - K2 * dHdxfH * dHdy * sinarg * dJ1 * dHdy
4386 - K2 * dHdxdy * dHdy * cosarg * Y0 - K2 * dHdx * dHdydy * cosarg * Y0 + K2 * dHdx * dHdy * cosarg * Y1 * dHdy
4387 + K2 * dHdxfHdy * dHdy * cosarg * Y1 + K2 * dHdxfH * dHdydy * cosarg * Y1
4388 + K2 * dHdxfH * dHdy * cosarg * dY1 * dHdy;
4389
4391
4392 if(param.perturbed) {
4393 vars.u[sample] = dPhidx;
4394 vars.v[sample] = dPhidy;
4395 vars.p[sample] = pp;
4396 } else {
4397 vars.u[sample] = dPhidx + u0;
4398 vars.v[sample] = dPhidy;
4399 vars.p[sample] = pp + p0;
4400 vars.rho[sample] = rho0 + pp / (c0 * c0);
4401 }
4402 transformBackToOriginalCoordinate2D(&vars.u[sample], &vars.v[sample]);
4403 }
4404}
◆ genQuadrupoleAnalytic3D()
template<MInt nDim>
|
private |
Definition at line 4409 of file acasolver.cpp.
4409 {
4410 MFloat x = coord[0];
4412 MFloat z = coord[2];
4414
4423
4424 MFloat k = omega / c0;
4425
4435
4439
4441 // Note: adapted from acoupot3d/acou_pot.cpp::calc_p/ux/uy/uz_analytic()
4444
4447
4451
4453 C1 * 3 * y * beta2 * (fd2 * fd3 - 5 * x * x * fd2 * fd2 * fd3) * cosarg
4454 - C1 * 3 * x * y * beta2 * fd2 * fd3 * sinarg * dargdx
4455 - C1 * 3 * k * y * (fd2 * fd2 - 4 * x * x * fd3 * fd3) * sinarg
4456 - C1 * 3 * k * x * y * fd2 * fd2 * cosarg * dargdx - C1 * k * C2 * y * (fd3 - 3 * x * x * fd3 * fd2) * cosarg
4458 + C1 * mach * k * y * fd3 * cosarg * dargdx - C1 * mach * k * C2 * y * 2 * x * fd2 * fd2 * cosarg
4459 - C1 * mach * k * C2 * y * fd2 * sinarg * dargdx;
4460
4463 - C1 * 3 * x * y * beta * beta * fd3 * fd2 * sinarg * dargdy
4465 - C1 * 3 * k * x * y * fd2 * fd2 * cosarg * dargdy
4466 - C1 * k * C2 * x * (fd3 - 3 * beta2 * y * y * fd2 * fd3) * cosarg + C1 * k * C2 * x * y * fd3 * sinarg * dargdy
4467 + C1 * mach * k * (fd3 - 3 * beta2 * y * y * fd2 * fd3) * sinarg + C1 * mach * k * y * fd3 * cosarg * dargdy
4469 - C1 * mach * k * C2 * y * fd2 * sinarg * dargdy;
4470
4472 -1.0 * C1 * 15 * x * y * z * beta2 * beta2 * fd2 * fd2 * fd3 * cosarg
4473 - C1 * 3 * x * y * beta2 * fd2 * fd3 * sinarg * dargdz + C1 * 12 * k * x * y * z * beta2 * fd3 * fd3 * sinarg
4474 - C1 * 3 * k * x * y * fd2 * fd2 * cosarg * dargdz + C1 * 3 * k * k * x * y * z * fd2 * fd3 * cosarg
4475 + C1 * k * C2 * x * y * fd3 * sinarg * dargdz - C1 * 3 * mach * k * y * z * beta2 * fd2 * fd3 * sinarg
4476 + C1 * mach * k * y * fd3 * cosarg * dargdz - C1 * 2 * mach * k * k * y * z * fd2 * fd2 * cosarg
4477 - C1 * mach * k * C2 * y * fd2 * sinarg * dargdz;
4478
4480
4481 if(param.perturbed) {
4482 vars.u[sample] = dPhidx;
4483 vars.v[sample] = dPhidy;
4484 vars.w[sample] = dPhidz;
4485 vars.p[sample] = pp;
4486 } else {
4487 vars.u[sample] = dPhidx + u0;
4488 vars.v[sample] = dPhidy;
4489 vars.w[sample] = dPhidz;
4490 vars.p[sample] = pp + p0;
4491 vars.rho[sample] = rho0 + pp / (c0 * c0);
4492 }
4493 transformBackToOriginalCoordinate3D(&vars.u[sample], &vars.v[sample], &vars.w[sample]);
4494 }
4495}
◆ genVortexConvectionAnalytic2D()
template<MInt nDim>
|
private |
Calculate solution vars for analytical inviscid convecting vortex.
- Date
- 11.05.2023
A vortex is convected in an inviscid flow with M_infty in x-direction. The analytical solution is known to not generate noise. Hence, all noise generate is of erroneous nature.
Definition at line 4506 of file acasolver.cpp.
4507 {
4512
4523 const MFloat pp = pow(1 + (m_gamma - 1) * exp(1 - r * r) / (8.0 * PI * PI), m_gamma / (m_gamma - 1.0)) - 1.0;
4524 if(param.perturbed) {
4525 vars.u[sample] = up;
4526 vars.v[sample] = vp;
4527 vars.p[sample] = pp;
4528 } else {
4529 vars.u[sample] = param.Ma * cs + up;
4530 vars.v[sample] = vp;
4531 vars.p[sample] = pp;
4533 }
4534 }
4535}
◆ getGlobalObserverCoordinates()
template<MInt nDim>
|
inlineprivate |
Definition at line 169 of file acasolver.h.
169 {
170 std::array<MFloat, nDim> coord;
172 coord[dir] = a_globalObserverCoordinate(globalObserverId, dir);
173 }
174 return coord;
175 };
MFloat a_globalObserverCoordinate(const MInt globalObserverId, const MInt dir)
Accessor for global observer coordinates.
Definition: acasolver.h:165
◆ getObserverDomainId()
template<MInt nDim>
|
inlineprivate |
Definition at line 177 of file acasolver.h.
177 {
180 if(globalObserverId < remainderObs * (minNoObsPerDomain + 1)) {
181 return globalObserverId / (minNoObsPerDomain + 1);
182 } else {
183 return (globalObserverId - remainderObs) / minNoObsPerDomain;
184 }
185 };
◆ getSurfaceDataFileName()
template<MInt nDim>
|
inlineprivate |
Return the name of a surface data \input[in] surfaceId id of the surface element \input[in] fileNo index of the file.
Definition at line 2913 of file acasolver.cpp.
2913 {
2914 TRACE();
2917 std::ostringstream os;
2918 os << m_inputDir << baseFileName << std::to_string(m_surfaceIds[surfaceId]) << "_" << std::setw(8)
2919 << std::setfill('0') << fileNo << baseFileNameSuffix;
2921}
◆ globalNoObservers()
Definition at line 153 of file acasolver.h.
◆ hasSurface()
Definition at line 132 of file acasolver.h.
◆ initConversionFactors()
Definition at line 4721 of file acasolver.cpp.
4721 {
4722 // TODO labels:ACA IMHO it does not make sense to have a different dimensions
4723 // or even quantities (wave number or frequency) in the output. The output
4724 // should have the same dimensions for a certain solver to avoid confusion
4725 // when switching between different input solver.
4726 std::stringstream ss;
4739 // .. and corresponding inverse factors
4741 // Default definition is correct -> dump out in ACA dimension
4749 }
4750 break;
4751 }
4753 // To reduce some confusion of the LB user. In LB all quantities are made
4754 // non-dimensional with dxLb (cell length at max level), density at
4755 // infinity, and xi_0=dxLb/dtLb. But, the time in the surface sampling
4756 // file is given by t_infty = t * u_infty / L_FV , L_FV=1.
4759 // TODO labels:ACA How to choose dxLb (question of input only)
4760 // const MFloat dxLb = Context::getSolverProperty<MFloat>("dxLb", m_solverId, AT_, 0);
4761 // m_input2Aca.length = dxLb;
4765 break;
4766 }
4767 default: {
4769 }
4770 }
4772 ss << " Conversion factors input2Aca | aca2Output:" << std::endl;
4778 printMessage(ss.str());
4779}
◆ initObserverPoints()
Initialize the observer points (load/generate coordinates and init collector)
Definition at line 4584 of file acasolver.cpp.
◆ initParallelization()
Initialize the parallelization, i.e., distribute all surface segments among ranks create communicators etc.
Definition at line 2672 of file acasolver.cpp.
2672 {
2673 TRACE();
2676
2681
2682 std::vector<MInt> totalNoSegments(noSurfaces());
2683
2684 // TODO labels:ACA maybe do this only on rank 0
2687 std::vector<MInt> dummyVec{};
2688 // Determine total number of segments for each surface
2689 // const MInt totalNoSegments = readNoSegmentsSurfaceGeneration(sId, dummyVec);
2690 totalNoSegments[sId] = readNoSegmentsSurfaceGeneration(sId, dummyVec);
2691 }
2692 } else {
2693 m_totalNoSamples = -1;
2694
2695 // Determine total number of samples and number of segments for each surface from input data files
2699
2700 // Check that number of samples matches for all surfaces
2702 "number of samples mismatch for different surfaces");
2704 }
2705
2707
2708 // Check noSamples/stride/offset
2710 TERMM_IF_COND(m_noSamples < 1 || m_sampleStride < 1 || m_sampleOffset < 0 || usedSampleRange > m_totalNoSamples,
2711 "Error: number of samples and stride need to be > 0, offset >= 0 and samples*stride+offset cannot "
2712 "exceed the total number of samples (noSamples="
2714 + std::to_string(m_sampleOffset) + ", samples*stride+offset=" + std::to_string(usedSampleRange)
2716
2717 // Only even number of samples
2721 std::cout << "Note: Number of Samples was uneven which is not feasible. The new number of samples is: "
2722 << m_noSamples << std::endl;
2723 }
2724 }
2725 }
2726
2727 // Store global number of segments for all surfaces
2728 std::copy(totalNoSegments.begin(), totalNoSegments.end(), &m_noSurfaceElementsGlobal[0]);
2729
2730 // Determine distribution of surface segments among all processes
2731 const MInt globalNoSegments = std::accumulate(totalNoSegments.begin(), totalNoSegments.end(), 0);
2732
2737 }
2741
2742 MInt localNoSegments = -1;
2743
2744 // Distribute all segments equally among all domains. This means that a process can have segments
2745 // of different surfaces, e.g., when using only a small number of processes. However, for larger
2746 // production cases it makes sense to limit each process to have only segments of a single surface
2747 // if the I/O of data of different surfaces should be fully parallel.
2749 // Average number of segments (rounded down)
2750 const MInt avgNoSegments = floor(static_cast<MFloat>(globalNoSegments) / static_cast<MFloat>(noDomains()));
2751 localNoSegments = avgNoSegments;
2752 // Equally distribute missing segments
2754 localNoSegments += 1;
2755 }
2756 } else { // Restrict each rank to have only segments of one surface to improve performance in production cases
2757 std::vector<MInt> noDomainsPerSurface(noSurfaces());
2758 // Average number of segments
2759 const MFloat avgNoSegments = static_cast<MFloat>(globalNoSegments) / static_cast<MFloat>(noDomains());
2760 // Distribute ranks on surfaces
2762 noDomainsPerSurface[i] = std::max(1, static_cast<MInt>(std::floor(totalNoSegments[i] / avgNoSegments)));
2763 }
2764 // Assign remaining ranks to the surfaces with the highest number of segments per rank
2765 while(noDomains() - std::accumulate(noDomainsPerSurface.begin(), noDomainsPerSurface.end(), 0) > 0) {
2766 MFloat maxWeight = 0.0;
2767 MInt maxIndex = -1;
2770 if(weight > maxWeight) {
2771 maxWeight = weight;
2772 maxIndex = i;
2773 }
2774 }
2775 noDomainsPerSurface[maxIndex]++;
2776 }
2777 const MInt noDomainsSum = std::accumulate(noDomainsPerSurface.begin(), noDomainsPerSurface.end(), 0);
2778 TERMM_IF_NOT_COND(noDomainsSum == noDomains(),
2779 "number of domains mismatch " + std::to_string(noDomainsSum) + " " + std::to_string(noDomains()));
2780
2781 // Determine the local surface
2782 MInt localSurfaceId = -1;
2783 MInt localSurfaceDomainId = -1;
2786 localSurfaceId = i;
2787 localSurfaceDomainId = domainId() - sumDomains;
2788 break;
2789 }
2790 sumDomains += noDomainsPerSurface[i];
2791 }
2792
2793 const MInt avgNoSegmentsSurface = floor(totalNoSegments[localSurfaceId] / noDomainsPerSurface[localSurfaceId]);
2794 // Determine local number of segments of this surface
2795 localNoSegments = avgNoSegmentsSurface;
2796 if(localSurfaceDomainId
2797 < totalNoSegments[localSurfaceId] - noDomainsPerSurface[localSurfaceId] * avgNoSegmentsSurface) {
2798 localNoSegments += 1;
2799 }
2800 }
2801
2802 // Communicate number of segments per domain
2803 std::vector<MInt> noSegmentsPerDomain(noDomains());
2806
2807 // Sanity check
2808 const MInt globalCount = std::accumulate(noSegmentsPerDomain.begin(), noSegmentsPerDomain.end(), 0);
2809 TERMM_IF_NOT_COND(globalCount == globalNoSegments, "globalNoSegments mismatch");
2810
2811 // Compute domain offsets
2812 std::vector<MInt> domainSegmentOffsets(std::max(noDomains() + 1, 1));
2813 domainSegmentOffsets[0] = 0;
2815 domainSegmentOffsets[i] = domainSegmentOffsets[i - 1] + noSegmentsPerDomain[i - 1];
2816 }
2817
2819
2820 MInt surfaceOffset = 0;
2821 // Check which surfaces are present on the local rank
2826
2827 MInt noLocalSegments = 0;
2828 MInt localSegmentsOffset = 0;
2829
2830 if(start < surfaceOffset && end >= surfaceOffset) {
2831 // First local segment belongs to a previous surface but last local segment belongs to current
2832 // or a subsequent surface
2833 noLocalSegments = std::min(totalNoSegments[sId], end - surfaceOffset);
2834 localSegmentsOffset = 0;
2835 } else if(start >= surfaceOffset && start < nextSurfaceOffset) {
2836 // First local segment belongs to this surface
2837 noLocalSegments = std::min(nextSurfaceOffset - start, end - start);
2838 localSegmentsOffset = start - surfaceOffset;
2839 }
2840
2841 m_noSurfaceElements[sId] = noLocalSegments;
2842 m_surfaceElementOffsets[sId] = localSegmentsOffset;
2843
2844 // Store in communication buffer
2845 noSegments_(domainId(), sId) = noLocalSegments;
2846
2847 // Set new offset for next surface
2848 surfaceOffset = nextSurfaceOffset;
2849 }
2850
2851 MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &noSegments_[0], noSurfaces(), type_traits<MInt>::mpiType(),
2853
2857
2858 // Create MPI communicators for each surface
2861 MInt noDomains_ = 0;
2862
2863 // Determine list of subdomains with segments of this surface
2865 if(noSegments_(d, sId) > 0) {
2866 domains[noDomains_] = d;
2867 noDomains_++;
2868 }
2869 }
2870
2871 m_log << "Create MPI communicator for surface #" << m_surfaceIds[sId] << " with " << noDomains_
2872 << " ranks:" << std::endl;
2873#ifndef NDEBUG
2875 m_log << " * local rank #" << d << " global rank #" << domains[d] << " with " << noSegments_(domains[d], sId)
2876 << " segments" << std::endl;
2877 }
2878#endif
2879
2880 // Obtain new group
2881 MPI_Group globalGroup, localGroup;
2883 MPI_Group_incl(globalGroup, noDomains_, &domains[0], &localGroup, AT_);
2884
2885 // Create new communicator and clean up
2887 "m_mpiCommSurface[" + std::to_string(sId) + "]");
2888
2889 MPI_Group_free(&globalGroup, AT_);
2890 MPI_Group_free(&localGroup, AT_);
2891
2893 // If surface is active on current domain, determine local rank and number of domains
2896
2898 } else {
2899 // Reset for inactive surface
2900 m_domainIdSurface[sId] = -1;
2901 m_noDomainsSurface[sId] = -1;
2902 }
2903 }
2904
2906}
void readNoSegmentsAndSamples(const MInt surfaceId, MInt &noSegments, MInt &noSamples, MString &inputFileName)
Read the number of segments and samples for the given surface id.
Definition: acasolver.cpp:2994
MInt m_sampleOffset
Offset in input data files from which to start reading samples.
Definition: acasolver.h:341
std::vector< MString > m_surfaceInputFileName
Surface input file names that were used for sampling.
Definition: acasolver.h:332
std::vector< MInt > m_surfaceElementOffsets
Local surface element offsets for each surface.
Definition: acasolver.h:319
MInt readNoSegmentsSurfaceGeneration(const MInt sId, std::vector< MInt > &noSegments, const MInt lengthCheck=-1)
Determine the number of segments to be generated for a surface, returns the total number of segments ...
Definition: acasolver.cpp:3541
std::vector< MInt > m_noSurfaceElementsGlobal
Total number of surface elements for each surface.
Definition: acasolver.h:316
MInt noSurfaceElements(const MInt surfaceId)
Return number of (local) surfaces elements for a surface.
Definition: acasolver.h:135
int MPI_Comm_create(MPI_Comm comm, MPI_Group group, MPI_Comm *newcomm, const MString &name, const MString &varname)
same as MPI_Comm_create, but updates the number of MPI communicators
Definition: mpioverride.cpp:249
int MPI_Group_incl(MPI_Group group, int n, const int ranks[], MPI_Group *newgroup, const MString &name)
same as MPI_Group_incl
Definition: mpioverride.cpp:1615
int MPI_Comm_group(MPI_Comm comm, MPI_Group *group, const MString &name, const MString &varname)
same as MPI_Comm_group
Definition: mpioverride.cpp:398
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
int MPI_Group_free(MPI_Group *group, const MString &name)
same as MPI_Group_free
Definition: mpioverride.cpp:1649
◆ initSolver()
Implements Solver.
Definition at line 28 of file acasolver.cpp.
28 {
29 TRACE();
30
31 // Initialize all solver-wide timers
32 initTimers();
33 // Input/output
34 setInputOutputProperties();
35 // Numerical method
36 setNumericalProperties();
37}
◆ initTimers()
Definition at line 41 of file acasolver.cpp.
41 {
42 TRACE();
43
44 // Invalidate timer ids
46
47 // Create timer group & timer for solver, and start the timer
48 NEW_TIMER_GROUP_NOCREATE(m_timerGroup, "AcaSolver (solverId = " + std::to_string(m_solverId) + ")");
51
52 // Create regular solver-wide timers
54
55 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::InitParallelization], "initParallelization", m_timers[Timers::Run]);
56 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::ReadSurfaceData], "readSurfaceData", m_timers[Timers::Run]);
57
59 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::CalcSourceTerms], "calcSourceTerms", m_timers[Timers::Run]);
62 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::CalcSurfaceIntegrals], "calcSurfaceIntegrals", m_timers[Timers::Run]);
63 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::CombineSurfaceIntegrals], "combineSurfaceIntegrals", m_timers[Timers::Run]);
64 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::BacktransformObservers], "backtransformObservers", m_timers[Timers::Run]);
65 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::SaveObserverSignals], "saveObserverSignals", m_timers[Timers::Run]);
66 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Postprocessing], "postprocessing", m_timers[Timers::Run]);
67}
◆ interpolateFromSourceToObserverTime()
template<MInt nDim>
|
private |
Definition at line 2123 of file acasolver.cpp.
2123 {
2125
2126 // Loop over all the observer timeSteps
2128 // Set the observer time
2130
2131 // initialize surface integral
2132 MFloat sumFwhPP = 0.0;
2133
2134 // Surface element loop
2136 // initialize interpolation value
2137 MFloat interpolateFwhPP = 0.0;
2138 MInt t_lower = 0;
2140
2141 // Read the observer lower/upper bound
2144
2145 // Search in FWH_TIME::timeObs
2146 if(timeObserver < timeObserver_lower || timeObserver > timeObserver_upper) {
2147 // In case the timeObserver is out of bound
2148 interpolateFwhPP = 0.0;
2150 // In case the timeObserver is at the lower bound
2153 // In case the timeObserver is at the upper bound
2155 } else {
2156 // In the case of: timeObserver_lower < timeObserver < timeObserver_upper
2157 // Use binary search algorithm to search
2158
2159 MInt t_range = t_upper - t_lower;
2160 do {
2161 // Compute t_mid
2162 const MFloat t_mid = (t_range % 2 == 0) ? ((t_lower + t_upper) / 2) : ((t_lower + t_upper - 1) / 2);
2164
2166 // In case the timeObserver is at the middle
2168 break;
2169 }
2170
2171 // Update t_upper and t_lower
2172 if(timeObserver < timeObserver_mid) {
2173 t_upper = t_mid;
2174 } else if(timeObserver > timeObserver_mid) {
2175 t_lower = t_mid;
2176 }
2177 t_range = t_upper - t_lower;
2178
2179 // For the smallest un-dividable region, do linear interpolation to get the observer fwhPP
2180 if(t_range == 1) {
2185 const MFloat fwhPP_slop = (fwhPP_upper - fwhPP_lower) / (timeObserver_upper - timeObserver_lower);
2186 interpolateFwhPP = fwhPP_lower + (timeObserver - timeObserver_lower) * fwhPP_slop;
2187 break;
2188 }
2189 } while(t_range > 1);
2190 } // End of the search in FWH_TIME::timeObs
2191
2192 sumFwhPP += interpolateFwhPP;
2193 } // End of the surface loop
2194
2195 p_perturbedFWH[t] = inverseFourPI * sumFwhPP;
2196 }); // End of the observer time loop (parallel)
2197}
◆ isActive()
Return if the solver is active on this rank; needs to be implemented in derived solver since access to gridproxy not possible here
Reimplemented from Solver.
Definition at line 103 of file acasolver.h.
103{ return true; } // All ranks are active
◆ loadObserverSignals()
Definition at line 2477 of file acasolver.cpp.
2477 {
2478 TRACE();
2481
2484
2485 std::stringstream ss;
2487 printMessage(ss.str());
2489
2490 // read times and frequencies on rank 0 and broadcast
2492 file.setOffset(noSamples(), 0);
2493 } else {
2494 file.setOffset(0, 0);
2495 }
2498
2499 MPI_Bcast(&m_times[0], noSamples(), type_traits<MFloat>::mpiType(), 0, mpiComm(), AT_, "&m_times[0]");
2500 MPI_Bcast(&m_frequencies[0], noSamples(), type_traits<MFloat>::mpiType(), 0, mpiComm(), AT_, "&m_frequencies[0]");
2501
2502 // Read observer data in time domain
2503 file.setOffset(obsCount, obsOffset, 2);
2505
2506 file.readArray(&pressure_time[0], "pressure_time");
2510 }
2511 }
2512
2513 // Read observer data in frequency domain
2514 file.setOffset(obsCount, obsOffset, 2);
2515 ScratchSpace<MFloat> pressure_frequency(std::max(obsCount, 1), 2 * noSamples(), AT_, "pressure_frequency");
2516 file.readArray(&pressure_frequency[0], "pressure_frequency");
2517
2522 }
2523 }
2524
2526}
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
◆ localSurfaceOffset()
Definition at line 141 of file acasolver.h.
141 {
143 }
◆ noInternalCells()
◆ noObservers()
Definition at line 149 of file acasolver.h.
◆ noSamples()
Definition at line 146 of file acasolver.h.
◆ noSurfaceElements()
Definition at line 135 of file acasolver.h.
◆ noSurfaces()
Definition at line 129 of file acasolver.h.
◆ noVariables()
◆ observerOffset()
Definition at line 151 of file acasolver.h.
◆ postprocessObserverSignals()
Definition at line 2531 of file acasolver.cpp.
2531 {
2532 TRACE();
2535
2536 // Loop over all postprocessing operations
2539 switch(operation) {
2542 break;
2543 }
2546 break;
2547 }
2550 break;
2551 }
2554 break;
2555 }
2557 TERMM_IF_NOT_COND(m_generateSurfaceData, "Only useful for analytical testcases with generated surface data.");
2558 calcObsPressureAnalytic();
2559 break;
2560 }
2561 default: {
2562 std::cerr << " "
2563 << "skippinginvalid postprocessing operation: " << operation << std::endl;
2564 break;
2565 }
2566 }
2567 }
2568
2570 post->init(noObservers(), globalNoObservers(), observerOffset(), noSamples(), m_globalObserverCoordinates.data(),
2574 }
2575 post->finish();
2576 }
2577
2579}
ACA post processing class for overall sound pressure calculation.
Definition: acapostprocessing.hpp:10
ACA post processing class for absoulte pressure calculation.
Definition: acapostprocessing.hpp:68
ACA post processing class for calculation of root mean square pressure.
Definition: acapostprocessing.hpp:127
ACA post processing class for sound pressure level calculation.
Definition: acapostprocessing.hpp:208
std::vector< std::unique_ptr< AcaPostProcessing > > m_post
Definition: acasolver.h:272
◆ postTimeStep()
◆ preTimeStep()
◆ printMessage()
Definition at line 203 of file acasolver.h.
◆ readNoSegmentsAndSamples()
template<MInt nDim>
|
private |
Definition at line 3000 of file acasolver.cpp.
◆ readNoSegmentsSurfaceGeneration()
template<MInt nDim>
|
private |
Definition at line 3543 of file acasolver.cpp.
◆ readSurfaceData()
Definition at line 2925 of file acasolver.cpp.
2925 {
2926 TRACE();
2929
2930 // Set data sizes in surface data collector
2934
2935 // Reset surface data collector to total number of local surface elements
2937 TERMM_IF_COND(totalNoSurfaceElements_ < 1, "Less surface elements than computational ranks is critical.");
2939
2940 // Set storage size for times and frequencies
2943
2944 // Read in the surface data and store in m_surfaceData
2945 // or generate the surface elements and the corresponding analytical input data
2947 generateSurfaces();
2948 generateSurfaceData();
2949
2950 // Testing: store generated data
2953 storeSurfaceData(sId);
2954 }
2955 }
2957 // Nothing to do if just postprocessing is enabled
2958 } else {
2960 // Read data if surface present on local rank
2962 readSurfaceDataFromFile(sId);
2963
2964 // Testing: store data for validating the I/O
2966 storeSurfaceData(sId);
2967 }
2968 }
2969 }
2970
2971 { // Compare the loaded times on all ranks
2974
2975 MPI_Bcast(×Buffer[0], noSamples(), type_traits<MFloat>::mpiType(), 0, mpiComm(), AT_, "×Buffer[0]");
2976
2978 TERMM_IF_COND(!approx(timesBuffer[i], m_times[i], MFloatEps), "Error: times for surfaces do not match.");
2979 }
2980 }
2981
2982 // The simulation non-dimensionalises with different variables! For the FWH simulation one needs
2983 // to change how the numbers are non-dimensionalized
2984 changeDimensionsSurfaceData();
2985 computeDt();
2986 }
2987
2989}
void storeSurfaceData(const MInt surfaceId)
Store the surface data. Note: used for validating the data input and for the output of generated data...
Definition: acasolver.cpp:3428
void changeDimensionsSurfaceData()
Change of the non-dimensionalization of input data from the stagnation state (0) non-dim....
Definition: acasolver.cpp:4784
void readSurfaceDataFromFile(const MInt surfaceId)
Read the data for one surface.
Definition: acasolver.cpp:3064
MBool hasSurface(const MInt surfaceId)
Return if this rank has data of this surface.
Definition: acasolver.h:132
void reset()
Reset, re-create data structures with given capacity, and set size to zero.
Definition: acasurfacedatacollector.h:232
void setNoComplexVars(const MInt noComplexVars_)
Set number of complex variables.
Definition: acasurfacedatacollector.h:180
void setNoSamples(const MInt noSamples_)
Set number of samples.
Definition: acasurfacedatacollector.h:189
◆ readSurfaceDataFromFile()
Definition at line 3064 of file acasolver.cpp.
3064 {
3065 TRACE();
3067
3071
3073
3074 // Create surface elements and store its data in m_surfaceData
3079 MFloat totalSurfaceMeasure = 0.0;
3080 MFloat totalSurfaceMeasureWeighted = 0.0;
3081
3085 // Add a new surface element
3087 // Check that the total size of m_surfaceData is correct
3089
3092
3093 // Store segment area/length (weighted with surface averaging factor)
3095 totalSurfaceMeasure += l_measure[i];
3096 totalSurfaceMeasureWeighted += (weight * l_measure[i]);
3097
3098 // Set coordinates and normal vector
3100 surfCoordinates[j] = l_coordinates[i * nDim + j];
3101 surfNormal[j] = l_normalVector[i * nDim + j];
3102 }
3103 }
3104
3105 // Compute total area/length of surface
3108 MPI_Allreduce(MPI_IN_PLACE, &totalSurfaceMeasureWeighted, 1, type_traits<MFloat>::mpiType(), MPI_SUM,
3111 std::cout << "Total " << measureName << " of surface #" << surfaceId << ": " << totalSurfaceMeasure
3112 << " (weighted: " << totalSurfaceMeasureWeighted << ")" << std::endl;
3113 }
3114 };
3115
3118 + Context::getSolverProperty<MString>(
3121
3125
3126 // Read surface area/length
3127 file.setOffset(count, offset);
3128 file.readArray(&measure[0], measureName);
3129
3130 // Read coordinates and normal vectors
3131 file.setOffset(count, offset, 2);
3132 file.readArray(&coordinates[0], "coordinates");
3133 file.readArray(&normalVector[0], "normalVector");
3134
3135 // Create surface elements and store its data in m_surfaceData
3136 createSurfaceElements(surfaceId, measure, coordinates, normalVector);
3137
3138 // Read times only on surface local rank 0 and broadcast to other ranks in this communicator
3141 // Read full array, then copy selected range to m_times
3142 file.setOffset(m_totalNoSamples, 0);
3143 file.readArray(×Buffer[0], "time");
3144
3145 // Check the times with those of other surfaces already loaded on this rank (before overwriting them), later
3146 // cross-check the times with all surfaces
3149 TERMM_IF_COND(!approx(timesBuffer[m_sampleOffset + i * m_sampleStride], m_times[i], MFloatEps),
3150 "Error: local times for surfaces do not match.");
3151 }
3152 }
3153
3154 // Store times in m_times
3157 } else {
3160 }
3161 }
3162 } else {
3163 file.setOffset(0, 0);
3164 file.readArray(×Buffer[0], "time");
3165 }
3166
3167 MPI_Bcast(&m_times[0], noSamples(), type_traits<MFloat>::mpiType(), 0, m_mpiCommSurface[surfaceId], AT_,
3168 "&m_times[0]");
3170
3171 // Set 2D offset (second dimension is the time axis which is set internally to read all samples)
3172 // TODO labels:ACA,IO directly read only selected sample range once this is supported by the ParallelIO class
3173 file.setOffset(count, offset, 2);
3174
3175 // Data buffer
3177 // Read datasets and store in m_surfaceData
3180
3182 std::cout << "Surface #" << m_surfaceIds[surfaceId] << ": Loading variable: " << varName << std::endl;
3183 }
3184
3185 file.readArray(&varData[0], varName);
3186
3192 }
3193 }
3194 }
3195 } else { // ! m_useMergedInputFile -------------------------------------------
3198
3200 std::vector<MInt> sendCount(noDom);
3201 std::vector<MInt> sendCountNdim(noDom);
3202 std::vector<MInt> displs(noDom);
3203 std::vector<MInt> displsNdim(noDom);
3204 {
3205 std::fill(sendCount.begin(), sendCount.end(), 0);
3206 std::fill(displs.begin(), displs.end(), 0);
3207 sendCount[m_domainIdSurface[surfaceId]] = count;
3208 displs[m_domainIdSurface[surfaceId]] = offset;
3209 MPI_Allreduce(MPI_IN_PLACE, sendCount.data(), sendCount.size(), MPI_INT, MPI_SUM, m_mpiCommSurface[surfaceId],
3210 AT_, "MPI_IN_PLACE", "sendCount");
3211 MPI_Allreduce(MPI_IN_PLACE, displs.data(), displs.size(), MPI_INT, MPI_SUM, m_mpiCommSurface[surfaceId], AT_,
3212 "MPI_IN_PLACE", "displs");
3214 sendCountNdim[i] = sendCount[i] * nDim;
3215 displsNdim[i] = displs[i] * nDim;
3216 }
3217 }
3219
3220 //--Read sortIndex from first input file
3221 std::vector<MInt> sortIndex;
3222 {
3224 //-Open the correct sample file
3227 //-First read the famous sortIndex_ from first file
3230 sortIndex_.resize(noTotalCount);
3231 file.setOffset(noTotalCount, 0);
3232 file.readArray(sortIndex_.data(), varName);
3233 } else {
3234 sortIndex_.resize(0);
3235 file.setOffset(0, 0);
3236 file.readArray(sortIndex_.data(), varName);
3237 }
3238 };
3239 readSortIndex(m_inputFileIndexStart, sortIndex);
3240 // Sanity check: Is sort index the same for all input files?
3241 std::vector<MInt> tmpSortIndex(sortIndex.size());
3243 readSortIndex(fileNo, tmpSortIndex);
3244 if(!std::equal(sortIndex.begin(), sortIndex.end(), tmpSortIndex.begin())) {
3245 std::stringstream ss;
3246 ss << "SortIndex varies between ";
3248 ss << " and ";
3249 ss << getSurfaceDataFileName(surfaceId, fileNo);
3250 TERMM(1, ss.str());
3251 }
3252 }
3253 }
3254
3255 //--Read geometry data from separate file
3256 {
3260
3261 //-Read everything from first rank and sort it
3265 std::ostringstream fileName;
3266 fileName << m_inputDir << baseGeometryFileName << std::to_string(m_surfaceIds[surfaceId]) << baseFileNameSuffix;
3271 // Read surface area/length
3272 file.setOffset(noTotalCount, 0);
3273 file.readArray(readBuffer.data(), varName);
3275 sendBuffer(i) = readBuffer(sortIndex[i]);
3276 }
3277 MPI_Scatterv(sendBuffer.data(), sendCount.data(), displs.data(), MPI_DOUBLE, measure.data(), count, MPI_DOUBLE,
3279
3280 // Read coordinates and normal vectors
3281 file.setOffset(noTotalCount * nDim, 0);
3282 file.readArray(readBuffer.data(), "centroid");
3285 sendBuffer(i * nDim + j) = readBuffer(sortIndex[i] * nDim + j);
3286 }
3287 }
3288 MPI_Scatterv(sendBuffer.data(), sendCountNdim.data(), displsNdim.data(), MPI_DOUBLE, coordinates.data(),
3290 file.readArray(readBuffer.data(), "normalVector");
3293 sendBuffer(i * nDim + j) = readBuffer(sortIndex[i] * nDim + j);
3294 }
3295 }
3296 MPI_Scatterv(sendBuffer.data(), sendCountNdim.data(), displsNdim.data(), MPI_DOUBLE, normalVector.data(),
3298 //-Scatter data (send)
3299 } else { // not root rank -> receive sorted data
3300 //-Scatter data (receive)
3302 file.setOffset(0, 0);
3303 file.readArray(dummySendBuffer, varName);
3304 MPI_Scatterv(dummySendBuffer, sendCount.data(), displs.data(), MPI_DOUBLE, measure.data(), count, MPI_DOUBLE, 0,
3306 file.readArray(dummySendBuffer, "centroid");
3307 MPI_Scatterv(dummySendBuffer, sendCountNdim.data(), displsNdim.data(), MPI_DOUBLE, coordinates.data(),
3309 file.readArray(dummySendBuffer, "normalVector");
3310 MPI_Scatterv(dummySendBuffer, sendCountNdim.data(), displsNdim.data(), MPI_DOUBLE, normalVector.data(),
3312 }
3313
3314 // Create surface elements and store its data in m_surfaceData
3315 createSurfaceElements(surfaceId, measure, coordinates, normalVector);
3316 }
3317
3318 //--Read times only on surface local rank 0 and broadcast to other ranks in this communicator
3319 {
3321 // const MString varName = "timeStep";
3324 MInt offsetInBuffer = 0;
3326 // get fileName of fileNo th sample
3329
3330 std::vector<MLong> dims = file.getArrayDims(varName);
3332
3333 file.setOffset(samplesInFile, 0);
3334 file.readArray(×Buffer[offsetInBuffer], varName);
3335
3336 offsetInBuffer += samplesInFile;
3337 }
3338 // Check the times with those of other surfaces already loaded on this rank (before overwriting them), later
3339 // cross-check the times with all surfaces
3342 TERMM_IF_COND(!approx(timesBuffer[m_sampleOffset + i * m_sampleStride], m_times[i], MFloatEps),
3343 "Error: local times for surfaces do not match.");
3344 }
3345 }
3346
3347 // Store times in m_times
3350 } else {
3353 }
3354 }
3355 } else { // not root process
3357 // get fileName of fileNo th sample
3359 // read only dummy to ensure valid 'parallel' reading
3361 file.setOffset(0, 0);
3362 file.readArray(×Buffer[0], varName);
3363 }
3364 }
3365 MPI_Bcast(&m_times[0], noSamples(), type_traits<MFloat>::mpiType(), 0, m_mpiCommSurface[surfaceId], AT_,
3366 "&m_times[0]");
3368 }
3369
3370 //--Read data from each file and sort it correctly into my buffers
3371 {
3373 MInt offsetInBuffer = 0;
3375 //-Open the correct sample file
3378
3379 //-Read variable data
3380 {
3381 // Read data
3383 file.setOffset(count, offset, 3);
3384
3385 std::vector<MLong> dims = file.getArrayDims(varName);
3387
3389 file.readArray(varData.data(), varName);
3390
3391 // Create mapping from input variables in file to order of stored variable
3392 std::vector<MInt> varOrder;
3393 {
3394 std::map<MInt, MString> fileVars{};
3396 MString var;
3397 file.getAttribute(&var, "var_" + std::to_string(varIndex), varName);
3398 fileVars[varIndex] = var;
3399 }
3401 varOrder.push_back(m_inputVars[fileVars[varIndex]]);
3402 }
3403 }
3404
3405 // Store in local data
3411 // TODO labels:ACA add possibility to set a strideA
3412 // TODO labels:ACA add offset, e.g. skip first X samples in files
3414 }
3415 }
3416 }
3417 offsetInBuffer += samplesInFile;
3418 }
3419 }
3420 }
3421 }
3422}
MString getSurfaceDataFileName(const MInt surfaceId, const MInt fileNo)
IO methods.
Definition: acasolver.cpp:2913
int MPI_Scatterv(const void *sendbuf, const int sendcount[], const int displs[], MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, int root, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Scatterv
Definition: mpioverride.cpp:1074
◆ saveObserverSignals()
Definition at line 2388 of file acasolver.cpp.
2388 {
2389 TRACE();
2391
2393
2395 const MString fileName = outputDir() + m_outputFilePrefix + "observerData" + ParallelIo::fileExt();
2396
2398
2399 ParallelIo::size_type dimSizesCoordinates[] = {globalNoObservers(), nDim};
2400 file.defineArray(PIO_FLOAT, "coordinates", 2, &dimSizesCoordinates[0]);
2401
2404
2406 file.defineArray(PIO_FLOAT, "pressure_time", 2, &dimSizesVarTime[0]);
2407
2409 file.defineArray(PIO_FLOAT, "pressure_frequency", 2, &dimSizesVarFreq[0]);
2410
2411 // Store main parameters/configuration as attributes
2413 file.setAttribute(acousticMethod, "acousticMethod");
2417 }
2419 // Sample attributes
2424 // Windowing
2427 // Surface attributes
2429 file.setAttribute(m_surfaceInputFileName[sId], "inputFileName_" + std::to_string(m_surfaceIds[sId]));
2430 file.setAttribute(m_noSurfaceElementsGlobal[sId], "noSurfaceElements_" + std::to_string(m_surfaceIds[sId]));
2431 file.setAttribute(m_surfaceWeightingFactor[sId], "surfaceWeight_" + std::to_string(m_surfaceIds[sId]));
2432 }
2434
2435 // Store also observer coordinates to have everything in one file
2437 file.setOffset(globalNoObservers(), 0, 2);
2438 } else {
2439 file.setOffset(0, 0, 2);
2440 }
2442
2443 // write times and frequencies
2445 file.setOffset(noSamples(), 0);
2446 } else {
2447 file.setOffset(0, 0);
2448 }
2451
2453 ScratchSpace<MFloat> pressure_time(std::max(noObservers(), 1), noSamples(), AT_, "pressure_time");
2457 }
2458 }
2459 file.writeArray(&pressure_time[0], "pressure_time");
2460
2462 ScratchSpace<MFloat> pressure_frequency(std::max(noObservers(), 1), 2 * noSamples(), AT_, "pressure_frequency");
2467 }
2468 }
2469 file.writeArray(&pressure_frequency[0], "pressure_frequency");
2470
2472}
static const std::array< MString, 4 > windowNames
Definition: acasolver.h:544
◆ saveSolverSolution()
template<MInt nDim>
|
inlinevirtual |
Definition at line 117 of file acasolver.h.
117{};
◆ setInputOutputProperties()
Definition at line 110 of file acasolver.cpp.
◆ setNumericalProperties()
Definition at line 424 of file acasolver.cpp.
◆ solutionStep()
Reimplemented from Solver.
Definition at line 807 of file acasolver.cpp.
807 {
808 TRACE();
810
811 // Distribute surface elements among ranks and create separate communicators etc.
812 initParallelization();
813
814 initConversionFactors();
815
816 // read in data of all surfaces
817 // Note: for validating/testing instead of loading the data from disk you can generate the
818 // analytical input data
819 readSurfaceData();
820
821 // read in observer points or generate the according to some properties if specified
822 initObserverPoints();
823
825 // Compute source terms depending on method
826 calcSourceTerms();
827
828 // Compute observer pressure signal by using time or frequency FW-H formulation
830 // apply window to source terms and transform
831 windowAndTransformSources();
832
833 calcFrequencies();
834
835 // Calculate the observer signals in frequency domain
836 calculateObserverInFrequency();
837
838 // backtransformation of local observer signal from frequency into time domain
839 backtransformObservers();
841 // Surface Integration
842 calculateObserverInTime();
843
844 calcFrequencies();
845
846 // Windowing and FFT of the observer data
847 windowAndTransformObservers();
848 }
849
850 // change dimensions from ACA to desired output format
851 changeDimensionsObserverData();
852
853 // output of observer signals (frequency and time domain)
854 saveObserverSignals();
855
856 // perform additional post-processing operations with the observer signals, e.g., compute
857 // RMS-pressure, average over observers, ..., and output the data
858 postprocessObserverSignals();
859
860 // Calculates the accuracy, the difference between analytic and computed solution. Can be only be
861 // used with analytical data at the moment.
862 // computeAccuracy();
863 } else {
864 loadObserverSignals();
865
866 // Postprocessing with loaded observer data
867 postprocessObserverSignals();
868 }
869
871
873 return true;
874}
void loadObserverSignals()
Load the observer signals from a file for further postprocessing.
Definition: acasolver.cpp:2477
void initConversionFactors()
Initialize the conversion factors required to convert between.
Definition: acasolver.cpp:4721
void changeDimensionsObserverData()
Definition: acasolver.cpp:4816
void calculateObserverInFrequency()
Calculate the observer signals in frequency domain.
Definition: acasolver.cpp:2294
void windowAndTransformSources()
Apply window function to source terms and transform into frequency domain.
Definition: acasolver.cpp:1253
void saveObserverSignals()
Save the observer signals. Save observer data: 1. pressure data in time domain and 2....
Definition: acasolver.cpp:2388
void readSurfaceData()
Read the surface data (or generate for an analytical case)
Definition: acasolver.cpp:2925
void initParallelization()
Initialize parallelization (distribution of surface segments on all domains)
Definition: acasolver.cpp:2672
void backtransformObservers()
Backtransform the observer frequency signals into the time domain.
Definition: acasolver.cpp:2348
void windowAndTransformObservers()
Apply window function to the observer time series data and transform it into frequency domain.
Definition: acasolver.cpp:1354
void calculateObserverInTime()
Calculate the observer signals in time domain.
Definition: acasolver.cpp:2309
◆ solverConverged()
◆ storeSurfaceData()
Definition at line 3428 of file acasolver.cpp.
3428 {
3429 TRACE();
3432 + std::to_string(m_surfaceIds[surfaceId]) + ParallelIo::fileExt();
3433 std::stringstream ss;
3434 ss << " "
3435 << "Testing: Storing surface data in " << fileName;
3436 printMessage(ss.str());
3437
3439
3443
3445 file.defineArray(PIO_FLOAT, measureName, totalCount);
3446
3448
3449 ParallelIo::size_type dimSizesXD[] = {totalCount, nDim};
3450
3451 file.defineArray(PIO_FLOAT, "coordinates", 2, &dimSizesXD[0]);
3452 file.defineArray(PIO_FLOAT, "normalVector", 2, &dimSizesXD[0]);
3453
3454 ParallelIo::size_type dimSizesVar[] = {totalCount, noSamples()};
3455
3458 file.defineArray(PIO_FLOAT, varName, 2, &dimSizesVar[0]);
3459 }
3460
3462
3463 // Write surface area/length
3464 file.setOffset(count, offset);
3466
3467 // Write surface area/length
3470
3471 // Write coordinates and normal vectors
3472 file.setOffset(count, offset, 2);
3475
3476 file.setOffset(count, offset, 2);
3483 }
3484 }
3485
3486 file.writeArray(&varBuffer[0], var.first);
3487 }
3488}
◆ time()
◆ totalNoSurfaceElements()
Definition at line 138 of file acasolver.h.
◆ transformBackToOriginalCoordinate2D()
template<MInt nDim>
|
private |
Definition at line 2020 of file acasolver.cpp.
2020 {
2022
2023 // Back transform to the original coordinate system
2024 // Use negative angular position (-theta) in m_transformationMatrix calculation
2027
2030}
std::array< MFloat, nDim *nDim > m_transformationMatrix
Number of local observer.
Definition: acasolver.h:375
◆ transformBackToOriginalCoordinate3D()
template<MInt nDim>
|
private |
Definition at line 2033 of file acasolver.cpp.
2033 {
2035
2036 // Back transform to the original coordinate system
2037 // Use negative angular position (-alpha) and (-theta) in m_transformationMatrix calculation
2041
2042 *fX = fX_old * m_transformationMatrix[0] - fY_old * m_transformationMatrix[1] - fZ_old * m_transformationMatrix[2];
2043 *fY = -fX_old * m_transformationMatrix[3] + fY_old * m_transformationMatrix[4] + fZ_old * m_transformationMatrix[5];
2044 *fZ = -fX_old * m_transformationMatrix[6] + fY_old * m_transformationMatrix[7] + fZ_old * m_transformationMatrix[8];
2045}
◆ transformCoordinatesToFlowDirection2D()
template<MInt nDim>
|
private |
Definition at line 1994 of file acasolver.cpp.
1994 {
1996
1997 // Project the vector on the flow direction U_inf = (U_inf, 0)
2000
2003}
◆ transformCoordinatesToFlowDirection3D()
template<MInt nDim>
|
private |
Definition at line 2006 of file acasolver.cpp.
2006 {
2008
2009 // Project the vector on the flow direction U_inf = (U_inf, 0, 0)
2013
2014 *fX = fX_old * m_transformationMatrix[0] + fY_old * m_transformationMatrix[1] + fZ_old * m_transformationMatrix[2];
2015 *fY = fX_old * m_transformationMatrix[3] + fY_old * m_transformationMatrix[4] + fZ_old * m_transformationMatrix[5];
2016 *fZ = fX_old * m_transformationMatrix[6] + fY_old * m_transformationMatrix[7] + fZ_old * m_transformationMatrix[8];
2017}
◆ windowAndTransformObservers()
Definition at line 1354 of file acasolver.cpp.
1354 {
1355 TRACE();
1356
1357 // Compute the window function factors and the normalization factor
1358 std::vector<MFloat> window(noSamples());
1359 std::fill(window.begin(), window.end(), 1.0);
1360 calcWindowFactors(window.data());
1361
1364 // Sum up observer pressure over all samples
1365 MFloat sumObsPressure = 0.0;
1368 }
1370
1371 // Substract average from all samples to get the fluctuations
1372 // After substracting the mean apply the window function [see e.g. Mendez et al. or Lockard]
1373 // copy variables into complex variables array: re1, im1, re2, im2, re3, im3, ... with im_n = 0
1375 // Do nothing to "m_observerData.variables" output the original data
1379 }
1380
1381 // Use FFT or DFT to transform source terms (depending on possibility and the set property)
1382 checkNoSamplesPotencyOfTwo();
1384 FastFourierTransform(&m_observerData.complexVariables(observerId, 0), noSamples(), 1, nullptr, true);
1386 // 1 in function means forward; -1 is backward
1387 DiscreteFourierTransform(&m_observerData.complexVariables(observerId, 0), noSamples(), 1, nullptr, true);
1388 } else {
1389 TERMM(1, "Error Fourier-Transform: Check transformationType and noSamples.");
1390 }
1391
1392 // Normalize for conservation of energy
1396 }
1397 } // End of observer loop
1398}
void checkNoSamplesPotencyOfTwo()
Definition: acasolver.cpp:4864
void calcWindowFactors(MFloat *const p_window)
Compute the window function factors and the normalization factor.
Definition: acasolver.cpp:1195
◆ windowAndTransformSources()
Definition at line 1253 of file acasolver.cpp.
1253 {
1254 TRACE();
1257
1258 // Compute the window function factors and the normalization factor
1259 std::vector<MFloat> window(noSamples());
1260 std::fill(window.begin(), window.end(), 1.0);
1261 calcWindowFactors(window.data());
1262
1263 // Set source term indices depending on used method
1264 std::vector<MInt> VarDim;
1265 std::vector<MInt> VarDimC;
1267 if constexpr(nDim == 3) {
1270 } else {
1273 }
1275 if constexpr(nDim == 3) {
1278 } else {
1281 }
1282 } else {
1283 TERMM(1, "Unknown acoustic method.");
1284 }
1285
1286 MInt noSourceTerms = 0;
1287 if constexpr(nDim == 3) {
1288 noSourceTerms = 4;
1289 } else if constexpr(nDim == 2) {
1290 noSourceTerms = 3;
1291 }
1292
1296 // Sum up source term over all samples
1297 MFloat sumSource = 0.0;
1300 }
1302
1303 // Substract average from all samples to get the fluctuations
1304 // After substracting the mean apply the window function [see e.g. Mendez et al. or Lockard]
1308 }
1309 }
1310
1311 // Copy variables into complex variables array: Re1, Im1, Re2, Im2, Re3, Im3, ... with Im_n = 0
1317 }
1318 }
1319
1320 // Use FFT or DFT to transform source terms (depending on possibility and the set property)
1321 checkNoSamplesPotencyOfTwo();
1324 FastFourierTransform(&m_surfaceData.complexVariables(segmentId, VarDimC[i]), noSamples(), 1, nullptr, true);
1325 }
1328 // 1 in function means forward; -1 is backward
1329 DiscreteFourierTransform(&m_surfaceData.complexVariables(segmentId, VarDimC[i]), noSamples(), 1, nullptr, true);
1330 }
1331 } else {
1332 TERMM(1, "Error Fourier-Transform: Check transformationType and noSamples.");
1333 }
1334
1335 // Normalize for conservation of energy
1340 }
1341 }
1342
1343 // Calculate frequencies is moved out of segmentId Loop because it is the same for every
1344 // segmentId. Moved to calcSurfaceIntegral Calculate time: It is calculated at
1345 // "generateSurfaceData" m_times = dt * noSamples() Initialize observerdata array. Function
1346 // initObserverPoints
1347 } // End of loop over all segments
1348
1350}
Member Data Documentation
◆ m_aca2Output
template<MInt nDim>
|
private |
◆ m_acaPostprocessingFile
Definition at line 442 of file acasolver.h.
◆ m_acaPostprocessingMode
Definition at line 441 of file acasolver.h.
◆ m_acaResultsNondimMode
Definition at line 288 of file acasolver.h.
◆ m_acousticMethod
Definition at line 347 of file acasolver.h.
◆ m_allowMultipleSurfacesPerRank
Definition at line 313 of file acasolver.h.
◆ m_Analyticfreq
Definition at line 412 of file acasolver.h.
◆ m_cInfty
Definition at line 408 of file acasolver.h.
◆ m_domainIdSurface
Definition at line 326 of file acasolver.h.
◆ m_dt
Definition at line 418 of file acasolver.h.
◆ m_FastFourier
Definition at line 430 of file acasolver.h.
◆ m_frequencies
Definition at line 434 of file acasolver.h.
◆ m_fwhMeanStateType
Definition at line 349 of file acasolver.h.
◆ m_fwhTimeShiftType
Definition at line 348 of file acasolver.h.
◆ m_gamma
Definition at line 427 of file acasolver.h.
◆ m_generateObservers
Definition at line 388 of file acasolver.h.
◆ m_generateSurfaceData
Definition at line 392 of file acasolver.h.
◆ m_globalObserverCoordinates
Definition at line 368 of file acasolver.h.
◆ m_hasTimes
Definition at line 435 of file acasolver.h.
◆ m_input2Aca
template<MInt nDim>
|
private |
◆ m_inputDir
Definition at line 379 of file acasolver.h.
◆ m_inputFileIndexEnd
Definition at line 307 of file acasolver.h.
◆ m_inputFileIndexStart
Definition at line 306 of file acasolver.h.
◆ m_inputVars
Definition at line 355 of file acasolver.h.
◆ m_isFinished
Member variables Solver status
Definition at line 286 of file acasolver.h.
◆ m_lambdaMach
Definition at line 405 of file acasolver.h.
◆ m_lambdaZero
Definition at line 404 of file acasolver.h.
◆ m_Ma
Definition at line 423 of file acasolver.h.
◆ m_MaDim
Definition at line 425 of file acasolver.h.
◆ m_MaVec
Definition at line 421 of file acasolver.h.
◆ m_mpiCommSurface
Definition at line 322 of file acasolver.h.
◆ m_noComplexVariables
Definition at line 352 of file acasolver.h.
◆ m_noDomainsSurface
Definition at line 324 of file acasolver.h.
◆ m_noGlobalObservers
Definition at line 366 of file acasolver.h.
◆ m_noObservers
Definition at line 372 of file acasolver.h.
◆ m_noPostprocessingOps
Definition at line 438 of file acasolver.h.
◆ m_noSamples
Definition at line 335 of file acasolver.h.
◆ m_noSurfaceElements
Definition at line 310 of file acasolver.h.
◆ m_noSurfaceElementsGlobal
Definition at line 316 of file acasolver.h.
◆ m_noSurfaces
Definition at line 301 of file acasolver.h.
◆ m_noVariables
Definition at line 351 of file acasolver.h.
◆ m_observerData
template<MInt nDim>
|
private |
Definition at line 369 of file acasolver.h.
◆ m_observerFileName
Definition at line 386 of file acasolver.h.
◆ m_offsetObserver
Definition at line 371 of file acasolver.h.
◆ m_omega_factor
Definition at line 401 of file acasolver.h.
◆ m_outputFilePrefix
Definition at line 380 of file acasolver.h.
◆ m_post
template<MInt nDim>
|
private |
Definition at line 272 of file acasolver.h.
◆ m_postprocessingOps
Definition at line 439 of file acasolver.h.
◆ m_rhoInfty
Definition at line 409 of file acasolver.h.
◆ m_rmsP_Analytic
Definition at line 415 of file acasolver.h.
◆ m_sampleOffset
Definition at line 341 of file acasolver.h.
◆ m_sampleStride
Definition at line 339 of file acasolver.h.
◆ m_sourceParameters
template<MInt nDim>
|
private |
Definition at line 398 of file acasolver.h.
◆ m_sourceType
Definition at line 395 of file acasolver.h.
◆ m_storeSurfaceData
Enable storing of surface data (testing or output of generated data)
Definition at line 383 of file acasolver.h.
◆ m_surfaceData
template<MInt nDim>
|
private |
Definition at line 344 of file acasolver.h.
◆ m_surfaceElementOffsets
Definition at line 319 of file acasolver.h.
◆ m_surfaceIds
Definition at line 303 of file acasolver.h.
◆ m_surfaceInputFileName
Definition at line 332 of file acasolver.h.
◆ m_surfaceWeightingFactor
Definition at line 329 of file acasolver.h.
◆ m_symBc
Definition at line 558 of file acasolver.h.
◆ m_timerGroup
Definition at line 446 of file acasolver.h.
◆ m_timers
template<MInt nDim>
|
private |
Definition at line 448 of file acasolver.h.
◆ m_times
Definition at line 433 of file acasolver.h.
◆ m_totalNoSamples
Definition at line 337 of file acasolver.h.
◆ m_transformationMatrix
template<MInt nDim>
|
private |
Coordinate system transformation
Definition at line 375 of file acasolver.h.
◆ m_transformationType
Definition at line 363 of file acasolver.h.
◆ m_useMergedInputFile
Definition at line 305 of file acasolver.h.
◆ m_windowNormFactor
Definition at line 360 of file acasolver.h.
◆ m_windowType
Definition at line 358 of file acasolver.h.
◆ m_zeroFlowAngle
Definition at line 376 of file acasolver.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/ACA/acasolver.h
- /home/gitlab-runner/scratch/builds/oxpnswJ6/1/aia/m-AIA/m-AIA/src/ACA/acasolver.cpp
