/*==================================================================================== ! Purpose(s): ! ! Contain(s): cstress ! ! Author(s): Javier Lopez Lara (jav.lopez@unican.es) ! Gabriel Barajas (barajasg@unican.es) !====================================================================================*/ #include "CArrayMasks.h" #include "CVars.h" #include "CMisc.h" #include "CErrors.h" #include "math.h" #include #include #include "CMathMisc.h" #include #include "CCalcFricVel.h" #ifdef __INTEL_COMPILER #define PYM cvtilde #else #define PYM cvtilde_ #endif extern "C" void _FORTRAN PYM (_DOUBLE_ARRAY AAC, _DOUBLE_ARRAY AAR, _DOUBLE_ARRAY AAT, _DOUBLE_ARRAY AXCoefs, _DOUBLE_ARRAY AYCoefs, _DOUBLE_ARRAY AXi, _DOUBLE_ARRAY AYj, _DOUBLE_ARRAY ADelX, _DOUBLE_ARRAY ADelY, _DOUBLE_ARRAY AU, _DOUBLE_ARRAY AV, _DOUBLE_ARRAY AF, _DOUBLE_ARRAY ATauxx, _DOUBLE_ARRAY ATauxy, _DOUBLE_ARRAY ATauyy, _INT_ARRAY ANf, _INT_ARRAY ANpc, _INT_ARRAY Axnpc, EPorosityModel _ENUM APorosityModel, _DOUBLE_ARRAY AXPorosity, _DOUBLE_ARRAY AD50, _DOUBLE_ARRAY AGc, _DOUBLE_ARRAY Axa, _DOUBLE_ARRAY Axxb, _DOUBLE Axxt, _DOUBLE AGx, _DOUBLE AGy, _DOUBLE ADelT, _DOUBLE AMaxVelUComp, _DOUBLE AMaxVelVComp, _BOOL ASetVelToZero, _DOUBLE Afrsurf, _DOUBLE ARhoF, _DOUBLE ANu, _DOUBLE Asd, _INT Anks, _INT nrough, _INT AKb, _INT AKt, _INT AKl, _INT AKr, _INT AIMax, _INT AJMax, _INT AIdxLSp, _INT AIdxRSp) { /*!======================================================================= ! Purpose(s): ! ! Author(s): Javier Lopez Lara (jav.lopez@unican.es) ! Gabriel Barajas (barajasg@unican.es) !=======================================================================*/ static double *SpongeCoef=NULL; double newU, newV; int TotCells=AIMax*AJMax; bool KbRigid = (AKb==bcRigidFreeSlip || AKb==bcRigidNoSlip); bool KtRigid = (AKt==bcRigidFreeSlip || AKt==bcRigidNoSlip); bool KlRigid = (AKl==bcRigidFreeSlip || AKl==bcRigidNoSlip); bool KrRigid = (AKr==bcRigidFreeSlip || AKr==bcRigidNoSlip); // Roughness. double nks=Anks*1.0; double Ks=nks*Asd; bool IsRoughSurface=(nrough==1); // Computes the sponge coefficient. It just does this the first time (Actually, this can be done elsewhere). if (SpongeCoef==NULL) // If array memory hasn't been allocated yet... { SpongeCoef=(double *)malloc(AIMax*sizeof(double)); for (int i=0;i begining of the // sponge, 1.0 -> End). double SpX=0.0; if (AIdxLSp>0 && XC0 && XC>RSpX) SpX=(XC-RSpX)/RSpWidth; SpongeCoef[i]=1.0; if (SpX>0.0) { double spbb=0.4; double SpX3=SpX*SpX*SpX; SpongeCoef[i]=1.0-(SpX3*(spbb+(1.0-spbb)*SpX3)); } //if } //for } //if // Creates the mesh object. TMesh Mesh(AIMax, AJMax, AXCoefs, AYCoefs, AXi, AYj, ADelX, ADelY); // Creates the masks. TCMask wAC (Mesh, AAC); TRMask wAR (Mesh, AAR); TTMask wAT (Mesh, AAT); TCMask wF (Mesh, AF); TRMask wU (Mesh, AU); TTMask wV (Mesh, AV); TCMask wTauxx (Mesh, ATauxx); TTRMask wTauxy (Mesh, ATauxy); TCMask wTauyy (Mesh, ATauyy); TIntMask wNf (Mesh, ANf); TIntMask wPorIdx (Mesh, ANpc); TIntMask wXPorIdx (Mesh, Axnpc); Mesh.DefineArea(2, AIMax-1, 2, AJMax-1); Mesh.First2(); do { // Finds the nodes that will be considered in various derivatives. int NotEmptyCells=wNf.NodesNotEqualTo(ctEmptyCell); // Calculates which cells have porous media. int PorousCells; if (APorosityModel!=pmNone) PorousCells=wPorIdx.NodesNotEqualTo(1); else PorousCells=0; if (Mesh.IsInDefinedArea()) // If the cell is inside the area of calculus... { double SpCoef=SpongeCoef[Mesh.I-1]; double TauxyxT = T_Dx(wTauxy, NotEmptyCells); double TauxyyR = R_Dy(wTauxy, NotEmptyCells); double TauyyyT = T_Dy(wTauyy, NotEmptyCells); double TauxxxR = R_Dx(wTauxx, NotEmptyCells); double CDist, TRDist; TVector ObsSrfce; bool IsObsIntrfce=IsObstacleInterface(Mesh, wAC, wAR, wAT, KlRigid, KrRigid, KtRigid, KbRigid, ObsSrfce, CDist, TRDist); if (IsObsIntrfce) { assert(ObsSrfce.x*ObsSrfce.y==0.0); // Checks that ObsSrfce is parallel either to the X or Y axis. TVector VelTR(wU.ITR(), wV.ITR()); TVector ObsSrfceTan=NormalVector(ObsSrfce); double VTR=DotProduct(ObsSrfceTan, VelTR); double FricVel=CalcFrictionVelocity(E, VTR, TRDist, ANu, Kappa, IsRoughSurface, Ks); double SqrFricVel=FricVel*fabs(FricVel)*(ObsSrfce.x+ObsSrfce.y); // Trick: (ObsSrfce.x+ObsSrfce.y) is the sign of the meaningful component of the vector. if (ObsSrfce.y>0.0) TauxyyR=(wTauxy.c()-SqrFricVel)/(2*TRDist); else if (ObsSrfce.y<0.0) TauxyyR=(SqrFricVel-wTauxy.b())/(2*TRDist); else if (ObsSrfce.x>0.0) TauxyxT=(wTauxy.c()-SqrFricVel)/(2*TRDist); else if (ObsSrfce.x<0.0) TauxyxT=(SqrFricVel-wTauxy.l())/(2*TRDist); } double VisX = (TauxxxR+TauxyyR)/ARhoF; // <- LMask.r() = RMask.c() double VisY = (TauyyyT+TauxyxT)/ARhoF; // <- BMask.t() = TMask.c() // VTilde is calculated in a different way wether the cell is a porous one or not. int PorIdxC=wXPorIdx.c(); if (PorousCells==0) // If the cell has not porous material at any of the sides... { newU = SpCoef*wU.c()+ADelT*(VisX+AGx); newV = SpCoef*wV.c()+ADelT*(VisY+AGy); } else // If the cell has porous material at any of the sides... { // Gets the porosity index at every side of the cell. int PIdxL, PIdxR, PIdxT, PIdxB; IdxUnpack(PorIdxC, PIdxL, PIdxR, PIdxT, PIdxB); // Linear and nonlinear stress for porous media double kcR = max2(fabs(wU.c()),0.01L)*max2(1.0L,Axxt)/(AXPorosity[PIdxR]*AD50[PIdxR]); double bbR = Axxb[PIdxR]*(1.0+7.5/kcR); double kcT = max2(fabs(wV.c()),0.01L)*max2(1.0L,Axxt)/(AXPorosity[PIdxT]*AD50[PIdxT]); double bbT = Axxb[PIdxT]*(1.0+7.5/kcT); double ModUC=sqrt(wU.c()*wU.c()+wV.c()*wV.c()); newU = SpCoef*wU.c()+ADelT*(VisX+AGx-(Axa[PIdxR]+bbR*ModUC)*wU.c())/AGc[PIdxR]; newV = SpCoef*wV.c()+ADelT*(VisY+AGy-(Axa[PIdxT]+bbT*ModUC)*wV.c())/AGc[PIdxT]; //newU = SpCoef*wU.c()+ADelT*(VisX+AGx)/AGc[PIdxR]; //newV = SpCoef*wV.c()+ADelT*(VisY+AGy)/AGc[PIdxT]; } //if // This code is a patch to 'solve' the problem that occur when small drops are over a // horizontal obstacle. In these cases, velocities grow due to a not completely known // effect related to the advection mechanism. if (fabs(newU)>AMaxVelUComp) { //printf("U too big at (%d,%d)", Mesh.I, Mesh.J); if (ASetVelToZero) newU=0.0; else newU=copysign(AMaxVelUComp, newU); } if (fabs(newV)>AMaxVelVComp) { //printf("V too big at (%d,%d)", Mesh.I, Mesh.J); if (ASetVelToZero) newV=0.0; else newV=copysign(AMaxVelVComp, newV); } if (ARhoF*wF.IR()