problems.h

Go to the documentation of this file.

/* -----------------------------------------------------------------------------------------
  $Id: problems.h 645 2007-05-11 13:44:31Z svacek $
  $Author: svacek $
   ----------------------------------------------------------------------------------------- */
#include "math.h"
#include "umfpack/umfpack.h"
#include "fluidf.h"
#include "vector.h"
#include "matrix.h"
#include "sparse.h"
#include "datastr.h"
#ifndef __PROBLEMS_SVH
#define __PROBLEMS_SVH
#define PAR_BNORM 0
#define PAR_H 1
#define PAR_REYNOLDS 2
#define PAR_XI 3
#define PAR_DELTA 4
/* --------------------------------------------------------------------------- */
double TimeDifference(double U1,double U2,double U3,double U4);
/* --------------------------------------------------------------------------- */
void INSEVelocitiesPreparation(long idx,const matrix M,motion U1,const vector u1,motion U0,const vector u0,motion W,const vector w);
void RANSVelocitiesPreparation(long idx,const matrix M,motion U1,const vector u1,motion U0,const vector u0,motion W,const vector w,basefunction NU_TURBULENT,const vector nu_turbulent);
void VariablesPreparation(long idx,const matrix M,motion W,const vector w,basefunction VARIABLE1,const vector variable1);
void VariablesPreparation2(long idx,const matrix M,motion W,const vector w,basefunction VAR1,const vector var1,basefunction VAR2,const vector var2);
void VariablesPreparation4(long idx,const matrix M,motion W,const vector w,basefunction VAR1,const vector var1,basefunction VAR2,const vector var2,basefunction VAR3,const vector var3,basefunction VAR4,const vector var4);
/* --------------------------------------------------------------------------- */
void   ConstructWind(vector *wind,double xwind(double,double),double ywind(double,double));
double ComputeQComponent(long idx,const matrix M,basefunction U,const vector _w);
double ComputeQLinComponent(long idx,const matrix M,basefunction U,const vector _w);
void   PrepareLinearFunctionValues(long _iELE,const matrix M,const double *GlobalNodeValues,double *LinearFunctionQValues);
void   ComputeQMotionNoP(long idx,const matrix M,motion U,const vector _w);
void   ComputeQMotion(long idx,const matrix M,motion U,const vector _w);
void   PrepareBFValues(long _bfn_no,const matrix M,basefunction H);
double ComputeDomain();
/* --------------------------------------------------------------------------- */
void ConstructMatrices(sparse *A,sparse *B,sparse *BT);
void PrepareProblem(sparse *A,vector b);
void PrepareProblemStable(sparse *A,vector b,const vector wind,double _fce(double,double));
void StokesProblem(sparse *S,vector b,double visc);
void OseenProblemStable(sparse *S,vector b,const vector wind,double _fce(double,double));
void NS_SemiImplicit(sparse *S,vector b,double visc,double tau,const vector u1,const vector u0,const vector w);
void NS_CrankNicholson(sparse *S,vector b,double visc,double tau,const vector u,const vector sol);
void ConvectionDiffusionReactionProblem(sparse *A,vector b,double _visc,double _react,const vector SOLaprox,const vector wind,double _fce(double,double));
void ConvectionDiffusionReactionArtificialDiffusion(sparse *A,vector b,double _visc,double _react,const vector wind,double _fce(double,double),int TYP);
/* --------------------------------------------------------------------------- */
/* Drag & lift evaluation                                                      */
/* --------------------------------------------------------------------------- */
void AerodynamicalHistoryExtrapolation();
void InviscidForcesEval();
void AerodynamicalForcesEval();
void AerodynamicalForces();
/* ------------------------------------------------------------------------------------- */
/* stabilizace */
double  StabFunctionXi(double _PecletNumber);
double  StabParameterH(long iELE,const matrix M,double _b1,double _b2,double _bnorm);
double StabParameterH_min(long iELE);
double StabParameterH_max(long iELE);
double StabParameterH_mean(long iELE);
int IsIsotropicElement(long iELE);
double  ComputeLocalStreamNorm(const double *_u,const double *_v,double *_b1,double *_b2);
double  ComputeStabParam(long iELE,double visc,double tau,const matrix M,const double *_u,const double *_v);
double  ComputeStabParamDiv(long iELE,double visc,double tau,const matrix M,const double *_u,const double *_v,double *TAU_K);
double ComputeStabParamTwoPhase(long iELE,double visc,double rho,double tau,const matrix M,const double *_u,const double *_v,double *TAU_K);
double  ComputeStabParamTurbulence(long iELE,double visc,double tau,double _react,const matrix M,const double *_u,const double *_v);
void ComputeStabParamADC(long iELE,double _visc,double _react,double _rez,double _grad,const matrix M,const double *_u,const double *_v, const double GRADX, const double GRADY,double *OUT_DELTA_K, double *OUT_DELTA_ADC);
/* ------------------------------------------------------------------------------------ */
/* ALE - method (vynucene kmitani)                                                      */
/* ------------------------------------------------------------------------------------ */
double DivOfGridVelocity(const long i);
double DivOfGridVelocity0(const long i);
double DivOfGridVelocity1(const long i);
double DivOfGridVelocity2(const long i);
double DivOfGridVelocity3(const long i);
double DivOfGridVelocity4(const long i);
void   ReinitPosition();
double VibrationsIn_ALPHA_H(double _t);
double VibrationsIn_ALPHA_H_BETA(double _t);
void UpdateGridAnalytical();
void MoveInitialGrid();
void MoveGrid_AlphaH(double _Alpha,double _hh);
void UpdateGrid_AlphaH(double _Alpha,double _hh);
void UpdateGridVelocity(vector unew);
void SubtractLocalDomainVelocity(long iELE,motion W);
void InitGridElasticExtension();
void MoveGridElasticExtension();
void UpdateGridElasticExtension();
void MoveAirfoilBoundaryWithFlap(const double WingAngle,const double Vertical,const double TailRelAngle);
void MoveAirfoilBoundary(const double WingAngle,const double Vertical);
/* ------------------------------------------------------------------------------------ */
/** @brief
 * NS_NONNEWTONIAN: 
 * Discretization of Non-Newtonian NS equations with 3-point 2nd order time difference. Stabilization with SUPG method.
 * @param[out] S    linear system of eqs.
 * @param[out] b    right hand side of the system.
 * @param[in] visc  nondimensional viscosity.
 * @param[in] tau   time step.
 * @param[in] u2    vector with UVP at T - approximation.
 * @param[in] u1    vector with UVP at T-  tau.
 * @param[in] u0    vector with UVP at T-2 tau.
 * @param[in] ww    convective velocity.
*/
void NS_NONNEWTONIAN(sparse *S,vector b,double visc,double tau,const vector u2,const vector u1,const vector u0,const vector ww);
void NS_SUPG1(sparse *S,vector b,double visc,double tau,const vector u1,const vector ww);
/* ------------------------------------------------------------------------------------ */
/** @brief
 * NS_SUPG: 
 * Discretization of NS equations with 3-point 2nd order time difference. Stabilization with SUPG method.
 * @param[out] S    linear system of eqs.
 * @param[out] b    right hand side of the system.
 * @param[in] visc  nondimensional viscosity.
 * @param[in] tau   time step.
 * @param[in] u1    unknowns vector with UVP at T-  tau.
 * @param[in] u0    unknowns vector with UVP at T-2 tau.
 * @param[in] ww    convective velocity.
*/
void NS_SUPG(sparse *S,vector b,double visc,double tau,const vector u1,const vector u0,const vector ww);
/* ------------------------------------------------------------------------------------ */
/** @brief
 * NS_GALS:
 * Discretization of NS equations with 3-point 2nd order time difference. Stabilization with GaLS method.
 * @param[out] S    linear system of eqs.
 * @param[out] b    right hand side of the system.
 * @param[in] visc  nondimensional viscosity.
 * @param[in] tau   time step.
 * @param[in] u1    unknowns vector with UVP at T-  tau.
 * @param[in] u0    unknowns vector with UVP at T-2 tau.
 * @param[in] ww    convective velocity.
*/
void NS_GALS(sparse *S,vector b,double visc,double tau,const vector u1,const vector u0,const vector ww);
/* ------------------------------------------------------------------------------------ */
/** @brief
 * NS_SUBGRID:
 * Discretization of NS equations with 3-point 2nd order time difference. Stabilization with SUBGRID-SCALE method.
 * @param[out] S    linear system of eqs.
 * @param[out] b    right hand side of the system.
 * @param[in] visc  nondimensional viscosity.
 * @param[in] tau   time step.
 * @param[in] u1    unknowns vector with UVP at T-  tau.
 * @param[in] u0    unknowns vector with UVP at T-2 tau.
 * @param[in] ww    convective velocity.
*/
void NS_SUBGRID(sparse *S,vector b,double visc,double tau,const vector u1,const vector u0,const vector ww);
/* ------------------------------------------------------------------------------------ */
/** @brief
 * NS_MODGLS: 
 * Discretization of NS equations with 3-point 2nd order time difference. Stabilization with Wang GLS method.
 * @param[out] S    linear system of eqs.
 * @param[out] b    right hand side of the system.
 * @param[in] visc  nondimensional viscosity.
 * @param[in] tau   time step.
 * @param[in] u1    unknowns vector with UVP at T-  tau.
 * @param[in] u0    unknowns vector with UVP at T-2 tau.
 * @param[in] ww    convective velocity.
*/
void NS_MODGLS(sparse *S,vector b,double visc,double tau,const vector u1,const vector u0,const vector ww);
void NS_PMGRID_Coarse(sparse *S,vector b,double visc,double tau,const vector u1,const vector u0,const vector ww,const vector uLES);
void NS_PMGRID_Fine(sparse *S,vector b,double visc,double tau,const vector u1,const vector u0,const vector ww,const vector uLES);
void DiscretizeNSSystem(sparse *S,vector b,double visc,double tau,const vector u1,const vector u0,const vector ww);
/* ------------------------------------------------------------------------------------ */
double ConservativityTest(double FluidModel(const long,const long,const double),void FunctionGridMotion());
double AbstractFluidModel(const long iFSI,const long _MaxOfIter,const double _eps);
double AbstractStructureModel(const long iFSI);
double NonNewtonianFluidModel(const long iFSI,const long _NumberOfIter,const double _eps);
double LaminarFluidModel(const long iFSI,const long _NumberOfIter,const double _eps);
double SolveFSI_TimeStep(double FluidModel(const long,const long,const double),double StructuralModel(const long));
double SolveLaminar(const long _NumberOfIter,const double _eps);
double SolveNonNewtonian(const long _NumberOfIter,const double _eps);
double SolvePMGRID(const long _NumberOfIter,const double _eps);
double SolveInviscid(const long _NumberOfIter,const double _eps);
double InviscidFluidModel(const long iFSI,const long _NumberOfIter,const double _eps);
/* ------------------------------------------------------------------------------------- */
double UpdateLastStepODE(double *_x,const int _NN);
void   ShiftLastStepsODE(double *_x,const int _NN);
void ODE_RHS_Nonlinear(double _t,double *_f,const double *_x,double CL,double CM);
void ODE_RHS_Linear(double _t,double *_f,const double *_x,double CL,double CM);
void ODE_RHS2(double _t,double *_f,const double *_x,double CL,double CM);
double ODEPredictorRK4(double _t);
double ODECorrectorRK4(double _t);
double ODE3_PredictorRK4(double _t);
double ODE3_CorrectorRK4(double _t);
void   ODE3_RHS_NonlinearForces(double _t,double *_f,const double *_x,double CL,double CM,double CFlap);
void   ODE3_RHS_LinearForces(double _t,double *_f,const double *_x,double CL,double CM,double CFlap);
void   ODE3_RHS(double _t,double *_f,const double *_x,double CL,double CM,double CFlap);
void   AerodynCoeficients(double *MLift,double *MTMoment,double *MLiftSlope,double *MTMomentSlope);
void   AerodynCoeficients3(double *MLift,double *MTMoment,double *MLiftSlope,double *MTMomentSlope,double *FlapM,double *FlapMSlope);
double StructureModel_Fixed(const long ifsi);
double StructureModel_Prescribed2DOF(const long ifsi);
double StructureModel_Prescribed3DOF(const long ifsi);
double StructureModel_Aeroelastic2DOF(const long ifsi);
double StructureModel_Aeroelastic3DOF(const long ifsi);
/* ------------------------------------------------------------------------------------ */
void MultiplyDD(vector out,const vector x);
void PrecondDD(vector out);
void ConstructPC(const sparse S);
void DestructPC(); 
/* ------------------------------------------------------------------------------------ */
void Project2HDIV(vector ww);
void HDivProjection(sparse *S,vector b,const vector ww);
/* ------------------------------------------------------------------------------------ */
void AddDoNothingBC(const vector ww);
void AddDonothingBCWithPressure(const vector ww,vector rhs);
void AddDonothingBCWithPressureCoarse(const vector ww,vector rhs);
void AddDonothingBCWithPressureFine(const vector ww,vector rhs);
void AddInviscidPressure();
void AddDonothingBCDensity(const vector ww,vector rhs);
void NormalVelocityZero();
/* ------------------------------------------------------------------------------------ */
void Project2HDIV(vector ww);
void HDivProjection(sparse *S,vector b,const vector ww);
void ComputeLinearProduction(const vector ww);
void ComputeLinearProductionPatches(const vector ww);
void ComputeLinearProductionPatches1(const vector ww);
void ComputeLinearProductionPatchesTURB(const vector ww,const vector varA, const vector varB);
void UpdateDensityViscosity();
void UpdateDensity();
void MeanValueReconstruction(vector _FEfunction,const vector _FVfunction);
/* ------------------------------------------------------------------------------------ */
double ComputeStabParamTurbulence(long iELE,double visc,double tau,double _react,const matrix M,const double *_u,const double *_v);
double ComputeStabParamTurbulenceSST_OMG(long iELE,double _visc,double _tau,double _react,const matrix M,const double *_u,const double *_v);
double ComputeStabParamTurbulenceSST_K(long iELE,double _visc,double _tau,double _react,const matrix M,const double *_u,const double *_v);
void   ComputeStabParamADC(long iELE,double _visc,double _react,double _rez,double _grad,const matrix M,const double *_u,const double *_v,const double GRADX, const double GRADY,double *OUT_DELTA_K, double *OUT_DELTA_ADC);
/* ------------------------------------------------------------------------------------ */
void DynamicElasticityFillLinearEqs(double Lame_Mu_Const,double Lame_Lambda_Const,long idxELE,double _volume,long idxJ,short sign_j,vector b);
void DynamicElasticityFillBodyForces(long _iELE,double _volume,long idxJ,short sign_j,double *_QVAL1,double *_QVAL2);
void DynamicLinearElasticity(sparse *S,vector b,double YOUNGS_MODULE,double POISSONS_RATIO);
void DynamicMassMatrix(sparse *M,const double _density);
void NewmarkMethodStiffnessMatrix(sparse *S,const double YOUNGS_MODULE,const double POISSONS_RATIO, const double density, double tau);
void NewmarkMethodDynamicLoads(vector b,const double YOUNGS_MODULE,const double POISSONS_RATIO, const double density,const double tau,const vector XVELOCITY,const vector XACCELERATION);
/* ------------------------------------------------------------------------------------ */
void StaticElasticityFillLinearEqs(double Lame_Mu_Const,double Lame_Lambda_Const,long idxELE,double _volume,long idxJ,short sign_j,vector b);
void StaticElasticityFillBodyForces(long _iELE,double _volume,long idxJ,short sign_j,double *_QVAL1,double *_QVAL2);
void StaticLinearElasticity(sparse *S,vector b,double YOUNGS_MODULE,double POISSONS_RATIO);
void StaticMassMatrix(sparse *M,const double _density);
/* ------------------------------------------------------------------------------------ */
void MeshMotionFillLinearEqs(double Lame_Mu_Const,double Lame_Lambda_Const,long idxELE,double _volume,long idxJ,short sign_j);
void ExtendALEVelocityFillStiffnessMatrix(sparse *S);
void ExtendALEVelocityFillTension(vector b);
/* ------------------------------------------------------------------------------------ */
void InitGridElasticExtension();
void MoveGridElasticExtension();
void UpdateGridElasticExtension();
void NS_SUPG_reziduum(vector rez,double visc,double tau,const vector u1,const vector u0,const vector ww);
void NS_GLS_reziduum(vector rez,const vector u2,const vector u1,const vector u0,const vector ww);
void NS_GLS_matrixmultiplication(vector mult,const vector u2,const vector u1,const vector u0,const vector ww);
void NS_EULER(sparse *S,vector b,double tau,const vector u1,const vector u0,const vector ww);
/* ------------------------------------------------------------------------------------ */
void SolveTwoPhaseFlow();
void InitTwoPhase();
void StepTwoPhaseFlow();
void LevelSetReinitializationA();
void LevelSetReinitialization();
void SurfaceIntersection(long iELE,double *xSEG,double *ySEG);
void LevelSetReconstruction(const vector _LevelSetFunction);
void LevelSetReconstructionA(const vector _LevelSetFunction);
void Density_ViscosityReconstruction(vector _P1Density,vector _P1Viscosity,const vector _P1LevelSetFunction);
double TWOPHASE_TransportExplicit(double deltatau,const vector ww);
void NS_TWOPHASE(sparse *S,vector b,double tau,const vector u1,const vector u0,const vector ww);
double LevelSetFunctionEvaluate(long iELE,int mQuad);
double LevelSetFunctionEvaluate1(long iELE,int mQuad);
void TWOPHASE_RHSVelocityNS(long iELE,double tau,double _rhoA,const double DELTA_K,double *QVAL_U,double *QVAL_V);
void TWOPHASE_FillMomentumEqs(long idxELE,double tau,double _rhoA,double _volume,long idxJ,short sign_j,double DELTA_K,double TAU_K);
void TWOPHASE_FillPressureGradNS(long idxELE,double _rhoA,double _volume,long idxJ,short sign_j,const double DELTA_K);
void TWOPHASE_FillContinuityPressureStabilization(long idxELE,double _volume,long idxJ,short sign_j,const double DELTA_K);
double TWOPHASE_RHSPressureNS(long iELE,double tau);
double TWOPHASE_TransportFEMImplicit(double deltatau,const vector u2,const vector u1);

void ComputeSurfaceHeight(double *AverageHeight,double *LowestHeight,double *HighestHeight,double *LeftSrf,double *RightSrf);
void ComputeEnergies();





#endif

---