diff --git a/integration/genAbs/allMake b/integration/genAbs/allMake new file mode 100755 index 0000000000000000000000000000000000000000..defec71d8cbb1978c8b10f993f102a4546c73702 --- /dev/null +++ b/integration/genAbs/allMake @@ -0,0 +1,33 @@ +#!/bin/bash + +if [ $WM_PROJECT == "foam" ]; then + ofversion=`echo $WM_PROJECT_VERSION | sed -e 's/\.x/-9/' -e 's/\./\'$'\n/g' -e 's/-/\'$'\n/g' | grep "[0-9]" | head -2 | tr -d '\n'` +else + ofversion=`echo $WM_PROJECT_VERSION"-0" | sed -e 's/\.x/-9/' -e 's/\./\'$'\n/g' -e 's/-/\'$'\n/g' -e 's/+/\'$'\n/g' -e 's/v/\'$'\n/g' | grep "[0-9]" | head -3 | tr -d '\n'` +fi + +#IHC_dbg +echo $ofversion +#---- + +export OF_VERSION=$ofversion + +wclean all > /dev/null + +wmake libso waveGeneration + +if (( $? )) ; then + echo "\n\nWave generation boundary conditions (V2.0-ESI) compilation failed" + exit; else + echo -e "\n\nIH Wave generationV2.0 boundary conditions (V2.0-ESI) compiled successfully for $WM_PROJECT $WM_PROJECT_VERSION\n\n\n"; +fi + +wmake libso waveAbsorption + +if (( $? )) ; then + echo "\n\nWave absorption boundary conditions (V2.0-ESI) compilation failed" + exit; else + echo -e "\n\nIH Wave absorption boundary conditions (V2.0-ESI) compiled successfully for $WM_PROJECT $WM_PROJECT_VERSION\n\n\n"; +fi + +wclean all > /dev/null diff --git a/integration/genAbs/common/calculateWaterLevel/calculatedLevelRegular.H b/integration/genAbs/common/calculateWaterLevel/calculatedLevelRegular.H new file mode 100644 index 0000000000000000000000000000000000000000..3c4f92b7b9a6ceed3afa40393e120053437e1635 --- /dev/null +++ b/integration/genAbs/common/calculateWaterLevel/calculatedLevelRegular.H @@ -0,0 +1,88 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + + if ( waveTheory_ == "StokesI" ) + { + forAll(calculatedLevel, itS1) + { + calculatedLevel[itS1] = RealwaterDepth_ + + timeMult * StokesIFun :: eta + ( + waveHeight_, + waveKx, + xGroup[itS1], + waveKy, + yGroup[itS1], + waveOmega, + currTime, + wavePhase_ + ); + } + } + else if ( waveTheory_ == "StokesII" ) + { + forAll(calculatedLevel, itS2) + { + calculatedLevel[itS2] = RealwaterDepth_ + + timeMult * StokesIIFun :: eta + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[itS2], + waveKy, + yGroup[itS2], + waveOmega, + currTime, + wavePhase_ + ); + } + } + else if ( waveTheory_ == "StokesV" ) + { + forAll(calculatedLevel, it1) + { + calculatedLevel[it1] = RealwaterDepth_ + + timeMult * stokesVFun :: eta + ( + RealwaterDepth_, + waveKx, + waveKy, + lambdaStokesV_, + wavePeriod_, + xGroup[it1], + yGroup[it1], + currTime, + wavePhase_ + ); + } + } + else if ( waveTheory_ == "Cnoidal" ) + { + forAll(calculatedLevel, it3) + { + calculatedLevel[it3] = RealwaterDepth_ + + timeMult * cnoidalFun :: eta + ( + waveHeight_, + mCnoidal_, + waveKx, + waveKy, + wavePeriod_, + xGroup[it3], + yGroup[it3], + currTime + ); + } + } + else + { + FatalError << "Wave theory not supported, use:\n" + << "StokesI, StokesII, StokesV, Cnoidal, SolitaryBoussinesq" + << exit(FatalError); + } diff --git a/integration/genAbs/common/calculateWaterLevel/calculatedLevelRegularNormal.H b/integration/genAbs/common/calculateWaterLevel/calculatedLevelRegularNormal.H new file mode 100644 index 0000000000000000000000000000000000000000..09fa19c5edfaf1d3adf7326b8129a34a86e131c7 --- /dev/null +++ b/integration/genAbs/common/calculateWaterLevel/calculatedLevelRegularNormal.H @@ -0,0 +1,76 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + + if ( waveTheory_ == "StokesI" ) + { + calculatedLevel = RealwaterDepth_ + + timeMult * StokesIFun :: eta + ( + waveHeight_, + waveKx, + xGroup[1], + waveKy, + yGroup[1], + waveOmega, + currTime, + wavePhase_ + ); + } + else if ( waveTheory_ == "StokesII" ) + { + calculatedLevel = RealwaterDepth_ + + timeMult * StokesIIFun :: eta + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[1], + waveKy, + yGroup[1], + waveOmega, + currTime, + wavePhase_ + ); + } + else if ( waveTheory_ == "StokesV" ) + { + calculatedLevel = RealwaterDepth_ + + timeMult * stokesVFun :: eta + ( + RealwaterDepth_, + waveKx, + waveKy, + lambdaStokesV_, + wavePeriod_, + xGroup[1], + yGroup[1], + currTime, + wavePhase_ + ); + } + else if ( waveTheory_ == "Cnoidal" ) + { + calculatedLevel = RealwaterDepth_ + + timeMult * cnoidalFun :: eta + ( + waveHeight_, + mCnoidal_, + waveKx, + waveKy, + wavePeriod_, + xGroup[1], + yGroup[1], + currTime + ); + } + else + { + FatalError << "Wave theory not supported, use:\n" + << "StokesI, StokesII, StokesV, Cnoidal, SolitaryBoussinesq." + << exit(FatalError); + } diff --git a/integration/genAbs/common/calculateWaterLevel/calculatedLevelSolitary.H b/integration/genAbs/common/calculateWaterLevel/calculatedLevelSolitary.H new file mode 100644 index 0000000000000000000000000000000000000000..10ee599eb046953c5357a032ca154d32855efc2f --- /dev/null +++ b/integration/genAbs/common/calculateWaterLevel/calculatedLevelSolitary.H @@ -0,0 +1,30 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + + if ( waveTheory_ == "Boussinesq" ) + { + forAll(calculatedLevel, it2) + { + calculatedLevel[it2] = RealwaterDepth_ + + BoussinesqFun :: eta + ( + waveHeight_, + RealwaterDepth_, + xGroup[it2], + yGroup[it2], + waveAngle, + currTime, + X0 + ); + } + } + else + { + FatalError << "Wave theory not supported, use:\n" + << "Boussinesq" << exit(FatalError); + } diff --git a/integration/genAbs/common/checks/checkInputErrorsRegular.H b/integration/genAbs/common/checks/checkInputErrorsRegular.H new file mode 100644 index 0000000000000000000000000000000000000000..293ad83116ce6d5428ba8916807a777b620002fd --- /dev/null +++ b/integration/genAbs/common/checks/checkInputErrorsRegular.H @@ -0,0 +1,159 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + +if (nPaddles_ < 1) +{ + FatalError << "Check nPaddles value." << exit(FatalError); +} + +if ( nPaddles_ > 1 ) +{ + nPaddles_ = decreaseNPaddles( nPaddles_, patchD, dMin, dSpan ); + reduce(nPaddles_, minOp<label>()); +} + +if ( waveHeight_ < 0.0 ) +{ + FatalError << "Check wave height value." << exit(FatalError); +} + +if ( waveTheory_ == "aaa" ) +{ + FatalError << "Wave theory not specified." << exit(FatalError); +} +else if ( waveTheory_ == "StokesI" || waveTheory_ == "StokesII" ) +{ + if ( wavePeriod_ <= 0.0 ) + { + FatalError << "Check wave period value." << exit(FatalError); + } + + if ( waveLength_ <= 0.0 ) + { + word generation = "a"; + + waveLength_ = StokesIFun::waveLength (RealwaterDepth_, wavePeriod_); + + waveK = 2.0*PI()/waveLength_; + + if ( waveK*RealwaterDepth_ > PI() ) + { + generation = "Deep"; + } + else if ( waveK*RealwaterDepth_ < PI()/10.0 ) + { + generation = "Shallow"; + } + else + { + generation = "Intermediate"; + } + + Info << "\nIH Wave Generation BC" << endl; + Info << " Wave theory: " << waveTheory_ << endl; + Info << " H: " << waveHeight_ << endl; + Info << " T: " << wavePeriod_ << endl; + Info << " h: " << RealwaterDepth_ << endl; + Info << " waveLength: " << waveLength_ << endl; + Info << " Direction: " << waveDir_ << "º" << endl; + Info << " Generation in: " << generation << " waters." << endl; + Info << " Relative depth (kh): " << waveK*RealwaterDepth_ + << "\n\n" << endl; + } +} +else if ( waveTheory_ == "StokesV" ) +{ + if ( wavePeriod_ <= 0.0 ) + { + FatalError << "Check wave period value." << exit(FatalError); + } + + if ( lambdaStokesV_ <= 0.0 ) + { + word generation = "a"; + + scalar f1; + scalar f2; + + stokesVFun :: StokesVNR ( waveHeight_, + RealwaterDepth_, + wavePeriod_, + &waveK, + &lambdaStokesV_, + &f1, + &f2 + ); + + Info << "f1 residual " << f1 << endl; + Info << "f2 residual " << f2 << endl; + + if ( f1 > 0.001 || f2 > 0.001 ) + { + FatalError << "No convergence for Stokes V wave theory.\n" + << exit(FatalError); + } + + waveLength_ = 2.0*PI()/waveK; + + if ( waveK*RealwaterDepth_ > PI() ) + { + generation = "Deep"; + } + else if ( waveK*RealwaterDepth_ < PI()/10.0 ) + { + generation = "Shallow"; + } + else + { + generation = "Intermediate"; + } + + Info << "\nIH Wave Generation BC" << endl; + Info << " Wave theory: " << waveTheory_ << endl; + Info << " H: " << waveHeight_ << endl; + Info << " T: " << wavePeriod_ << endl; + Info << " h: " << RealwaterDepth_ << endl; + Info << " L: " << waveLength_ << endl; + Info << " Lambda: " << lambdaStokesV_ << endl; + Info << " Direction: " << waveDir_ << "º" << endl; + Info << " Generation in: " << generation << " waters." << endl; + Info << " Relative depth (kh): " << waveK*RealwaterDepth_ + << "\n\n" << endl; + } +} +else if ( waveTheory_ == "Cnoidal" ) +{ + if ( wavePeriod_ <= 0.0 ) + { + FatalError << "Check wave period value." << exit(FatalError); + } + + if ( mCnoidal_ <= 0.0 ) + { + cnoidalFun :: calculations ( waveHeight_, + RealwaterDepth_, + wavePeriod_, + &mCnoidal_, + &waveLength_); + + Info << "\nIH Wave Generation BC" << endl; + Info << " Wave theory: " << waveTheory_ << endl; + Info << " H: " << waveHeight_ << endl; + Info << " T: " << wavePeriod_ << endl; + Info << " h: " << RealwaterDepth_ << endl; + Info << " L: " << waveLength_ << endl; + Info << " m parameter: " << mCnoidal_ << endl; + Info << " Direction: " << waveDir_ << "º" << "\n\n" << endl; + } +} +else +{ + FatalError << "Wave theory not supported, use:\n" + << "StokesI, StokesII, StokesV, Cnoidal." << exit(FatalError); +} + diff --git a/integration/genAbs/common/checks/checkInputErrorsSolitary.H b/integration/genAbs/common/checks/checkInputErrorsSolitary.H new file mode 100644 index 0000000000000000000000000000000000000000..62061b6beb23104c154a15805ca0c119066eb02a --- /dev/null +++ b/integration/genAbs/common/checks/checkInputErrorsSolitary.H @@ -0,0 +1,41 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + + if (nPaddles_ < 1) + { + FatalError << "Check nPaddles value." << exit(FatalError); + } + + if ( nPaddles_ > 1 ) + { + nPaddles_ = decreaseNPaddles( nPaddles_, patchD, dMin, dSpan ); + reduce(nPaddles_, minOp<label>()); + } + + if ( waveHeight_ < 0.0 ) + { + FatalError << "Check wave height value." << exit(FatalError); + } + + if ( waveTheory_ == "aaa" ) + { + FatalError << "Wave theory not specified." << exit(FatalError); + } + else if ( waveTheory_ == "Boussinesq" ) + { + Info << "\nIH Wave Generation BC" << endl; + Info << "Wave theory: " << waveTheory_ << endl; + Info << "H: " << waveHeight_ << endl; + Info << "h: " << RealwaterDepth_ << endl; + Info << "Direction: " << waveDir_ << "º" << "\n\n" << endl; + } + else + { + FatalError << "Wave theory not supported, use:\n" + << "Boussinesq" << exit(FatalError); + } diff --git a/integration/genAbs/common/memberFun.H b/integration/genAbs/common/memberFun.H new file mode 100644 index 0000000000000000000000000000000000000000..34486ba69231a567e25986658f6a8f242adcb3bc --- /dev/null +++ b/integration/genAbs/common/memberFun.H @@ -0,0 +1,355 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + + //- Return Pi + scalar PI() + { + #if OFVERSION >= 200 + const scalar PI = constant::mathematical::pi; + #else + const scalar PI = mathematicalConstant::pi; + #endif + return PI; + } + + //- Return alphaName + word alphaName() + { + #if OFVERSION >= 230 + const word an = "alpha.water"; + #else + const word an = "alpha1"; + #endif + return an; + } + + //- Return prevalent direction - X or Y to sort groups of columns + scalarField patchDirection ( vector span, scalar* dMin, scalar* dSpan ) + { + const pointField& points = patch().patch().localPoints(); + if (span[0] >= span[1]) + { + *dMin = gMin(points.component(0)); + *dSpan = gMax(points.component(0)) - *dMin; + return patch().Cf().component(0); + } + else + { + *dMin = gMin(points.component(1)); + *dSpan = gMax(points.component(1)) - *dMin; + return patch().Cf().component(1); + } + } + + //- Return Zmax and Zmin of the patch faces + void faceBoundsZ(scalarField* zSup, scalarField* zInf) + { + const label nF = patch().faceCells().size(); + scalarField zMax = Foam::scalarField(nF, -9999.0); + scalarField zMin = Foam::scalarField(nF, 9999.0); + + const faceList& faces = this->patch().patch().localFaces(); + const List<vector>& points = this->patch().patch().localPoints(); + + forAll( faces, faceI ) + { + const face& f = faces[faceI]; + forAll(f, point) + { + scalar auxiliar = points[f[point]].component(2); + + zMax[faceI] = max(zMax[faceI], auxiliar); + zMin[faceI] = min(zMin[faceI], auxiliar); + } + } + + *zSup = zMax; + *zInf = zMin; + } + + //- Calculate water levels on each paddle (same zSpan) + scalarList calcWL + ( + scalarField alphaCell, + labelList cellGroup, + scalar zSpan + ) + { + scalarList groupTotalArea (nPaddles_,0.0); + scalarList groupWaterArea (nPaddles_,0.0); + scalarList heights (nPaddles_,0.0); + + const scalarField faceSurface = patch().magSf(); + + forAll( patch().faceCells(), index ) + { + groupTotalArea[cellGroup[index]-1] + += faceSurface[index]; + groupWaterArea[cellGroup[index]-1] + += faceSurface[index]*alphaCell[index]; + } + + for (label i=0; i<=gMax(cellGroup)-1; i++) + { + reduce(groupTotalArea[i], sumOp<scalar>()); + reduce(groupWaterArea[i], sumOp<scalar>()); + heights[i] = groupWaterArea[i]/groupTotalArea[i]*zSpan; + } + + return heights; + } + + //- Calculate mean velocities on each paddle + void calcUV + ( + scalarField alphaCell, + labelList cellGroup, + scalarField UxCell, + scalarList* Ux, + scalarField UyCell, + scalarList* Uy + ) + { + scalarList groupWaterArea (nPaddles_,0.0); + scalarList groupUx (nPaddles_,0.0); + scalarList groupUy (nPaddles_,0.0); + + const scalarField faceSurface = patch().magSf(); + + forAll( patch().faceCells(), index ) + { + groupWaterArea[cellGroup[index]-1] += + faceSurface[index]*alphaCell[index]; + + groupUx[cellGroup[index]-1] += + UxCell[index]*alphaCell[index]*faceSurface[index]; + groupUy[cellGroup[index]-1] += + UyCell[index]*alphaCell[index]*faceSurface[index]; + } + + for (label i=0; i<=nPaddles_-1; i++) + { + reduce(groupWaterArea[i], sumOp<scalar>()); + reduce(groupUx[i], sumOp<scalar>()); + reduce(groupUy[i], sumOp<scalar>()); + + groupUx[i] = groupUx[i]/groupWaterArea[i]; + groupUy[i] = groupUy[i]/groupWaterArea[i]; + } + + *Ux = groupUx; + *Uy = groupUy; + } + + //- Mean of a scalarList + scalar meanSL ( scalarList lst ) + { + scalar aux = 0.0; + + forAll(lst,i) + { + aux += lst[i]; + } + + return aux/scalar(lst.size()); + } + + //- Reduce number of paddles if too high (paddles without cells) + label decreaseNPaddles + ( + label np, + scalarField patchD, + scalar dMin, + scalar dSpan + ) + { + scalarList dBreakPoints(np+1, 0.0); + + while(1) + { + scalarList usedGroup (np,0.0); + + for (label i=0; i<=np; i++) + { + dBreakPoints[i] = dMin + dSpan/(np)*i; + } + + forAll(patch().faceCells(), patchCells) + { + for (label j=0; j<=np; j++) + { + if ( (patchD[patchCells]>=dBreakPoints[j]) && + (patchD[patchCells]<dBreakPoints[j+1]) ) + { + usedGroup[j] = 1.0; + continue; + } + } + } + + reduce(usedGroup, maxOp<scalarList>()); + + if ( gMin(usedGroup) < 0.5 ) + { + Info + << "Reduced nPaddles from " << np + << " to " << np-1 << endl; + np -= 1; + } + else + { + break; + } + } + return np; + } + + //- Own convenient definition of atan + scalar arcTan ( scalar x, scalar y ) + { + scalar A = 0.0; + + if( x > 0.0 ) + { + A = atan(y/x); + } + else if ( y >= 0.0 && x < 0.0 ) + { + A = PI() + atan(y/x); + } + else if ( y < 0.0 && x < 0.0 ) + { + A = -PI() + atan(y/x); + } + else if ( y > 0.0 && x == 0.0 ) + { + A = PI()/2.0; + } + else if ( y < 0.0 && x == 0.0 ) + { + A = -PI()/2.0; + } + else if ( y == 0.0 && x == 0.0 ) + { + A = 0; + } + + return A; + } + + //- Calculate mean horizontal angle for each paddle + scalarList meanPatchDirs ( labelList cellGroup ) + { + const vectorField nVecCell = patch().nf(); + + scalarList numAngles (nPaddles_, 0.0); + scalarList meanAngle (nPaddles_, 0.0); + + forAll(patch().faceCells(), patchCells) + { + meanAngle[cellGroup[patchCells]-1] + += arcTan( nVecCell[patchCells].component(0), + nVecCell[patchCells].component(1)); + numAngles[cellGroup[patchCells]-1] += 1.0; + } + + for (label i=0; i<=nPaddles_-1; i++) + { + reduce(meanAngle[i], sumOp<scalar>()); + reduce(numAngles[i], sumOp<scalar>()); + meanAngle[i] = meanAngle[i]/numAngles[i] + PI(); + } + return meanAngle; + } + + //- Calculate z-bounds for each paddle + scalarList zSpanList + ( + labelList cellGroup, + scalarField zInf, + scalarField zSup + ) + { + scalarList zMaxL (nPaddles_, -9999.0); + scalarList zMinL (nPaddles_, 9999.0); + + forAll(patch().faceCells(), patchCells) + { + zMaxL[cellGroup[patchCells]-1] = + max(zMaxL[cellGroup[patchCells]-1], zSup[patchCells] ); + zMinL[cellGroup[patchCells]-1] = + min(zMinL[cellGroup[patchCells]-1], zInf[patchCells] ); + } + + for (label i=0; i<=nPaddles_-1; i++) + { + reduce(zMaxL[i], maxOp<scalar>()); + reduce(zMinL[i], minOp<scalar>()); + } + return zMaxL-zMinL; + } + + //- In-line print for scalarLists + label inlinePrint ( std::string name, scalarList SL ) + { + Info << name << " " << SL.size() << "( "; + forAll(SL, i) + { + Info << SL[i] << " "; + } + Info << ")\n" << endl; + return 0; + } + + //- Simple linear interpolation + scalar interpolation + ( + scalar x1, + scalar x2, + scalar y1, + scalar y2, + scalar xInt + ) + { + scalar yInt = y1 + (y2-y1)/(x2-x1)*(xInt-x1); + return yInt; + } + + //- Limit angle between -pi/2 and pi/2 + scalar limAngle (scalar ang) + { + ang = abs(ang); + + while (ang >= 2.0*PI()) + { + ang -= 2.0*PI(); + } + + if ( ang >= PI()/2.0 && ang <= 3.0*PI()/2.0 ) + { + return PI()/2.0; + } + else + { + return ang; + } + } + + //- Interpolation BO style + scalar interpolatte + ( + scalar x1, + scalar x2, + scalar y1, + scalar y2, + scalar xInt + ) + { + scalar yInt = (y1 + y2)*0.5; + return yInt; + } diff --git a/integration/genAbs/common/waveFun.C b/integration/genAbs/common/waveFun.C new file mode 100644 index 0000000000000000000000000000000000000000..23be07e0a7b7e1d196822b8d055b0edb930f9858 --- /dev/null +++ b/integration/genAbs/common/waveFun.C @@ -0,0 +1,980 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + +#include "waveFun.H" +#include <math.h> +#include <stdlib.h> +#include <stdio.h> + +namespace otherFun +{ + #define PII 3.1415926535897932384626433832795028 + #define grav 9.81 + + double interpolation (double x1, double x2, double y1, double y2, double xInt) + { + + double yInt = y1 + (y2-y1)/(x2-x1)*(xInt-x1); + + return yInt; + } +} + +namespace StokesIFun +{ + #define PII 3.1415926535897932384626433832795028 + #define G 9.81 + + double waveLength (double h, double T) + { + double L0 = G*T*T/(2.0*PII); + double L = L0; + + for(int i=1; i<=100; i++) + { + L = L0*tanh(2.0*PII*h/L); + } + + return L; + } + + double eta (double H, double Kx, double x, double Ky, double y, double omega, double t, double phase) + { + double faseTot = Kx*x + Ky*y - omega*t + phase; + + double sup = H*0.5*cos(faseTot); + + return sup; + } + + double U (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z) + { + double k = sqrt(Kx*Kx + Ky*Ky); + double faseTot = Kx*x + Ky*y - omega*t + phase; + + double velocity = H*0.5*omega*cos(faseTot)*cosh(k*z)/sinh(k*h); + + return velocity; + } + + double W (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z) + { + double k = sqrt(Kx*Kx + Ky*Ky); + double faseTot = Kx*x + Ky*y - omega*t + phase; + + double velocity = H*0.5*omega*sin(faseTot)*sinh(k*z)/sinh(k*h); + + return velocity; + } +} + +namespace StokesIIFun +{ + #define PII 3.1415926535897932384626433832795028 + #define G 9.81 + + double eta (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase) + { + double k = sqrt(Kx*Kx + Ky*Ky); + double sigma = tanh(k*h); + double faseTot = Kx*x + Ky*y - omega*t + phase; + + double sup = H*0.5*cos(faseTot) + k*H*H/4.0*(3.0-sigma*sigma)/(4.0*pow(sigma,3))*cos(2.0*faseTot); + + return sup; + } + + double U (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z) + { + double k = sqrt(Kx*Kx + Ky*Ky); + double faseTot = Kx*x + Ky*y - omega*t + phase; + + double velocity = H*0.5*omega*cos(faseTot)*cosh(k*z)/sinh(k*h) + 3.0/4.0*H*H/4.0*omega*k*cosh(2.0*k*z)/pow(sinh(k*h),4)*cos(2.0*faseTot); + + return velocity; + } + + double W (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z) + { + double k = sqrt(Kx*Kx + Ky*Ky); + double faseTot = Kx*x + Ky*y - omega*t + phase; + + double velocity = H*0.5*omega*sin(faseTot)*sinh(k*z)/sinh(k*h) + 3.0/4.0*H*H/4.0*omega*k*sinh(2.0*k*z)/pow(sinh(k*h),4)*sin(2.0*faseTot); + + return velocity; + } +} + +namespace Elliptic +{ + #define PII 3.1415926535897932384626433832795028 + #define ITER 25 + + void ellipticIntegralsKE (double m, double* K, double* E) + { + long double a,aOld,g,gOld,aux,sum; + + if ( m == 0.0 ) + { + *K = PII/2.; + *E = PII/2.; + return; + } + + a = 1.0L; + g = sqrt(1.0L - m); + aux = 1.0L; + sum = 2.0L - m; + + while (1) + { + gOld = g; + aOld = a; + a = 0.5L * (gOld + aOld); + g = gOld * aOld; + aux += aux; + sum -= aux * (a * a - g); + + if ( fabs(aOld - gOld) <= (aOld * 1.e-22) ) + { + break; + } + + g = sqrt(g); + } + + *K = (double) (PII/2. / a); + *E = (double) (PII/4. / a * sum); + return; + } + + double JacobiAmp (double u, double m) + { + long double a[ITER+1], g[ITER+1], c[ITER+1]; + long double aux, amp; + int n; + + m = fabsl(m); + + if ( m == 0.0 ) + { + return u; + } + + if ( m == 1. ) + { + return 2. * atan( exp(u) ) - PII/2.; + } + + a[0] = 1.0L; + g[0] = sqrtl(1.0L - m); + c[0] = sqrtl(m); + + aux = 1.0L; + + for (n = 0; n < ITER; n++) + { + if ( fabsl(a[n] - g[n]) < (a[n] * 1.e-22) ) + { + break; + } + + aux += aux; + a[n+1] = 0.5L * (a[n] + g[n]); + g[n+1] = sqrtl(a[n] * g[n]); + c[n+1] = 0.5L * (a[n] - g[n]); + } + + amp = aux * a[n] * u; + + for (; n > 0; n--) + { + amp = 0.5L * ( amp + asinl( c[n] * sinl(amp) / a[n]) ); + } + + return (double) amp; + } + + void JacobiSnCnDn (double u, double m, double* Sn, double* Cn, double* Dn) + { + double amp = Elliptic::JacobiAmp( u, m ); + + *Sn = sin( amp ); + *Cn = cos( amp ); + *Dn = sqrt(1.0 - m * sin( amp ) * sin( amp )); + return; + } + +} + +namespace cnoidalFun +{ + #define PII 3.1415926535897932384626433832795028 + #define G 9.81 + + double eta (double H, double m, double kx, double ky, double T, double x, double y, double t) + { + double K, E; + Elliptic::ellipticIntegralsKE(m, &K, &E); + + double uCnoidal = K/PII*(kx*x + ky*y - 2.0*PII*t/T); + + double sn, cn, dn; + Elliptic::JacobiSnCnDn(uCnoidal, m, &sn, &cn, &dn); + + double etaCnoidal = H*((1.0-E/K)/m - 1.0 + pow(cn,2)); + + return etaCnoidal; + } + + double etaCnoidal1D (double H, double m, double t, double T) + { + double K, E; + Elliptic::ellipticIntegralsKE(m, &K, &E); + + double uCnoidal = -2.0*K*(t/T); + + double sn, cn, dn; + Elliptic::JacobiSnCnDn(uCnoidal, m, &sn, &cn, &dn); + + double etaCnoidal = H*((1.0-E/K)/m - 1.0 + pow(cn,2)); + + return etaCnoidal; + } + + double etaMeanSq (double H, double m, double T) + { + double eta = 0; + double etaSumSq = 0; + + for(int i=0; i<1000; i++) + { + eta = etaCnoidal1D(H, m, i*T/(1000.0), T); + etaSumSq += eta*eta; + } + + etaSumSq /= 1000.0; + return etaSumSq; + } + + double d1EtaDx (double H, double m, double uCnoidal, double L, double K, double E) + { + double dudx = 0; + dudx = 2.0*K/L; + + double sn, cn, dn; + Elliptic::JacobiSnCnDn(uCnoidal, m, &sn, &cn, &dn); + + double deriv = -2.0*H*cn*dn*sn*dudx; + + return deriv; + } + + double d2EtaDx (double H, double m, double uCnoidal, double L, double K, double E) + { + double dudxx = 0; + dudxx = 4.0*K*K/L/L; + + double sn, cn, dn; + Elliptic::JacobiSnCnDn(uCnoidal, m, &sn, &cn, &dn); + + double deriv = 2.0*H*(dn*dn*sn*sn - cn*cn*dn*dn + m*cn*cn*sn*sn)*dudxx; + + return deriv; + } + + double d3EtaDx (double H, double m, double uCnoidal, double L, double K, double E) + { + double dudxxx = 0; + dudxxx = 8.0*K*K*K/L/L/L; + + double sn, cn, dn; + Elliptic::JacobiSnCnDn(uCnoidal, m, &sn, &cn, &dn); + + double deriv = 8.0*H*( cn*sn*dn*dn*dn*(-4.0 -2.0*m) + 4.0*m*cn*sn*sn*sn*dn -2.0*m*cn*cn*cn*sn*dn )*dudxxx; + + return deriv; + } + + int calculations (double H, double d, double T, double* mOut, double* LOut) + { + double mTolerance = 0.0001; + double mElliptic = 0.5; + double LElliptic = 0; + double phaseSpeed = 0; + + double mError = 0.0; + double mMinError = 999; + + while (mElliptic < 1.0) + { + double KElliptic, EElliptic; + Elliptic::ellipticIntegralsKE(mElliptic, &KElliptic, &EElliptic); + + LElliptic = KElliptic*sqrt(16.0*pow(d,3)*mElliptic/(3.0*H)); + + phaseSpeed = sqrt(G*d)*(1.0 - H/d/2.0 + H/d/mElliptic*(1.0-3.0/2.0*EElliptic/KElliptic)); + + mError = fabs(T-LElliptic/phaseSpeed); + + if (mError <= mMinError) + { + *mOut = mElliptic; + *LOut = LElliptic; + mMinError = mError; + } + + mElliptic += mTolerance; + } + + return 0; + } + + double U (double H, double h, double m, double kx, double ky, double T, double x, double y, double t, double z) + { + double K, E; + Elliptic::ellipticIntegralsKE(m, &K, &E); + + double uCnoidal = K/PII*(kx*x + ky*y - 2.0*PII*t/T); + double k = sqrt(kx*kx + ky*ky); + double L = 2.0*PII/k; + double c = L/T; + + double etaCN = eta(H, m, kx, ky, T, x, y, t); + double etaXX = d2EtaDx(H, m, uCnoidal, L, K, E); + double etaMS = etaMeanSq(H, m, T); + + double velocity = c*etaCN/h - c*(etaCN*etaCN/h/h + etaMS*etaMS/h/h) + 1.0/2.0*c*h*(1.0/3.0 - z*z/h/h)*etaXX; + + return velocity; + } + + double W (double H, double h, double m, double kx, double ky, double T, double x, double y, double t, double z) + { + double K, E; + Elliptic::ellipticIntegralsKE(m, &K, &E); + + double uCnoidal = K/PII*(kx*x + ky*y - 2.0*PII*t/T); + double k = sqrt(kx*kx + ky*ky); + double L = 2.0*PII/k; + double c = L/T; + + double etaCN = eta(H, m, kx, ky, T, x, y, t); + double etaX = d1EtaDx(H, m, uCnoidal, L, K, E); + double etaXXX = d3EtaDx(H, m, uCnoidal, L, K, E); + + + double velocity = -c*z*( etaX/h*(1.0-2.0*etaCN/h) + 1.0/6.0*h*(1.0-z*z/h/h)*etaXXX ) ; + + return velocity; + } +} + +namespace stokesVFun +{ + #define PII 3.1415926535897932384626433832795028 + #define G 9.81 + + double A11 (double h, double k) + { + double s = sinh(k*h); + + double A = 1.0/s; + + return A; + } + + double A13 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double A = + -pow(c,2)*(5.0*pow(c,2)+1.0) + / + (8.0*pow(s,5)); + + return A; + } + + double A15 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double A = + -(1184.0*pow(c,10)-1440.0*pow(c,8)-1992.0*pow(c,6)+2641.0*pow(c,4)-249.0*pow(c,2)+18) + / + (1536.0*pow(s,11)); + + return A; + } + + double A22 (double h, double k) + { + double s = sinh(k*h); + + double A = + 3.0 + / + (8.0*pow(s,4)); + + return A; + } + + double A24 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double A = + (192.0*pow(c,8)-424.0*pow(c,6)-312.0*pow(c,4)+480.0*pow(c,2)-17) + / + (768.0*pow(s,10)); + + return A; + } + + double A33 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double A = + (13.0-4.0*pow(c,2)) + / + (64.0*pow(s,7)); + + return A; + } + + double A35 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double A = + (512.0*pow(c,12)+4224.0*pow(c,10)-6800.0*pow(c,8)-12808.0*pow(c,6)+16704.0*pow(c,4)-3154.0*pow(c,2)+107.0) + / + (4096.0*pow(s,13)*(6.0*pow(c,2)-1.0)); + + return A; + } + + double A44 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double A = + (80.0*pow(c,6)-816.0*pow(c,4)+1338.0*pow(c,2)-197.0) + / + (1536.0*pow(s,10)*(6.0*pow(c,2)-1.0)); + + return A; + } + + double A55 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double A = + -(2880.0*pow(c,10)-72480.0*pow(c,8)+324000.0*pow(c,6)-432000.0*pow(c,4)+163470.0*pow(c,2)-16245.0) + / + (61440.0*pow(s,11)*(6.0*pow(c,2)-1.0)*(8.0*pow(c,4)-11.0*pow(c,2)+3.0)); + + return A; + } + + double B22 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double B = + (2.0*pow(c,2)+1)*c + / + (4.0*pow(s,3)); + + return B; + } + + double B24 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double B = + (272.0*pow(c,8)-504.0*pow(c,6)-192.0*pow(c,4)+322.0*pow(c,2)+21.0)*c + / + (384.0*pow(s,9)); + + return B; + } + + double B33 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double B = + (8.0*pow(c,6)+1.0)*3.0 + / + (64.0*pow(s,6)); + + return B; + } + + double B33k (double h, double k) // d B33 / d k + { + double s = sinh(k*h); + double c = cosh(k*h); + + double sk = h*s; + double ck = h*c; + + double B = + 9.0*pow(c,5)*ck/(4.0*pow(s,6))- + (9.0*(8.0*pow(c,6)+1.0))/(32.0*pow(s,7))*sk; + + return B; + } + + double B35 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double B = + (88128.0*pow(c,14)-208224.0*pow(c,12)+70848.0*pow(c,10)+54000.0*pow(c,8)-21816.0*pow(c,6)+6264.0*pow(c,4)-54.0*pow(c,2)-81) + / + (12288.0*pow(s,12)*(6.0*pow(c,2)-1.0)); + + return B; + } + + double B35k (double h, double k) // d B35 / d k + { + double s = sinh(k*h); + double c = cosh(k*h); + + double sk = h*s; + double ck = h*c; + + double B = + (14.0*88128.0*pow(c,13)*ck-12.0*208224.0*pow(c,11)*ck + +10.0*70848.0*pow(c,9)*ck+8.0*54000.0*pow(c,7)*ck + -6.0*21816.0*pow(c,5)*ck+4.0*6264.0*pow(c,3)*ck + -2.0*54.0*c*ck) + /(12288.0*pow(s,12)*(6.0*pow(c,2)-1.0)) + -(88128.0*pow(c,14)-208224.0*pow(c,12)+70848.0*pow(c,10)+54000.0*pow(c,8) + -21816.0*pow(c,6)+6264.0*pow(c,4)-54.0*pow(c,2)-81.0)*12.0 + /(12288.0*pow(s,13)*(6.0*pow(c,2)-1.0))*sk + -(88128.0*pow(c,14)-208224.0*pow(c,12)+70848.0*pow(c,10)+54000.0*pow(c,8) + -21816.0*pow(c,6)+6264.0*pow(c,4)-54.0*pow(c,2)-81.0)*12.0*c*ck + /(12288.0*pow(s,12)*pow((6.0*pow(c,2)-1.0),2)); + + return B; + } + + double B44 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double B = + (768.0*pow(c,10)-448.0*pow(c,8)-48.0*pow(c,6)+48.0*pow(c,4)+106.0*pow(c,2)-21.0)*c + / + (384.0*pow(s,9)*(6.0*pow(c,2)-1.0)); + + return B; + } + + double B55 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double B = + (192000.0*pow(c,16)-262720.0*pow(c,14)+83680.0*pow(c,12)+20160.0*pow(c,10)-7280.0*pow(c,8)+7160.0*pow(c,6)-1800.0*pow(c,4)-1050.0*pow(c,2)+225.0) + / + (12288.0*pow(s,10)*(6.0*pow(c,2)-1.0)*(8.0*pow(c,4)-11.0*pow(c,2)+3.0)); + + return B; + } + + double B55k (double h, double k) // d B55 / d k + { + double s = sinh(k*h); + double c = cosh(k*h); + + double sk = h*s; + double ck = h*c; + + double B = + (16.0*192000.0*pow(c,15)*ck-14.0*262720.0*pow(c,13)*ck + +12.0*83680.0*pow(c,11)*ck+10.0*20160.0*pow(c,9)*ck + -8.0*7280.0*pow(c,7)*ck+6.0*7160.0*pow(c,5)*ck + -4.0*1800.0*pow(c,3)*ck-2.0*1050.0*pow(c,1)*ck) + /(12288.0*pow(s,10)*(6.0*pow(c,2)-1.0)*(8.0*pow(c,4)-11.0*pow(c,2)+3.0)) + -(192000.0*pow(c,16)-262720.0*pow(c,14)+83680.0*pow(c,12)+20160.0*pow(c,10) + -7280.0*pow(c,8)+7160.0*pow(c,6)-1800.0*pow(c,4)-1050.0*pow(c,2)+225.0)*10.0 + /(12288.0*pow(s,11)*(6.0*pow(c,2)-1.0)*(8.0*pow(c,4)-11.0*pow(c,2)+3.0))*sk + -(192000.0*pow(c,16)-262720.0*pow(c,14)+83680.0*pow(c,12)+20160.0*pow(c,10) + -7280.0*pow(c,8)+7160.0*pow(c,6)-1800.0*pow(c,4)-1050.0*pow(c,2)+225.0) + *12.0*c*ck + /(12288.0*pow(s,10)*pow((6.0*pow(c,2)-1.0),2)*(8.0*pow(c,4)-11.0*pow(c,2)+3.0)) + -(192000.0*pow(c,16)-262720.0*pow(c,14)+83680.0*pow(c,12)+20160.0*pow(c,10) + -7280.0*pow(c,8)+7160.0*pow(c,6)-1800.0*pow(c,4)-1050.0*pow(c,2)+225.0) + *(32.0*pow(c,3)-22.0*c)*ck + /(12288.0*pow(s,10)*(6.0*pow(c,2)-1.0)*pow((8.0*pow(c,4)-11.0*pow(c,2)+3.0),2)); + + return B; + } + + double C1 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double C = + (8.0*pow(c,4)-8.0*pow(c,2)+9.0) + / + (8.0*pow(s,4)); + + return C; + } + + double C1k (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double sk = h*s; + double ck = h*c; + + double C = + (4.0*8.0*pow(c,3)*ck-2.0*8.0*c*ck)/(8.0*pow(s,4)) + -(8.0*pow(c,4)-8.0*pow(c,2)+9.0)*4.0*sk/(8.0*pow(s,5)); + + return C; + } + + double C2 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double C = + (3840.0*pow(c,12)-4096.0*pow(c,10) + 2592.0*pow(c,8)-1008.0*pow(c,6)+5944.0*pow(c,4)-1830.0*pow(c,2)+147.0) // - 2592 + / + (512.0*pow(s,10)*(6.0*pow(c,2)-1.0)); + + return C; + } + + double C2k (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double sk = h*s; + double ck = h*c; + + double C = + (12.0*3840.0*pow(c,11)*ck-10.0*4096.0*pow(c,9)*ck + +8.0*2592.0*pow(c,7)*ck-6.0*1008.0*pow(c,5)*ck+ + 4.0*5944.0*pow(c,3)*ck-2.0*1830.0*c*ck) + /(512.0*pow(s,10)*(6.0*pow(c,2)-1.0)) + -(3840.0*pow(c,12)-4096.0*pow(c,10)+2592.0*pow(c,8)-1008.0*pow(c,6)+ + 5944.0*pow(c,4)-1830.0*pow(c,2)+147.0)*10.0*sk/ + (512.0*pow(s,11)*(6.0*pow(c,2)-1.0)) + -(3840.0*pow(c,12)-4096.0*pow(c,10)+2592.0*pow(c,8)-1008.0*pow(c,6)+ + 5944.0*pow(c,4)-1830.0*pow(c,2)+147.0)*12.0*c*ck + /(512.0*pow(s,10)*pow((6.0*pow(c,2)-1.0),2)); + + return C; + } + + double C3 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double C = + (-1.0) + / + (4.0*s*c); + + return C; + } + + double C4 (double h, double k) + { + double s = sinh(k*h); + double c = cosh(k*h); + + double C = + (12.0*pow(c,8)+36.0*pow(c,6)-162.0*pow(c,4)+141.0*pow(c,2)-27.0) + / + (192.0*c*pow(s,9)); + + return C; + } + + int StokesVNR (double H, double d, double T, double* kOut, double* LambdaOut, double* f1Out, double* f2Out ) + { + double f1 = 1; + double f2 = 1; + + double k = 2.0*PII/(sqrt(G*d)*T); + double Lambda = H/2.0*k; + + double Bmat11, Bmat12, Bmat21, Bmat22; + double b33, b35, b55, b33k, b35k, b55k; + double c1, c2, c1k, c2k; + double kPr, LambdaPr; + + int n = 0; + + while ( (fabs(f1)>1.0e-12 || fabs(f2)>1.0e-12) && n<10000 ) + { + b33 = B33(d, k); + b35 = B35(d, k); + b55 = B55(d, k); + c1 = C1(d, k); + c2 = C2(d, k); + + b33k = B33k(d, k); + b35k = B35k(d, k); + b55k = B55k(d, k); + c1k = C1k(d, k); + c2k = C2k(d, k); + + Bmat11 = 2.0*PII/(pow(k,2)*d)*(Lambda+pow(Lambda,3)*b33+pow(Lambda,5)*(b35+b55)) + -2.0*PII/(k*d)*(pow(Lambda,3)*b33k+pow(Lambda,5)*(b35k+b55k)); + + Bmat12 = -2.0*PII/(k*d)* + (1.0+3.0*pow(Lambda,2)*b33+5.0*pow(Lambda,4)*(b35+b55)); + + Bmat21 = -d/(2.0*PII)*tanh(k*d)*(1.0+pow(Lambda,2)*c1+pow(Lambda,4)*c2) + -k*d/(2.0*PII)*(1.0-pow((tanh(k*d)),2))*d* + (1.0+pow(Lambda,2)*c1+pow(Lambda,4)*c2) + -k*d/(2.0*PII)*tanh(k*d) + *(pow(Lambda,2)*c1k+pow(Lambda,4)*c2k); + + Bmat22 = -k*d/(2.0*PII)*tanh(k*d) + *(2.0*Lambda*c1+4.0*pow(Lambda,3)*c2); + + f1 = PII*H/d-2.0*PII/(k*d)*(Lambda+pow(Lambda,3)*b33+pow(Lambda,5)*(b35+b55)); + + f2 = (2.0*PII*d)/(G*pow(T,2))-k*d/(2.0*PII)*tanh(k*d) + *(1.0+pow(Lambda,2)*c1+pow(Lambda,4)*c2); + + LambdaPr = (f1*Bmat21-f2*Bmat11)/(Bmat11*Bmat22-Bmat12*Bmat21); + kPr = (f2*Bmat12-f1*Bmat22)/(Bmat11*Bmat22-Bmat12*Bmat21); + + Lambda += LambdaPr; + k += kPr; + + n++; + } + + *kOut = k; + *LambdaOut = Lambda; + + *f1Out = fabs(f1); + *f2Out = fabs(f2); + + return 1; + } + + double eta (double h, double kx, double ky, double lambda, double T, double x, double y, double t, double phase) + { + double k = sqrt(kx*kx + ky*ky); + + double b22 = B22(h, k); + double b24 = B24(h, k); + double b33 = B33(h, k); + double b35 = B35(h, k); + double b44 = B44(h, k); + double b55 = B55(h, k); + + double amp1 = lambda/k; + double amp2 = (b22*pow(lambda,2)+b24*pow(lambda,4))/k; + double amp3 = (b33*pow(lambda,3)+b35*pow(lambda,5))/k; + double amp4 = b44*pow(lambda,4)/k; + double amp5 = b55*pow(lambda,5)/k; + + double theta = kx*x + ky*y - 2.0*PII/T*t + phase; + + double C = amp1*cos(theta) + + amp2*cos(2*theta) + + amp3*cos(3*theta) + + amp4*cos(4*theta) + + amp5*cos(5*theta); + + return C; + } + + double U (double d, double kx, double ky, double lambda, double T, double x, double y, double t, double phase, double z) + { + double k = sqrt(kx*kx + ky*ky); + + double a11 = A11(d, k); + double a13 = A13(d, k); + double a15 = A15(d, k); + double a22 = A22(d, k); + double a24 = A24(d, k); + double a33 = A33(d, k); + double a35 = A35(d, k); + double a44 = A44(d, k); + double a55 = A55(d, k); + + double a1u=2.0*PII/T/k*(lambda*a11+pow(lambda,3)*a13+pow(lambda,5)*a15); + double a2u=2.0*2.0*PII/T/k*(pow(lambda,2)*a22+pow(lambda,4)*a24); + double a3u=3.0*2.0*PII/T/k*(pow(lambda,3)*a33+pow(lambda,5)*a35); + double a4u=4.0*2.0*PII/T/k*(pow(lambda,4)*a44); + double a5u=5.0*2.0*PII/T/k*(pow(lambda,5)*a55); + + double velU = 0; + + double theta = kx*x + ky*y - 2.0*PII/T*t + phase; + + velU = a1u*cosh(k*z)*cos(theta) + + a2u*cosh(2.0*k*z)*cos(2.0*(theta)) + + a3u*cosh(3.0*k*z)*cos(3.0*(theta)) + + a4u*cosh(4.0*k*z)*cos(4.0*(theta)) + + a5u*cosh(5.0*k*z)*cos(5.0*(theta)); + + return velU; + } + + double W (double d, double kx, double ky, double lambda, double T, double x, double y, double t, double phase, double z) + { + double k = sqrt(kx*kx + ky*ky); + + double a11 = A11(d, k); + double a13 = A13(d, k); + double a15 = A15(d, k); + double a22 = A22(d, k); + double a24 = A24(d, k); + double a33 = A33(d, k); + double a35 = A35(d, k); + double a44 = A44(d, k); + double a55 = A55(d, k); + + double a1u=2.0*PII/T/k*(lambda*a11+pow(lambda,3)*a13+pow(lambda,5)*a15); + double a2u=2.0*2.0*PII/T/k*(pow(lambda,2)*a22+pow(lambda,4)*a24); + double a3u=3.0*2.0*PII/T/k*(pow(lambda,3)*a33+pow(lambda,5)*a35); + double a4u=4.0*2.0*PII/T/k*(pow(lambda,4)*a44); + double a5u=5.0*2.0*PII/T/k*(pow(lambda,5)*a55); + + double velV = 0; + + double theta = kx*x+ky*y-2.0*PII/T*t+phase; + + velV = a1u*sinh(k*z)*sin(theta) + + a2u*sinh(2.0*k*z)*sin(2.0*(theta)) + + a3u*sinh(3.0*k*z)*sin(3.0*(theta)) + + a4u*sinh(4.0*k*z)*sin(4.0*(theta)) + + a5u*sinh(5.0*k*z)*sin(5.0*(theta)); + + return velV; + } +} + +namespace BoussinesqFun +{ + #define PII 3.1415926535897932384626433832795028 + #define G 9.81 + + double eta (double H, double h, double x, double y, double theta, double t, double X0) + { + double C = sqrt(G*(H+h)); + double ts = 3.5*h/sqrt(H/h); + double aux = sqrt(3.0*H/(4.0*h))/h; + double Xa = -C * t + ts - X0 + x * cos(theta) + y * sin(theta); + + double sup = H * 1.0/pow(cosh( aux * Xa ),2); + + return sup; + } + + double Deta1 (double H, double h, double x, double y, double theta, double t, double X0) + { + double C = sqrt(G*(H+h)); + double ts = 3.5*h/sqrt(H/h); + double aux = sqrt(3.0*H/(4.0*h))/h; + double Xa = -C * t + ts - X0 + x * cos(theta) + y * sin(theta); + + double deta = 8.0*aux*h * exp(2.0*aux*Xa) * (1.0-exp(2.0*aux*Xa)) + / pow(1.0+exp(2.0*aux*Xa),3); + + return deta; + } + + double Deta2 (double H, double h, double x, double y, double theta, double t, double X0) + { + double C = sqrt(G*(H+h)); + double ts = 3.5*h/sqrt(H/h); + double aux = sqrt(3.0*H/(4.0*h))/h; + double Xa = -C * t + ts - X0 + x * cos(theta) + y * sin(theta); + + double deta = 16.0*pow(aux,2)*h * exp(2.0*aux*Xa) * (exp(4.0*aux*Xa) + - 4.0*exp(2.0*aux*Xa)+1.0) / pow(1.0+exp(2.0*aux*Xa),4); + + return deta; + } + + double Deta3 (double H, double h, double x, double y, double theta, double t, double X0) + { + double C = sqrt(G*(H+h)); + double ts = 3.5*h/sqrt(H/h); + double aux = sqrt(3.0*H/(4.0*h))/h; + double Xa = -C * t + ts - X0 + x * cos(theta) + y * sin(theta); + + double deta = -32.0*pow(aux,3)*h * exp(2.0*aux*Xa) * (exp(6.0*aux*Xa) + - 11.0*exp(4.0*aux*Xa) + 11.0*exp(2.0*aux*Xa)-1.0) / pow(1.0+exp(2.0*aux*Xa),5); + + return deta; + } + + double U (double H, double h, double x, double y, double theta, double t, double X0, double z) + { + double C = sqrt(G*(H+h)); + double etaSolit = eta( H, h, x, y, theta, t, X0); + double Detas2 = Deta2( H, h, x, y, theta, t, X0); + + double vel = C*etaSolit/h + *(1.0 - etaSolit/(4.0*h) + + pow(h,2)/(3.0*etaSolit) + *(1.0 - 3.0/2.0*pow(z/h,2)) + *Detas2); + + return vel; + } + + double W (double H, double h, double x, double y, double theta, double t, double X0, double z) + { + double C = sqrt(G*(H+h)); + double etaSolit = eta( H, h, x, y, theta, t, X0); + double Detas1 = Deta1( H, h, x, y, theta, t, X0); + double Detas3 = Deta3( H, h, x, y, theta, t, X0); + + double vel = -C*z/h + *((1.0 - etaSolit/(2.0*h))*Detas1 + + pow(h,2)/3.0 + *(1.0 - 1.0/2.0*pow(z/h,2)) + *Detas3); + + return vel; + } +} + diff --git a/integration/genAbs/common/waveFun.H b/integration/genAbs/common/waveFun.H new file mode 100644 index 0000000000000000000000000000000000000000..53406b64c53f52b557a41cab9a8c826379a07754 --- /dev/null +++ b/integration/genAbs/common/waveFun.H @@ -0,0 +1,95 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + +#ifndef waveFun_H +#define waveFun_H + +namespace otherFun +{ + double interpolation (double x1, double x2, double y1, double y2, double xInt); +} + +namespace StokesIFun +{ + double waveLength (double h, double T); + double eta (double H, double Kx, double x, double Ky, double y, double omega, double t, double phase); + double U (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z); + double W (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z); +} + +namespace StokesIIFun +{ + double eta (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase); + double U (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z); + double W (double H, double h, double Kx, double x, double Ky, double y, double omega, double t, double phase, double z); + double timeLag (double H, double h, double Kx, double x, double Ky, double y, double T, double phase); +} + +namespace Elliptic +{ + void ellipticIntegralsKE (double m, double* K, double* E); + void JacobiSnCnDn (double u, double m, double* Sn, double* Cn, double* Dn); +} + +namespace cnoidalFun +{ + double eta (double H, double m, double kx, double ky, double T, double x, double y, double t); + double etaCnoidal1D (double H, double m, double t, double T, double K, double E); + double etaMeanSq (double H, double m, double T); + double d1EtaDx (double H, double m, double uCnoidal, double L, double K, double E); + double d2EtaDx (double H, double m, double uCnoidal, double L, double K, double E); + double d3EtaDx (double H, double m, double uCnoidal, double L, double K, double E); + int calculations (double H, double d, double T, double* mOut, double* LOut); + double U (double H, double h, double m, double kx, double ky, double T, double x, double y, double t, double z); + double W (double H, double h, double m, double kx, double ky, double T, double x, double y, double t, double z); +} + +namespace stokesVFun +{ + double A11 (double h, double k); + double A13 (double h, double k); + double A15 (double h, double k); + double A22 (double h, double k); + double A24 (double h, double k); + double A33 (double h, double k); + double A35 (double h, double k); + double A44 (double h, double k); + double A55 (double h, double k); + double B22 (double h, double k); + double B24 (double h, double k); + double B33 (double h, double k); + double B33k (double h, double k); + double B35 (double h, double k); + double B35k (double h, double k); + double B44 (double h, double k); + double B55 (double h, double k); + double B55k (double h, double k); + double C1 (double h, double k); + double C1k (double h, double k); + double C2 (double h, double k); + double C2k (double h, double k); + double C3 (double h, double k); + double C4 (double h, double k); + int StokesVNR (double H, double d, double T, double* kOut, double* LambdaOut, double* f1Out, double* f2Out ); + int StokesExtendedSolver (double H, double d, double T, double* kOut, double* LambdaOut, double* LambdaErrOut ); + double eta (double h, double kx, double ky, double lambda, double T, double x, double y, double t, double phase); + double U (double d, double kx, double ky, double lambda, double T, double x, double y, double t, double phase, double z); + double W (double d, double kx, double ky, double lambda, double T, double x, double y, double t, double phase, double z); +} + +namespace BoussinesqFun +{ + double eta (double H, double h, double x, double y, double theta, double t, double X0); + double Deta1 (double H, double h, double x, double y, double theta, double t, double X0); + double Deta2 (double H, double h, double x, double y, double theta, double t, double X0); + double Deta3 (double H, double h, double x, double y, double theta, double t, double X0); + double U (double H, double h, double x, double y, double theta, double t, double X0, double z); + double W (double H, double h, double x, double y, double theta, double t, double X0, double z); +} + +#endif diff --git a/integration/genAbs/waveAbsorption/IH_3D_3DAbsorption_InletVelocity/IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.C b/integration/genAbs/waveAbsorption/IH_3D_3DAbsorption_InletVelocity/IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.C new file mode 100644 index 0000000000000000000000000000000000000000..ee86b0517037b7c192c6c77917f5b2480b08ba78 --- /dev/null +++ b/integration/genAbs/waveAbsorption/IH_3D_3DAbsorption_InletVelocity/IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.C @@ -0,0 +1,339 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2006-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + +#include "IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.H" +#include "volFields.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" +#include "surfaceFields.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam:: +IH_3D_3DAbsorption_InletVelocityFvPatchVectorField:: +IH_3D_3DAbsorption_InletVelocityFvPatchVectorField +( + const fvPatch& p, + const DimensionedField<vector, volMesh>& iF +) +: + fixedValueFvPatchField<vector>(p, iF), + nPaddles_(1), + leftORright_(-1), + waterDepth_(-1), + initialDepthABS_(-1), + RealwaterDepth_(-1), + allCheck_(true), + waveDictName_("IHWavesDict") +{} + + +Foam:: +IH_3D_3DAbsorption_InletVelocityFvPatchVectorField:: +IH_3D_3DAbsorption_InletVelocityFvPatchVectorField +( + const IH_3D_3DAbsorption_InletVelocityFvPatchVectorField& ptf, + const fvPatch& p, + const DimensionedField<vector, volMesh>& iF, + const fvPatchFieldMapper& mapper +) +: + fixedValueFvPatchField<vector>(ptf, p, iF, mapper), + nPaddles_(ptf.nPaddles_), + leftORright_(ptf.leftORright_), + waterDepth_(ptf.waterDepth_), + initialDepthABS_(ptf.initialDepthABS_), + RealwaterDepth_(ptf.RealwaterDepth_), + allCheck_(ptf.allCheck_), + waveDictName_(ptf.waveDictName_) +{} + + +Foam:: +IH_3D_3DAbsorption_InletVelocityFvPatchVectorField:: +IH_3D_3DAbsorption_InletVelocityFvPatchVectorField +( + const fvPatch& p, + const DimensionedField<vector, volMesh>& iF, + const dictionary& dict +) +: + fixedValueFvPatchField<vector>(p, iF, dict), + nPaddles_(dict.lookupOrDefault<label>("nPaddles", 1)), + leftORright_(dict.lookupOrDefault<scalar>("leftORright",-1)), + waterDepth_(dict.lookupOrDefault<scalar>("waterDepth",-1 )), + initialDepthABS_(dict.lookupOrDefault<scalar>("initialDepthABS",-1)), + RealwaterDepth_(dict.lookupOrDefault<scalar>("RealwaterDepth",-1)), + allCheck_(dict.lookupOrDefault<bool>("allCheck", true )), + waveDictName_(dict.lookupOrDefault<word>("waveDict","IHWavesDict")) +{} + + +Foam:: +IH_3D_3DAbsorption_InletVelocityFvPatchVectorField:: +IH_3D_3DAbsorption_InletVelocityFvPatchVectorField +( + const IH_3D_3DAbsorption_InletVelocityFvPatchVectorField& ptf +) +: + fixedValueFvPatchField<vector>(ptf), + nPaddles_(ptf.nPaddles_), + leftORright_(ptf.leftORright_), + waterDepth_(ptf.waterDepth_), + initialDepthABS_(ptf.initialDepthABS_), + RealwaterDepth_(ptf.RealwaterDepth_), + allCheck_(ptf.allCheck_), + waveDictName_(ptf.waveDictName_) +{} + + +Foam:: +IH_3D_3DAbsorption_InletVelocityFvPatchVectorField:: +IH_3D_3DAbsorption_InletVelocityFvPatchVectorField +( + const IH_3D_3DAbsorption_InletVelocityFvPatchVectorField& ptf, + const DimensionedField<vector, volMesh>& iF +) +: + fixedValueFvPatchField<vector>(ptf, iF), + nPaddles_(ptf.nPaddles_), + leftORright_(ptf.leftORright_), + waterDepth_(ptf.waterDepth_), + initialDepthABS_(ptf.initialDepthABS_), + RealwaterDepth_(ptf.RealwaterDepth_), + allCheck_(ptf.allCheck_), + waveDictName_(ptf.waveDictName_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::IH_3D_3DAbsorption_InletVelocityFvPatchVectorField::updateCoeffs() +{ + if (updated()) + { + return; + } + + // Variables & constants + const vector cMin = gMin(patch().patch().localPoints()); + const vector cMax = gMax(patch().patch().localPoints()); + const vector cSpan = cMax - cMin; + const scalar zSpan = cSpan[2]; + scalar dMin = 0.0; + scalar dSpan = 0.0; + const scalarField patchD = patchDirection( cSpan, &dMin, &dSpan ); + const volScalarField& alpha = + db().lookupObject<volScalarField>(alphaName()); + const volVectorField& U = db().lookupObject<volVectorField>("U"); + const fvMesh& mesh = alpha.mesh(); + const word& patchName = this->patch().name(); + const label patchID = mesh.boundaryMesh().findPatchID(patchName); + const label nF = patch().faceCells().size(); + const scalarField alphaCell = + alpha.boundaryField()[patchID].patchInternalField(); + const vectorField UCell = U.boundaryField()[patchID].patchInternalField(); + scalarField patchUc = Foam::scalarField(nF, 0.0); + scalarField patchVc = Foam::scalarField(nF, 0.0); + scalarField patchWc = Foam::scalarField(nF, 0.0); + const scalarField patchHeight = patch().Cf().component(2); + const scalar g = 9.81; + + // Calculate Z bounds of the faces + scalarField zSup, zInf; + faceBoundsZ( &zSup, &zInf ); + + // Define dictionary + IOdictionary IHWavesDict + ( + IOobject + ( + waveDictName_, + this->db().time().constant(), + this->db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ) + ); + + scalar currTime = this->db().time().value(); + + // Check for errors - Just the first time + if (allCheck_) + { + waterDepth_ = (IHWavesDict.lookupOrDefault<scalar>("waterDepth",-1.0)); + + initialDepthABS_ = + (IHWavesDict.lookupOrDefault<scalar>("initialDepthABS",-1.0)); + + // Check if the value of nPaddles is correct for the number of columns + if (nPaddles_ < 1) + { + FatalError << "Check nPaddles value." << exit(FatalError); + } + + if ( nPaddles_ > 1 ) + { + nPaddles_ = decreaseNPaddles( nPaddles_, patchD, dMin, dSpan ); + reduce(nPaddles_, minOp<label>()); + } + } + + // Grouping part + labelList faceGroup = Foam::labelList(nF, 0); + scalarList dBreakPoints = Foam::scalarList(nPaddles_+1, dMin); + scalarList xGroup = Foam::scalarList(nPaddles_, 0.0); + scalarList yGroup = Foam::scalarList(nPaddles_, 0.0); + + for (label i=0; i<nPaddles_; i++) + { + // Breakpoints, X & Y centre of the paddles + dBreakPoints[i+1] = dMin + dSpan/(nPaddles_)*(i+1); + xGroup[i] = cMin[0] + cSpan[0]/(2.0*nPaddles_) + + cSpan[0]/(nPaddles_)*i; + yGroup[i] = cMin[1] + cSpan[1]/(2.0*nPaddles_) + + cSpan[1]/(nPaddles_)*i; + } + + forAll(patchD, patchCells) + { + for (label i=0; i<nPaddles_; i++) + { + if ( (patchD[patchCells]>=dBreakPoints[i]) + && (patchD[patchCells]<dBreakPoints[i+1]) ) + { + faceGroup[patchCells] = i+1; + continue; + } + } + } + + if (allCheck_) + { + if (RealwaterDepth_ == -1.0) + { + if (waterDepth_ == -1.0) + { + RealwaterDepth_ = + calcWL( alphaCell, faceGroup, zSpan )[0]; + + if (initialDepthABS_ !=-1.0) + { + RealwaterDepth_ = RealwaterDepth_ + + initialDepthABS_; + } + } + else if ( waterDepth_ != -1.0 ) + { + RealwaterDepth_ = waterDepth_; + } + } + + allCheck_ = false; + } + + // Calculate water measured levels + scalarList measuredLevelsABS = calcWL( alphaCell, faceGroup, zSpan ); + + if (initialDepthABS_ !=-1.0) + { + forAll(measuredLevelsABS, iterMin) + { + measuredLevelsABS[iterMin] = measuredLevelsABS[iterMin] + + initialDepthABS_; + } + } + + forAll(patchHeight, cellIndex) + { + if ( zInf[cellIndex] >= measuredLevelsABS[faceGroup[cellIndex]-1] ) + { + patchUc[cellIndex] = 0.0; + patchVc[cellIndex] = 0.0; + patchWc[cellIndex] = 0.0; + } + else + { + patchUc[cellIndex] = + ( -measuredLevelsABS[faceGroup[cellIndex]-1] + + RealwaterDepth_ ) + * sqrt(g/measuredLevelsABS[faceGroup[cellIndex]-1]); + + patchUc[cellIndex] = leftORright_ * patchUc[cellIndex]; + patchVc[cellIndex] = 0.0; + patchWc[cellIndex] = 0.0; + } + } + + const vectorField n1 = Foam::vectorField(nF, vector(1.0, 0.0, 0.0)); + const vectorField n2 = Foam::vectorField(nF, vector(0.0, 1.0, 0.0)); + const vectorField n3 = Foam::vectorField(nF, vector(0.0, 0.0, 1.0)); + + operator == (n1*patchUc + n2*patchVc + n3*patchWc); + + fixedValueFvPatchField<vector>::updateCoeffs(); +} + + +void Foam::IH_3D_3DAbsorption_InletVelocityFvPatchVectorField:: +write(Ostream& os) const +{ + fvPatchField<vector>::write(os); + + os.writeKeyword("waveDictName") << waveDictName_ << + token::END_STATEMENT << nl; + + os.writeKeyword("leftORright") << leftORright_ << + token::END_STATEMENT << nl; + + os.writeKeyword("nPaddles") << nPaddles_ << token::END_STATEMENT << nl; + + os.writeKeyword("RealwaterDepth") << RealwaterDepth_ << + token::END_STATEMENT << nl; + + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + makePatchTypeField + ( + fvPatchVectorField, + IH_3D_3DAbsorption_InletVelocityFvPatchVectorField + ); +} + + +// ************************************************************************* // diff --git a/integration/genAbs/waveAbsorption/IH_3D_3DAbsorption_InletVelocity/IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.H b/integration/genAbs/waveAbsorption/IH_3D_3DAbsorption_InletVelocity/IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.H new file mode 100644 index 0000000000000000000000000000000000000000..6fea88efe7b36b236ec287af23e82c3c2543e110 --- /dev/null +++ b/integration/genAbs/waveAbsorption/IH_3D_3DAbsorption_InletVelocity/IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.H @@ -0,0 +1,199 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2006-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::IH_3D_3DAbsorption_InletVelocityFvPatchVectorField + +Description + Wave absorption boundary condition based on shallow water theory and on a + 3D approach. Works both in 2D and 3D and for waves out of the shallow water + regime. + + Example of the boundary condition specification: + @verbatim + outlet + { + type IH_3D_3DAbsorption_InletVelocityV2; + nPaddles 1; + value uniform (0 0 0); + leftORright -1.0; + } + @endverbatim + + \heading Function object usage + \table + Property | Description | Required | Default value + type | type: IH_3D_3DAbsorption_InletVelocityV2 | yes + nPaddles | number of wavepaddles for absorption | yes | 1 + leftORright | Define location of Boundary condition: Left(1) or Right (-1) | yes | -1 + \endtable + +Note + - + +SourceFiles + IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.C + +\*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + +#ifndef IH_3D_3DAbsorption_InletVelocityFvPatchVectorField_H +#define IH_3D_3DAbsorption_InletVelocityFvPatchVectorField_H + +#include "fixedValueFvPatchFields.H" +#include "mathematicalConstants.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +/*---------------------------------------------------------------------------*\ + Class IH_3D_3DAbsorption_InletVelocityFvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class IH_3D_3DAbsorption_InletVelocityFvPatchVectorField +: + public fixedValueFvPatchVectorField +{ + // Private data + + //- Number of paddles + label nPaddles_; + + //Aborption: Left(1) or Right(-1) + scalar leftORright_; + + //- Theoretical Water depth (meters) + scalar waterDepth_; + + //- Numerical Water depth (meters) + scalar RealwaterDepth_; + + //- BC check and read values just the first time + bool allCheck_; + + //- Initial Depth on BC for absorption (meters) + scalar initialDepthABS_; + + //- Dictionary name + word waveDictName_; + +public: + + //- Runtime type information + TypeName("IH_3D_3DAbsorption_InletVelocityV2"); + + + // Constructors + + //- Construct from patch and internal field + IH_3D_3DAbsorption_InletVelocityFvPatchVectorField + ( + const fvPatch&, + const DimensionedField<vector, volMesh>& + ); + + //- Construct from patch, internal field and dictionary + IH_3D_3DAbsorption_InletVelocityFvPatchVectorField + ( + const fvPatch&, + const DimensionedField<vector, volMesh>&, + const dictionary& + ); + + //- Construct by mapping given + // IH_3D_3DAbsorption_InletVelocityFvPatchVectorField + // onto a new patch + IH_3D_3DAbsorption_InletVelocityFvPatchVectorField + ( + const IH_3D_3DAbsorption_InletVelocityFvPatchVectorField&, + const fvPatch&, + const DimensionedField<vector, volMesh>&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + IH_3D_3DAbsorption_InletVelocityFvPatchVectorField + ( + const IH_3D_3DAbsorption_InletVelocityFvPatchVectorField& + ); + + //- Construct and return a clone + virtual tmp<fvPatchVectorField> clone() const + { + return tmp<fvPatchVectorField> + ( + new IH_3D_3DAbsorption_InletVelocityFvPatchVectorField(*this) + ); + } + + //- Construct as copy setting internal field reference + IH_3D_3DAbsorption_InletVelocityFvPatchVectorField + ( + const IH_3D_3DAbsorption_InletVelocityFvPatchVectorField&, + const DimensionedField<vector, volMesh>& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp<fvPatchVectorField> clone + ( + const DimensionedField<vector, volMesh>& iF + ) const + { + return tmp<fvPatchVectorField> + ( + new IH_3D_3DAbsorption_InletVelocityFvPatchVectorField + (*this, iF) + ); + } + + + // Member functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Write + virtual void write(Ostream&) const; + + // Other common member functions + #include "memberFun.H" + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/integration/genAbs/waveAbsorption/Make/files b/integration/genAbs/waveAbsorption/Make/files new file mode 100644 index 0000000000000000000000000000000000000000..7f9bff3dfbfc2fd0a49b3addf5ea87e906a443ce --- /dev/null +++ b/integration/genAbs/waveAbsorption/Make/files @@ -0,0 +1,3 @@ +IH_3D_3DAbsorption_InletVelocity/IH_3D_3DAbsorption_InletVelocityFvPatchVectorField.C + +LIB = $(FOAM_USER_LIBBIN)/libIHwaveAbsorptionV2esi diff --git a/integration/genAbs/waveAbsorption/Make/options b/integration/genAbs/waveAbsorption/Make/options new file mode 100644 index 0000000000000000000000000000000000000000..be8b9342202495df469eb38e5191a153e27f60f9 --- /dev/null +++ b/integration/genAbs/waveAbsorption/Make/options @@ -0,0 +1,7 @@ +EXE_INC = \ + -DOFVERSION=$(OF_VERSION) \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I../common + +LIB_LIBS = \ + -lfiniteVolume diff --git a/integration/genAbs/waveAbsorption/localMake b/integration/genAbs/waveAbsorption/localMake new file mode 100755 index 0000000000000000000000000000000000000000..aedb1c100c89e09e8db4e48b788933e6bd68453d --- /dev/null +++ b/integration/genAbs/waveAbsorption/localMake @@ -0,0 +1,25 @@ +#!/bin/bash + +if [ $WM_PROJECT == "foam" ]; then + ofversion=`echo $WM_PROJECT_VERSION | sed -e 's/\.x/-9/' -e 's/\./\'$'\n/g' -e 's/-/\'$'\n/g' | grep "[0-9]" | head -2 | tr -d '\n'` +else + ofversion=`echo $WM_PROJECT_VERSION"-0" | sed -e 's/\.x/-9/' -e 's/\./\'$'\n/g' -e 's/-/\'$'\n/g' -e 's/+/\'$'\n/g' -e 's/v/\'$'\n/g' | grep "[0-9]" | head -3 | tr -d '\n'` +fi + +#IHC_dbg +echo $ofversion +#---- + +export OF_VERSION=$ofversion + +wclean + +wmake libso + +if (( $? )) ; then + echo "\n\nIH Wave absorption boundary conditions (V2.0) compilation failed" + exit; else + echo -e "\n\nIH Wave absorption boundary conditions (V2.0) compiled successfully for $WM_PROJECT $WM_PROJECT_VERSION\n\n\n"; +fi + +wclean diff --git a/integration/genAbs/waveGeneration/IH_Waves_InletAlpha/IH_Waves_InletAlphaFvPatchScalarField.C b/integration/genAbs/waveGeneration/IH_Waves_InletAlpha/IH_Waves_InletAlphaFvPatchScalarField.C new file mode 100644 index 0000000000000000000000000000000000000000..51052f0b79c0db68bf2f46e3f58fba62f95db81a --- /dev/null +++ b/integration/genAbs/waveGeneration/IH_Waves_InletAlpha/IH_Waves_InletAlphaFvPatchScalarField.C @@ -0,0 +1,460 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2006-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + +#include "IH_Waves_InletAlphaFvPatchScalarField.H" +#include "volFields.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" +#include "surfaceFields.H" +#include "Random.H" + +#include "waveFun.H" + +//C++ +#include <stdio.h> +#include <unistd.h> +#include <errno.h> +#include "SquareMatrix.H" +#include "vector.H" +#include "Matrix.H" + +#include "mpi.h" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam:: +IH_Waves_InletAlphaFvPatchScalarField:: +IH_Waves_InletAlphaFvPatchScalarField +( + const fvPatch& p, + const DimensionedField<scalar, volMesh>& iF +) +: + fixedValueFvPatchField<scalar>(p, iF), + wavePeriod_(-1), + waveHeight_(-1), + waveLength_(-1), + waterDepth_(-1), + initialDepth_(-1), + wavePhase_(3.0*PI()/2.0), + lambdaStokesV_(-1), + mCnoidal_(-1), + nPaddles_(1), + tSmooth_(-1), + waveDictName_("IHWavesDict"), + waveType_("aaa"), + waveTheory_("aaa"), + allCheck_(true), + waveDir_(0), + RealwaterDepth_(-1) +{} + + +Foam:: +IH_Waves_InletAlphaFvPatchScalarField:: +IH_Waves_InletAlphaFvPatchScalarField +( + const IH_Waves_InletAlphaFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField<scalar, volMesh>& iF, + const fvPatchFieldMapper& mapper +) +: + fixedValueFvPatchField<scalar>(ptf, p, iF, mapper), + wavePeriod_(ptf.wavePeriod_), + waveHeight_(ptf.waveHeight_), + waveLength_(ptf.waveLength_), + waterDepth_(ptf.waterDepth_), + initialDepth_(ptf.initialDepth_), + wavePhase_(ptf.wavePhase_), + lambdaStokesV_(ptf.lambdaStokesV_), + mCnoidal_(ptf.mCnoidal_), + nPaddles_(ptf.nPaddles_), + tSmooth_(ptf.tSmooth_), + waveDictName_(ptf.waveDictName_), + waveType_(ptf.waveType_), + waveTheory_(ptf.waveTheory_), + allCheck_(ptf.allCheck_), + waveDir_(ptf.waveDir_), + RealwaterDepth_(ptf.RealwaterDepth_) +{} + + +Foam:: +IH_Waves_InletAlphaFvPatchScalarField:: +IH_Waves_InletAlphaFvPatchScalarField +( + const fvPatch& p, + const DimensionedField<scalar, volMesh>& iF, + const dictionary& dict +) +: + fixedValueFvPatchField<scalar>(p,iF,dict), + wavePeriod_(dict.lookupOrDefault<scalar>("wavePeriod",-1)), + waveHeight_(dict.lookupOrDefault<scalar>("waveHeight",-1)), + waveLength_(dict.lookupOrDefault<scalar>("waveLength",-1)), + waterDepth_(dict.lookupOrDefault<scalar>("waterDepth",-1)), + initialDepth_(dict.lookupOrDefault<scalar>("initialDepth",-1)), + wavePhase_(dict.lookupOrDefault<scalar>("wavePhase",3.0*PI()/2.0)), + lambdaStokesV_(dict.lookupOrDefault<scalar>("lambdaStokesV",-1)), + mCnoidal_(dict.lookupOrDefault<scalar>("mCnoidal",-1)), + nPaddles_(dict.lookupOrDefault<label>("nPaddles",1)), + tSmooth_(dict.lookupOrDefault<scalar>("tSmooth",-1)), + waveDictName_(dict.lookupOrDefault<word>("waveDict","IHWavesDict")), + waveType_(dict.lookupOrDefault<word>("waveType","aaa")), + waveTheory_(dict.lookupOrDefault<word>("waveTheory","aaa")), + allCheck_(dict.lookupOrDefault<bool>("allCheck",true)), + waveDir_(dict.lookupOrDefault<scalar>("waveDir",0)), + RealwaterDepth_(dict.lookupOrDefault<scalar>("RealwaterDepth",-1)) +{} + + +Foam:: +IH_Waves_InletAlphaFvPatchScalarField:: +IH_Waves_InletAlphaFvPatchScalarField +( + const IH_Waves_InletAlphaFvPatchScalarField& ptf +) +: + fixedValueFvPatchField<scalar>(ptf), + wavePeriod_(ptf.wavePeriod_), + waveHeight_(ptf.waveHeight_), + waveLength_(ptf.waveLength_), + waterDepth_(ptf.waterDepth_), + initialDepth_(ptf.initialDepth_), + wavePhase_(ptf.wavePhase_), + lambdaStokesV_(ptf.lambdaStokesV_), + mCnoidal_(ptf.mCnoidal_), + nPaddles_(ptf.nPaddles_), + tSmooth_(ptf.tSmooth_), + waveDictName_(ptf.waveDictName_), + waveType_(ptf.waveType_), + waveTheory_(ptf.waveTheory_), + allCheck_(ptf.allCheck_), + waveDir_(ptf.waveDir_), + RealwaterDepth_(ptf.RealwaterDepth_) +{} + + +Foam:: +IH_Waves_InletAlphaFvPatchScalarField:: +IH_Waves_InletAlphaFvPatchScalarField +( + const IH_Waves_InletAlphaFvPatchScalarField& ptf, + const DimensionedField<scalar, volMesh>& iF +) +: + fixedValueFvPatchField<scalar>(ptf, iF), + wavePeriod_(ptf.wavePeriod_), + waveHeight_(ptf.waveHeight_), + waveLength_(ptf.waveLength_), + waterDepth_(ptf.waterDepth_), + initialDepth_(ptf.initialDepth_), + wavePhase_(ptf.wavePhase_), + lambdaStokesV_(ptf.lambdaStokesV_), + mCnoidal_(ptf.mCnoidal_), + nPaddles_(ptf.nPaddles_), + tSmooth_(ptf.tSmooth_), + waveDictName_(ptf.waveDictName_), + waveType_(ptf.waveType_), + waveTheory_(ptf.waveTheory_), + allCheck_(ptf.allCheck_), + waveDir_(ptf.waveDir_), + RealwaterDepth_(ptf.RealwaterDepth_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::IH_Waves_InletAlphaFvPatchScalarField::updateCoeffs() +{ + if (updated()) + { + return; + } + + // Aux. values + scalar auxiliar = 0; + scalar auxiliarTotal = 0; + scalar X0 = 0; + scalarField patchXsolit; + + const vector cMin = gMin(patch().patch().localPoints()); + const vector cMax = gMax(patch().patch().localPoints()); + const vector cSpan = cMax - cMin; + const scalar zSpan = cSpan[2]; + + scalar dMin = 0.0; + scalar dSpan = 0.0; + const scalarField patchD = patchDirection( cSpan, &dMin, &dSpan ); + + // Variables & constants + const volScalarField& alpha = + db().lookupObject<volScalarField>(alphaName()); + const volVectorField& U = db().lookupObject<volVectorField>("U"); + const fvMesh& mesh = alpha.mesh(); + const word& patchName = this->patch().name(); + const label patchID = mesh.boundaryMesh().findPatchID(patchName); + const label nF = patch().faceCells().size(); + labelList cellGroup = Foam::labelList(nF, 1); + const scalarField alphaCell = + alpha.boundaryField()[patchID].patchInternalField(); + const vectorField UCell = U.boundaryField()[patchID].patchInternalField(); + scalarField patchVOF = alphaCell; + const labelList celdas = patch().faceCells(); + const scalarField patchHeight = patch().Cf().component(2); + + // Calculate Z bounds of the faces + scalarField zSup, zInf; + faceBoundsZ( &zSup, &zInf ); + + // Wave variables + scalar waveOmega; + scalarList waveOmegas; + scalar waveK; + scalarList waveKs; + scalar waveAngle; + scalar waveKx; + scalarList waveKxs; + scalar waveKy; + + // Define dictionary + IOdictionary IHWavesDict + ( + IOobject + ( + waveDictName_, + this->db().time().constant(), + this->db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ) + ); + + scalar currTime = this->db().time().value(); + + // Check for errors - Just the first time + if (allCheck_) + { + waveType_ = (IHWavesDict.lookupOrDefault<word>("waveType","aaa")); + tSmooth_ = (IHWavesDict.lookupOrDefault<scalar>("tSmooth",-1.0)); + waterDepth_ = + (IHWavesDict.lookupOrDefault<scalar>("waterDepth",-1.0)); + initialDepth_ = + (IHWavesDict.lookupOrDefault<scalar>("initialDepth",-1.0)); + nPaddles_ = (IHWavesDict.lookupOrDefault<label>("nPaddles",1)); + } + + // Grouping part + scalarList dBreakPoints = Foam::scalarList(nPaddles_+1, dMin); + scalarList xGroup = Foam::scalarList(nPaddles_, 0.0); + scalarList yGroup = Foam::scalarList(nPaddles_, 0.0); + + for (label i=0; i<nPaddles_; i++) + { + // Breakpoints, X & Y centre of the paddles + dBreakPoints[i+1] = dMin + dSpan/(nPaddles_)*(i+1); + xGroup[i] = cMin[0] + cSpan[0]/(2.0*nPaddles_) + + cSpan[0]/(nPaddles_)*i; + yGroup[i] = cMin[1] + cSpan[1]/(2.0*nPaddles_) + + cSpan[1]/(nPaddles_)*i; + } + + forAll(patchD, patchCells) + { + for (label i=0; i<nPaddles_; i++) + { + if ( (patchD[patchCells]>=dBreakPoints[i]) + && (patchD[patchCells]<dBreakPoints[i+1]) ) + { + cellGroup[patchCells] = i+1; + continue; + } + } + } + + // Check for errors - Just the first time + if (allCheck_) + { + if (RealwaterDepth_ == -1.0) + { + if (waterDepth_ == -1.0) + { + RealwaterDepth_ = + calcWL( alphaCell, cellGroup, zSpan )[0]; + if (initialDepth_ !=-1.0) + { + RealwaterDepth_ = RealwaterDepth_ + + initialDepth_; + } + } + else if ( waterDepth_ != -1.0 ) + { + RealwaterDepth_ = waterDepth_; + } + } + } + + if (allCheck_) + { + waveTheory_ = (IHWavesDict.lookupOrDefault<word>("waveTheory","aaa")); + waveHeight_ = (IHWavesDict.lookupOrDefault<scalar>("waveHeight",-1)); + wavePeriod_ = (IHWavesDict.lookupOrDefault<scalar>("wavePeriod",-1)); + waveDir_ = (IHWavesDict.lookupOrDefault<scalar>("waveDir",0)); + wavePhase_ = + (IHWavesDict.lookupOrDefault<scalar>("wavePhase",3.0*PI()/2.0)); + + if ( waveType_ == "regular" ) + { + waveLength_ = StokesIFun::waveLength (RealwaterDepth_, wavePeriod_); + } + } + + if ( waveType_ == "regular" ) + { + waveOmega = (2.0*PI())/wavePeriod_; + waveK = 2.0*PI()/waveLength_; + + waveAngle = waveDir_*PI()/180.0; + waveKx = waveK*cos(waveAngle); + waveKy = waveK*sin(waveAngle); + } + else if ( waveType_ == "solitary" ) + { + waveAngle = waveDir_*PI()/180.0; + patchXsolit = patch().Cf().component(0)*cos(waveAngle) + + patch().Cf().component(1)*sin(waveAngle); + X0 = gMin(patchXsolit); + } + + scalar timeMult = 1.0; + + if ( tSmooth_ > 0 && currTime < tSmooth_ ) + { + timeMult = currTime/tSmooth_; + } + + if (allCheck_) + { + if ( waveType_ == "regular" ) + { + #include "checkInputErrorsRegular.H" + } + else if ( waveType_ == "solitary" ) + { + #include "checkInputErrorsSolitary.H" + } + else + { + FatalError << "Wave type not supported, use:\n" + << "regular, solitary" << exit(FatalError); + } + allCheck_ = false; + } + + scalarList calculatedLevel (nPaddles_,0.0); + + if ( waveType_ == "regular" ) + { + if (waveDir_ == 0) + { + #include "calculatedLevelRegularNormal.H" + } + else + { + #include "calculatedLevelRegular.H" + } + } + else if ( waveType_ == "solitary" ) + { + #include "calculatedLevelSolitary.H" + } + + forAll(patchHeight, cellIndex) + { + auxiliarTotal = 0; + auxiliar = 0; + + if (zSup[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) + { + patchVOF[cellIndex] = 1.0; + } + else if ( zInf[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + { + patchVOF[cellIndex] = 0.0; + } + else + { + auxiliar = calculatedLevel[cellGroup[cellIndex]-1] + - zInf[cellIndex]; + auxiliarTotal = zSup[cellIndex]-zInf[cellIndex]; + auxiliarTotal = auxiliar/auxiliarTotal; + patchVOF[cellIndex] = auxiliarTotal; + } + } + + operator == (patchVOF); + + fixedValueFvPatchField<scalar>::updateCoeffs(); +} + + +void Foam::IH_Waves_InletAlphaFvPatchScalarField::write(Ostream& os) const +{ + fvPatchField<scalar>::write(os); + + os.writeKeyword("waveDictName") << waveDictName_ << + token::END_STATEMENT << nl; + + os.writeKeyword("RealwaterDepth") << RealwaterDepth_ << + token::END_STATEMENT << nl; + + os.writeKeyword("initialDepth") << initialDepth_ << + token::END_STATEMENT << nl; + + writeEntry("value", os); +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + makePatchTypeField + ( + fvPatchScalarField, + IH_Waves_InletAlphaFvPatchScalarField + ); + +} + + +// ************************************************************************* // diff --git a/integration/genAbs/waveGeneration/IH_Waves_InletAlpha/IH_Waves_InletAlphaFvPatchScalarField.H b/integration/genAbs/waveGeneration/IH_Waves_InletAlpha/IH_Waves_InletAlphaFvPatchScalarField.H new file mode 100644 index 0000000000000000000000000000000000000000..0e181a0fff9845213db6f3f9341a34da26dd3542 --- /dev/null +++ b/integration/genAbs/waveGeneration/IH_Waves_InletAlpha/IH_Waves_InletAlphaFvPatchScalarField.H @@ -0,0 +1,226 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2006-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::IH_Waves_InletAlphaFvPatchScalarField + +Description + Describes a volumetric/mass flow normal Scalar boundary condition by its + magnitude as an integral over its area. + + The basis of the patch (volumetric or mass) is determined by the + dimensions of the flux, phi. + The current density is used to correct the Alpha when applying the + mass basis. + + Example of the boundary condition specification: + @verbatim + inlet + { + type IH_Waves_InletAlphaV2; + waveDict IHWavesDict; + } + @endverbatim + + \heading Function object usage + \table + Property | Description | Required | Default value + type | type: IH_Waves_InletAlphaV2 | yes + waveDict | Dictionary where variables for generation/absorption are defined | yes | IHWavesDict + \endtable + +Note + - + +SourceFiles + IH_Waves_InletAlphaFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + +#ifndef IH_Waves_InletAlphaFvPatchScalarField_H +#define IH_Waves_InletAlphaFvPatchScalarField_H + +#include "fixedValueFvPatchFields.H" +#include "mathematicalConstants.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +/*---------------------------------------------------------------------------*\ + Class IH_Waves_InletAlphaFvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class IH_Waves_InletAlphaFvPatchScalarField +: + public fixedValueFvPatchScalarField +{ + // Private data + + //- Wave period (seconds) + scalar wavePeriod_; + + //- Wave height (meters) + scalar waveHeight_; + + //- Wave length (meters) + scalar waveLength_; + + //- Water depth (meters) + scalar waterDepth_; + + //- Wave phase (radians) + scalar wavePhase_; + + //- Lambda StokesV parameter + scalar lambdaStokesV_; + + //- m Cnoidal parameter + scalar mCnoidal_; + + //- Number of different paddles (for absorption) + label nPaddles_; + + //- Smoothing factor for the wave(s) (seconds) + scalar tSmooth_; + + //- Dictionary name + word waveDictName_; + + //- Regular or Irregular wave(s) + word waveType_; + + //- Name of the theory + word waveTheory_; + + //- BC check and read values just the first time + bool allCheck_; + + //- Direction of propagation (degrees, from X axis) + scalar waveDir_; + + //- Numerical Water depth (meters) + scalar RealwaterDepth_; + + //- Initial Depth on BC for wave generation (meters) + scalar initialDepth_; + +public: + + //- Runtime type information + TypeName("IH_Waves_InletAlphaV2"); + + + // Constructors + + //- Construct from patch and internal field + IH_Waves_InletAlphaFvPatchScalarField + ( + const fvPatch&, + const DimensionedField<scalar, volMesh>& + ); + + //- Construct from patch, internal field and dictionary + IH_Waves_InletAlphaFvPatchScalarField + ( + const fvPatch&, + const DimensionedField<scalar, volMesh>&, + const dictionary& + ); + + //- Construct by mapping given + // IH_Waves_InletAlphaFvPatchScalarField + // onto a new patch + IH_Waves_InletAlphaFvPatchScalarField + ( + const IH_Waves_InletAlphaFvPatchScalarField&, + const fvPatch&, + const DimensionedField<scalar, volMesh>&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + IH_Waves_InletAlphaFvPatchScalarField + ( + const IH_Waves_InletAlphaFvPatchScalarField& + ); + + //- Construct and return a clone + virtual tmp<fvPatchScalarField> clone() const + { + return tmp<fvPatchScalarField> + ( + new IH_Waves_InletAlphaFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + IH_Waves_InletAlphaFvPatchScalarField + ( + const IH_Waves_InletAlphaFvPatchScalarField&, + const DimensionedField<scalar, volMesh>& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp<fvPatchScalarField> clone + ( + const DimensionedField<scalar, volMesh>& iF + ) const + { + return tmp<fvPatchScalarField> + ( + new IH_Waves_InletAlphaFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Write + virtual void write(Ostream&) const; + + // Other common member functions + #include "memberFun.H" + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/IH_Waves_InletVelocityFvPatchVectorField.C b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/IH_Waves_InletVelocityFvPatchVectorField.C new file mode 100644 index 0000000000000000000000000000000000000000..a25bf8c254e6019c1ab85916de658c09a71b9091 --- /dev/null +++ b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/IH_Waves_InletVelocityFvPatchVectorField.C @@ -0,0 +1,495 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2006-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + +#include "IH_Waves_InletVelocityFvPatchVectorField.H" +#include "volFields.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" +#include "surfaceFields.H" + +#include "waveFun.H" + +#include <stdio.h> +#include <unistd.h> +#include <errno.h> +#include "SquareMatrix.H" +#include "vector.H" +#include "Matrix.H" + +#include "mpi.h" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam:: +IH_Waves_InletVelocityFvPatchVectorField:: +IH_Waves_InletVelocityFvPatchVectorField +( + const fvPatch& p, + const DimensionedField<vector, volMesh>& iF +) +: + fixedValueFvPatchField<vector>(p, iF), + wavePeriod_(-1), + waveHeight_(-1), + waveLength_(-1), + waterDepth_(-1), + initialDepth_(-1), + RealwaterDepth_(-1), + wavePhase_(3.0*PI()/2.0), + lambdaStokesV_(-1), + mCnoidal_(-1), + genAbs_(false), + nPaddles_(1), + tSmooth_(-1), + leftORright_(-1), + waveDictName_("IHWavesDict"), + waveType_("aaa"), + waveTheory_("aaa"), + allCheck_(true), + waveDir_(0) +{} + + +Foam:: +IH_Waves_InletVelocityFvPatchVectorField:: +IH_Waves_InletVelocityFvPatchVectorField +( + const IH_Waves_InletVelocityFvPatchVectorField& ptf, + const fvPatch& p, + const DimensionedField<vector, volMesh>& iF, + const fvPatchFieldMapper& mapper +) +: + fixedValueFvPatchField<vector>(ptf, p, iF, mapper), + wavePeriod_(ptf.wavePeriod_), + waveHeight_(ptf.waveHeight_), + waveLength_(ptf.waveLength_), + waterDepth_(ptf.waterDepth_), + initialDepth_(ptf.initialDepth_), + RealwaterDepth_(ptf.RealwaterDepth_), + wavePhase_(ptf.wavePhase_), + lambdaStokesV_(ptf.lambdaStokesV_), + mCnoidal_(ptf.mCnoidal_), + genAbs_(ptf.genAbs_), + nPaddles_(ptf.nPaddles_), + tSmooth_(ptf.tSmooth_), + leftORright_(ptf.leftORright_), + waveDictName_(ptf.waveDictName_), + waveType_(ptf.waveType_), + waveTheory_(ptf.waveTheory_), + allCheck_(ptf.allCheck_), + waveDir_(ptf.waveDir_) +{} + + +Foam:: +IH_Waves_InletVelocityFvPatchVectorField:: +IH_Waves_InletVelocityFvPatchVectorField +( + const fvPatch& p, + const DimensionedField<vector, volMesh>& iF, + const dictionary& dict +) +: + fixedValueFvPatchField<vector>(p, iF, dict), + wavePeriod_(dict.lookupOrDefault<scalar>("wavePeriod",-1)), + waveHeight_(dict.lookupOrDefault<scalar>("waveHeight",-1)), + waveLength_(dict.lookupOrDefault<scalar>("waveLength",-1)), + waterDepth_(dict.lookupOrDefault<scalar>("waterDepth",-1)), + initialDepth_(dict.lookupOrDefault<scalar>("initialDepth",-1)), + RealwaterDepth_(dict.lookupOrDefault<scalar>("RealwaterDepth",-1)), + wavePhase_(dict.lookupOrDefault<scalar>("wavePhase",3.0*PI()/2.0)), + lambdaStokesV_(dict.lookupOrDefault<scalar>("lambdaStokesV",-1)), + mCnoidal_(dict.lookupOrDefault<scalar>("mCnoidal",-1)), + genAbs_(dict.lookupOrDefault<bool>("genAbs",false)), + nPaddles_(dict.lookupOrDefault<label>("nPaddles",1)), + tSmooth_(dict.lookupOrDefault<scalar>("tSmooth",-1)), + leftORright_(dict.lookupOrDefault<scalar>("leftORright",-1)), + waveDictName_(dict.lookupOrDefault<word>("waveDict","IHWavesDict")), + waveType_(dict.lookupOrDefault<word>("waveType","aaa")), + waveTheory_(dict.lookupOrDefault<word>("waveTheory","aaa")), + allCheck_(dict.lookupOrDefault<bool>("allCheck",true)), + waveDir_(dict.lookupOrDefault<scalar>("waveDir",0)) +{} + + +Foam:: +IH_Waves_InletVelocityFvPatchVectorField:: +IH_Waves_InletVelocityFvPatchVectorField +( + const IH_Waves_InletVelocityFvPatchVectorField& ptf +) +: + fixedValueFvPatchField<vector>(ptf), + wavePeriod_(ptf.wavePeriod_), + waveHeight_(ptf.waveHeight_), + waveLength_(ptf.waveLength_), + waterDepth_(ptf.waterDepth_), + initialDepth_(ptf.initialDepth_), + RealwaterDepth_(ptf.RealwaterDepth_), + wavePhase_(ptf.wavePhase_), + lambdaStokesV_(ptf.lambdaStokesV_), + mCnoidal_(ptf.mCnoidal_), + genAbs_(ptf.genAbs_), + nPaddles_(ptf.nPaddles_), + tSmooth_(ptf.tSmooth_), + leftORright_(ptf.leftORright_), + waveDictName_(ptf.waveDictName_), + waveType_(ptf.waveType_), + waveTheory_(ptf.waveTheory_), + allCheck_(ptf.allCheck_), + waveDir_(ptf.waveDir_) +{} + + +Foam:: +IH_Waves_InletVelocityFvPatchVectorField:: +IH_Waves_InletVelocityFvPatchVectorField +( + const IH_Waves_InletVelocityFvPatchVectorField& ptf, + const DimensionedField<vector, volMesh>& iF +) +: + fixedValueFvPatchField<vector>(ptf, iF), + wavePeriod_(ptf.wavePeriod_), + waveHeight_(ptf.waveHeight_), + waveLength_(ptf.waveLength_), + waterDepth_(ptf.waterDepth_), + initialDepth_(ptf.initialDepth_), + RealwaterDepth_(ptf.RealwaterDepth_), + wavePhase_(ptf.wavePhase_), + lambdaStokesV_(ptf.lambdaStokesV_), + mCnoidal_(ptf.mCnoidal_), + genAbs_(ptf.genAbs_), + nPaddles_(ptf.nPaddles_), + tSmooth_(ptf.tSmooth_), + leftORright_(ptf.leftORright_), + waveDictName_(ptf.waveDictName_), + waveType_(ptf.waveType_), + waveTheory_(ptf.waveTheory_), + allCheck_(ptf.allCheck_), + waveDir_(ptf.waveDir_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::IH_Waves_InletVelocityFvPatchVectorField::updateCoeffs() +{ + if (updated()) + { + return; + } + + // Auxiliar variables + scalar auxiliar = 0; + scalar auxiliarTotal = 0; + scalar X0 = 0; + scalarField patchXsolit; + + // Variables stream function + scalar celerity = 0; + + // 3D Variables + const vector cMin = gMin(patch().patch().localPoints()); + const vector cMax = gMax(patch().patch().localPoints()); + const vector cSpan = cMax - cMin; + const scalar zSpan = cSpan[2]; + + scalar dMin = 0.0; + scalar dSpan = 0.0; + const scalarField patchD = patchDirection( cSpan, &dMin, &dSpan ); + + // Variables & constants + const volScalarField& alpha = + db().lookupObject<volScalarField>(alphaName()); + const fvMesh& mesh = alpha.mesh(); + const word& patchName = this->patch().name(); + const label patchID = mesh.boundaryMesh().findPatchID(patchName); + const label nF = patch().faceCells().size(); + labelList cellGroup = Foam::labelList(nF, 1); + const scalarField alphaCell = + alpha.boundaryField()[patchID].patchInternalField(); + scalarField patchU = Foam::scalarField(nF, 0.0); + scalarField patchUABS = Foam::scalarField(nF, 0.0); + scalarField patchV = Foam::scalarField(nF, 0.0); + scalarField patchVABS = Foam::scalarField(nF, 0.0); + scalarField patchW = Foam::scalarField(nF, 0.0); + const labelList celdas = patch().faceCells(); + const scalarField patchHeight = patch().Cf().component(2); + const scalar g = 9.81; + + // Calculate Z bounds of the faces + scalarField zSup, zInf; + faceBoundsZ( &zSup, &zInf ); + + // Waves variables + scalar waveOmega; + scalar waveK; + scalar waveAngle; + scalar waveKx; + scalar waveKy; + + // Define dictionary + IOdictionary IHWavesDict + ( + IOobject + ( + waveDictName_, + this->db().time().constant(), + this->db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ) + ); + + scalar currTime = this->db().time().value(); + + // Check for errors - Just the first time + if (allCheck_) + { + waveType_ = (IHWavesDict.lookupOrDefault<word>("waveType","aaa")); + tSmooth_ = (IHWavesDict.lookupOrDefault<scalar>("tSmooth",-1.0)); + genAbs_ = (IHWavesDict.lookupOrDefault<bool>("genAbs",false)); + waterDepth_ = (IHWavesDict.lookupOrDefault<scalar>("waterDepth",-1)); + initialDepth_ = + (IHWavesDict.lookupOrDefault<scalar>("initialDepth",-1)); + nPaddles_ = (IHWavesDict.lookupOrDefault<label>("nPaddles",1)); + } + + // Grouping part + scalarList dBreakPoints = Foam::scalarList(nPaddles_+1, dMin); + scalarList xGroup = Foam::scalarList(nPaddles_, 0.0); + scalarList yGroup = Foam::scalarList(nPaddles_, 0.0); + + for (label i=0; i<nPaddles_; i++) + { + // Breakpoints, X & Y centre of the paddles + dBreakPoints[i+1] = dMin + dSpan/(nPaddles_)*(i+1); + xGroup[i] = cMin[0] + cSpan[0]/(2.0*nPaddles_) + + cSpan[0]/(nPaddles_)*i; + yGroup[i] = cMin[1] + cSpan[1]/(2.0*nPaddles_) + + cSpan[1]/(nPaddles_)*i; + } + + forAll(patchD, patchCells) + { + for (label i=0; i<nPaddles_; i++) + { + if ( (patchD[patchCells]>=dBreakPoints[i]) + && (patchD[patchCells]<dBreakPoints[i+1]) ) + { + cellGroup[patchCells] = i+1; // Group of each face + continue; + } + } + } + + // Check for errors - Just the first time + if (allCheck_) + { + if (RealwaterDepth_ == -1.0) + { + if (waterDepth_ == -1.0) + { + RealwaterDepth_ = + calcWL( alphaCell, cellGroup, zSpan )[0]; + if (initialDepth_ !=-1.0) + { + RealwaterDepth_ = RealwaterDepth_ + + initialDepth_; + } + } + else if ( waterDepth_ != -1.0 ) + { + RealwaterDepth_ = waterDepth_; + } + } + } + + if (allCheck_) + { + waveTheory_ = (IHWavesDict.lookupOrDefault<word>("waveTheory","aaa")); + waveHeight_ = (IHWavesDict.lookupOrDefault<scalar>("waveHeight",-1)); + wavePeriod_ = (IHWavesDict.lookupOrDefault<scalar>("wavePeriod",-1)); + waveDir_ = (IHWavesDict.lookupOrDefault<scalar>("waveDir",0)); + genAbs_ = (IHWavesDict.lookupOrDefault<bool>("genAbs",false)); + wavePhase_ = + (IHWavesDict.lookupOrDefault<scalar>("wavePhase",3.0*PI()/2.0)); + + if ( waveType_ == "regular" ) + { + waveLength_ = StokesIFun::waveLength ( + RealwaterDepth_, wavePeriod_ ); + } + } + + if ( waveType_ == "regular" ) + { + waveOmega = (2.0*PI())/wavePeriod_; + waveK = 2.0*PI()/waveLength_; + + celerity = waveLength_/wavePeriod_; + + waveAngle = waveDir_*PI()/180.0; + waveKx = waveK*cos(waveAngle); + waveKy = waveK*sin(waveAngle); + } + else if ( waveType_ == "solitary" ) + { + waveAngle = waveDir_*PI()/180.0; + patchXsolit = patch().Cf().component(0)*cos(waveAngle) + + patch().Cf().component(1)*sin(waveAngle); + X0 = gMin(patchXsolit); + } + + scalar timeMult = 1.0; + + if ( tSmooth_ > 0 && currTime < tSmooth_ ) + { + timeMult = currTime/tSmooth_; + } + + if (allCheck_) + { + if ( waveType_ == "regular" ) + { + #include "checkInputErrorsRegular.H" + } + else if ( waveType_ == "solitary" ) + { + #include "checkInputErrorsSolitary.H" + } + else + { + FatalError << "Wave type not supported, use:\n" + << "regular, solitary."<< exit(FatalError); + } + + allCheck_ = false; + } + + scalarList calculatedLevel (nPaddles_,0.0); + + if ( waveType_ == "regular" ) + { + if (waveDir_ == 0) + { + #include "calculatedLevelRegularNormal.H" + } + else + { + #include "calculatedLevelRegular.H" + } + } + else if ( waveType_ == "solitary" ) + { + #include "calculatedLevelSolitary.H" + } + + scalarList measuredLevels (nPaddles_,0.0); + forAll(measuredLevels, iterMin) + { + measuredLevels[iterMin] = RealwaterDepth_; + } + + scalarList measuredLevelsGENAB = calcWL( alphaCell, cellGroup, zSpan ); + + if (initialDepth_ !=-1.0) + { + forAll(measuredLevelsGENAB, iterMin) + { + measuredLevelsGENAB[iterMin] = measuredLevelsGENAB[iterMin] + + initialDepth_; + } + } + + // Define heights as minimum of calculatedLevel and measuredLevels + scalarList heights (nPaddles_,0.0); + forAll(heights, iterMin) + { + heights[iterMin] = + min(calculatedLevel[iterMin],measuredLevels[iterMin]); + } + + forAll(patchHeight, cellIndex) + { + #include "velocityProfile.H" + + patchU[cellIndex] = patchU[cellIndex]*alphaCell[cellIndex]; + patchV[cellIndex] = patchV[cellIndex]*alphaCell[cellIndex]; + patchW[cellIndex] = patchW[cellIndex]*alphaCell[cellIndex]; + } + + const vectorField n1 = Foam::vectorField(nF, vector(1.0, 0.0, 0.0)); + const vectorField n2 = Foam::vectorField(nF, vector(0.0, 1.0, 0.0)); + const vectorField n3 = Foam::vectorField(nF, vector(0.0, 0.0, 1.0)); + + operator == (n1*patchU + n2*patchV + n3*patchW); + + fixedValueFvPatchField<vector>::updateCoeffs(); +} + + +void Foam::IH_Waves_InletVelocityFvPatchVectorField::write(Ostream& os) const +{ + fvPatchField<vector>::write(os); + + os.writeKeyword("waveDictName") << waveDictName_ + << token::END_STATEMENT << nl; + + os.writeKeyword("leftORright") << leftORright_ + << token::END_STATEMENT << nl; + + os.writeKeyword("RealwaterDepth") << RealwaterDepth_ + << token::END_STATEMENT << nl; + + os.writeKeyword("initialDepth") << initialDepth_ + << token::END_STATEMENT << nl; + + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + makePatchTypeField + ( + fvPatchVectorField, + IH_Waves_InletVelocityFvPatchVectorField + ); +} + + +// ************************************************************************* // diff --git a/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/IH_Waves_InletVelocityFvPatchVectorField.H b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/IH_Waves_InletVelocityFvPatchVectorField.H new file mode 100644 index 0000000000000000000000000000000000000000..d52984d738d0bf4c89344b5176b62967b682119d --- /dev/null +++ b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/IH_Waves_InletVelocityFvPatchVectorField.H @@ -0,0 +1,238 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2006-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::IH_Waves_InletVelocityFvPatchVectorField + +Description + Describes a volumetric/mass flow normal vector boundary condition by its + magnitude as an integral over its area. + + The basis of the patch (volumetric or mass) is determined by the + dimensions of the flux, phi. + The current density is used to correct the velocity when applying the + mass basis. + + Example of the boundary condition specification: + @verbatim + inlet + { + type IH_Waves_InletVelocityV2; + waveDict IHWavesDict; + value uniform (0 0 0); + leftORright 1.0; + + } + @endverbatim + + \heading Function object usage + \table + Property | Description | Required | Default value + type | type: IH_Waves_InletVelocityV2 | yes + waveDict | Dictionary where variables for generation/absorption are defined | yes | IHWavesDict + leftORright | Define location of Boundary condition: Left(1) or Right (-1) | yes | -1 + \endtable + +Note + - The value is positive inwards + - May not work correctly for transonic inlets + - Strange behaviour with potentialFoam since the U equation is not solved + +SourceFiles + IH_Waves_InletVelocityFvPatchVectorField.C + +\*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + +#ifndef IH_Waves_InletVelocityFvPatchVectorField_H +#define IH_Waves_InletVelocityFvPatchVectorField_H + +#include "fixedValueFvPatchFields.H" +#include "mathematicalConstants.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +/*---------------------------------------------------------------------------*\ + Class IH_Waves_InletVelocityFvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class IH_Waves_InletVelocityFvPatchVectorField +: + public fixedValueFvPatchVectorField +{ + // Private data + + //- Wave period (seconds) + scalar wavePeriod_; + + //- Wave height (meters) + scalar waveHeight_; + + //- Wave length (meters) + scalar waveLength_; + + //- Water depth (meters) + scalar waterDepth_; + + //- Wave phase (radians) + scalar wavePhase_; + + //- Lambda StokesV parameter + scalar lambdaStokesV_; + + //- m Cnoidal parameter + scalar mCnoidal_; + + //- Generation + Absorption at the same time in the boundary (on=1, off=0) + bool genAbs_; + + //- Number of different paddles (for absorption) + label nPaddles_; + + //- Smoothing factor for the wave(s) (seconds) + scalar tSmooth_; + + //Aborption Left(1) or Right(-1) + scalar leftORright_; + + //- Dictionary name + word waveDictName_; + + //- Regular or Irregular + word waveType_; + + //- Name of the theory + word waveTheory_; + + //- BC check and read values just the first time + bool allCheck_; + + //- Direction of propagation (degrees, from X axis) + scalar waveDir_; + + //- Initial water depth (referece) + scalar RealwaterDepth_; + + //- Water depth (meters) + scalar initialDepth_; + +public: + + //- Runtime type information + TypeName("IH_Waves_InletVelocityV2"); + + + // Constructors + + //- Construct from patch and internal field + IH_Waves_InletVelocityFvPatchVectorField + ( + const fvPatch&, + const DimensionedField<vector, volMesh>& + ); + + //- Construct from patch, internal field and dictionary + IH_Waves_InletVelocityFvPatchVectorField + ( + const fvPatch&, + const DimensionedField<vector, volMesh>&, + const dictionary& + ); + + //- Construct by mapping given + // IH_Waves_InletVelocityFvPatchVectorField + // onto a new patch + IH_Waves_InletVelocityFvPatchVectorField + ( + const IH_Waves_InletVelocityFvPatchVectorField&, + const fvPatch&, + const DimensionedField<vector, volMesh>&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + IH_Waves_InletVelocityFvPatchVectorField + ( + const IH_Waves_InletVelocityFvPatchVectorField& + ); + + //- Construct and return a clone + virtual tmp<fvPatchVectorField> clone() const + { + return tmp<fvPatchVectorField> + ( + new IH_Waves_InletVelocityFvPatchVectorField(*this) + ); + } + + //- Construct as copy setting internal field reference + IH_Waves_InletVelocityFvPatchVectorField + ( + const IH_Waves_InletVelocityFvPatchVectorField&, + const DimensionedField<vector, volMesh>& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp<fvPatchVectorField> clone + ( + const DimensionedField<vector, volMesh>& iF + ) const + { + return tmp<fvPatchVectorField> + ( + new IH_Waves_InletVelocityFvPatchVectorField(*this, iF) + ); + } + + + // Member functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Write + virtual void write(Ostream&) const; + + // Other common member functions + #include "memberFun.H" + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileBoussinesq.H b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileBoussinesq.H new file mode 100644 index 0000000000000000000000000000000000000000..fbbd38b91d44196e855fee21eaa52bf3f8b919f4 --- /dev/null +++ b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileBoussinesq.H @@ -0,0 +1,170 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + + if (zSup[cellIndex] < heights[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = BoussinesqFun :: U + ( + waveHeight_, + RealwaterDepth_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + waveAngle, + currTime, + X0, + patchHeight[cellIndex] + ); + + patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) * timeMult; + patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) * timeMult; + + patchW[cellIndex] = BoussinesqFun :: W + ( + waveHeight_, + RealwaterDepth_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + waveAngle, + currTime, + X0, + patchHeight[cellIndex] + ); + + patchW[cellIndex] = patchW[cellIndex] * timeMult; + + } + else if ( zInf[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = 0.0; + patchV[cellIndex] = 0.0; + patchW[cellIndex] = 0.0; + } + else + { + if ( calculatedLevel[cellGroup[cellIndex]-1] + < measuredLevels[cellGroup[cellIndex]-1] ) + { + if ((zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1])) + { + auxiliar = calculatedLevel[cellGroup[cellIndex]-1] + - zInf[cellIndex]; + auxiliarTotal = zSup[cellIndex]-zInf[cellIndex]; + auxiliarTotal = auxiliar/auxiliarTotal; + + patchU[cellIndex] = BoussinesqFun :: U + ( + waveHeight_, + RealwaterDepth_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + waveAngle, + currTime, + X0, + patchHeight[cellIndex] + ); + + patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] + * sin(waveAngle) * timeMult; + patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] + * cos(waveAngle) * timeMult; + + patchW[cellIndex] = BoussinesqFun :: W + ( + waveHeight_, + RealwaterDepth_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + waveAngle, + currTime, + X0, + patchHeight[cellIndex] + ); + + patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] + * timeMult; + } + } + else if ( calculatedLevel[cellGroup[cellIndex]-1] + > measuredLevels[cellGroup[cellIndex]-1] ) + { + if (zSup[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = BoussinesqFun :: U + ( + waveHeight_, + RealwaterDepth_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + waveAngle, + currTime, + X0, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) + * timeMult; + patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) + * timeMult; + + patchW[cellIndex] = BoussinesqFun :: W + ( + waveHeight_, + RealwaterDepth_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + waveAngle, + currTime, + X0, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchW[cellIndex] = patchW[cellIndex]*timeMult; + } + else if ((zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) ) + { + auxiliar = calculatedLevel[cellGroup[cellIndex]-1] + - zInf[cellIndex]; + auxiliarTotal = zSup[cellIndex]-zInf[cellIndex]; + auxiliarTotal = auxiliar/auxiliarTotal; + + patchU[cellIndex] = BoussinesqFun :: U + ( + waveHeight_, + RealwaterDepth_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + waveAngle, + currTime, + X0, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] + * sin(waveAngle) * timeMult; + patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] + * cos(waveAngle) * timeMult; + + patchW[cellIndex] = BoussinesqFun :: W + ( + waveHeight_, + RealwaterDepth_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + waveAngle, + currTime, + X0, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] + * timeMult; + } + } + } diff --git a/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileCnoidal.H b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileCnoidal.H new file mode 100644 index 0000000000000000000000000000000000000000..132f388d427efa1ce2e11d1316814203e2b8151d --- /dev/null +++ b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileCnoidal.H @@ -0,0 +1,184 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + + if (zSup[cellIndex] < heights[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = cnoidalFun :: U + ( + waveHeight_, + RealwaterDepth_, + mCnoidal_, + waveKx, + waveKy, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + patchHeight[cellIndex] + ); + + patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) * timeMult; + patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) * timeMult; + + patchW[cellIndex] = cnoidalFun :: W + ( + waveHeight_, + RealwaterDepth_, + mCnoidal_, + waveKx, + waveKy, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + patchHeight[cellIndex] + ); + + patchW[cellIndex] = patchW[cellIndex]*timeMult; + + } + else if ( zInf[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = 0.0; + patchV[cellIndex] = 0.0; + patchW[cellIndex] = 0.0; + } + else + { + if ( calculatedLevel[cellGroup[cellIndex]-1] + < measuredLevels[cellGroup[cellIndex]-1] ) + { + if ( (zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1])) + { + auxiliar = calculatedLevel[cellGroup[cellIndex]-1] - zInf[cellIndex]; + auxiliarTotal = zSup[cellIndex]-zInf[cellIndex]; + auxiliarTotal = auxiliar/auxiliarTotal; + + patchU[cellIndex] = cnoidalFun :: U + ( + waveHeight_, + RealwaterDepth_, + mCnoidal_, + waveKx, + waveKy, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + patchHeight[cellIndex] + ); + + patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] + * sin(waveAngle)*timeMult; + patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] + * cos(waveAngle)*timeMult; + + patchW[cellIndex] = cnoidalFun :: W + ( + waveHeight_, + RealwaterDepth_, + mCnoidal_, + waveKx, + waveKy, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + patchHeight[cellIndex] + ); + + patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] + * timeMult; + } + } + else if ( calculatedLevel[cellGroup[cellIndex]-1] + > measuredLevels[cellGroup[cellIndex]-1] ) + { + if (zSup[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = cnoidalFun :: U + ( + waveHeight_, + RealwaterDepth_, + mCnoidal_, + waveKx, + waveKy, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) + * timeMult; + patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) + * timeMult; + + patchW[cellIndex] = cnoidalFun :: W + ( + waveHeight_, + RealwaterDepth_, + mCnoidal_, + waveKx, + waveKy, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchW[cellIndex] = patchW[cellIndex]*timeMult; + } + else if ((zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) ) + { + auxiliar = calculatedLevel[cellGroup[cellIndex]-1] - zInf[cellIndex]; + auxiliarTotal = zSup[cellIndex]-zInf[cellIndex]; + auxiliarTotal = auxiliar/auxiliarTotal; + + patchU[cellIndex] = cnoidalFun :: U + ( + waveHeight_, + RealwaterDepth_, + mCnoidal_, + waveKx, + waveKy, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] + * sin(waveAngle) * timeMult; + patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] + * cos(waveAngle) * timeMult; + + patchW[cellIndex] = cnoidalFun :: W + ( + waveHeight_, + RealwaterDepth_, + mCnoidal_, + waveKx, + waveKy, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] + * timeMult; + } + } + } diff --git a/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileGENAB.H b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileGENAB.H new file mode 100644 index 0000000000000000000000000000000000000000..56ba8c1ce56ca194822058b576570bc1af3a9a70 --- /dev/null +++ b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileGENAB.H @@ -0,0 +1,20 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + + if ( genAbs_ ) + { + if ( zInf[cellIndex] >= measuredLevelsGENAB[cellGroup[cellIndex]-1] ) + { + patchUABS[cellIndex] = 0.0; + } + else + { + patchUABS[cellIndex] = ( calculatedLevel[cellGroup[cellIndex]-1] - measuredLevelsGENAB[cellGroup[cellIndex]-1] ) * sqrt(g/measuredLevelsGENAB[cellGroup[cellIndex]-1]); + } + } + patchU[cellIndex] = patchU[cellIndex] + leftORright_ * patchUABS[cellIndex]; diff --git a/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileStokesI.H b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileStokesI.H new file mode 100644 index 0000000000000000000000000000000000000000..b578c69c86cc1489b13dd702957fc1ecf713d20b --- /dev/null +++ b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileStokesI.H @@ -0,0 +1,187 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + + if (zSup[cellIndex] < heights[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = StokesIFun :: U + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + patchHeight[cellIndex] + ); + + patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) * timeMult; + patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) * timeMult; + + patchW[cellIndex] = StokesIFun :: W + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + patchHeight[cellIndex] + ); + + patchW[cellIndex] = patchW[cellIndex]*timeMult; + } + else if ( zInf[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = 0.0; + patchV[cellIndex] = 0.0; + patchW[cellIndex] = 0.0; + } + else + { + if ( calculatedLevel[cellGroup[cellIndex]-1] + < measuredLevels[cellGroup[cellIndex]-1] ) + { + if ( (zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) ) + { + auxiliar = calculatedLevel[cellGroup[cellIndex]-1] + - zInf[cellIndex]; + auxiliarTotal = zSup[cellIndex]-zInf[cellIndex]; + auxiliarTotal = auxiliar/auxiliarTotal; + + patchU[cellIndex] = StokesIFun :: U + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + patchHeight[cellIndex] + ); + + patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] + * sin(waveAngle)*timeMult; + patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] + * cos(waveAngle)*timeMult; + + patchW[cellIndex] = StokesIFun :: W + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + patchHeight[cellIndex] + ); + + patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] + * timeMult; + + } + } + else if ( calculatedLevel[cellGroup[cellIndex]-1] + > measuredLevels[cellGroup[cellIndex]-1] ) + { + if (zSup[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = StokesIFun :: U + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchV[cellIndex] = patchU[cellIndex] + * sin(waveAngle)*timeMult; + patchU[cellIndex] = patchU[cellIndex] + * cos(waveAngle)*timeMult; + + patchW[cellIndex] = StokesIFun :: W + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchW[cellIndex] = patchW[cellIndex]*timeMult; + + } + else if ((zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1])) + { + auxiliar = calculatedLevel[cellGroup[cellIndex]-1] + - zInf[cellIndex]; + auxiliarTotal = zSup[cellIndex]-zInf[cellIndex]; + auxiliarTotal = auxiliar/auxiliarTotal; + + patchU[cellIndex] = StokesIFun :: U + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] + * sin(waveAngle) * timeMult; + patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] + * cos(waveAngle) * timeMult; + + patchW[cellIndex] = StokesIFun :: W + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] + * timeMult; + } + } + } diff --git a/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileStokesII.H b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileStokesII.H new file mode 100644 index 0000000000000000000000000000000000000000..3313587706b0409c8bf6406184f4c5f57258ef61 --- /dev/null +++ b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileStokesII.H @@ -0,0 +1,186 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + + if (zSup[cellIndex] < heights[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = StokesIIFun :: U + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + patchHeight[cellIndex] + ); + + patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) * timeMult; + patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) * timeMult; + + patchW[cellIndex] = StokesIIFun :: W + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + patchHeight[cellIndex] + ); + + patchW[cellIndex] = patchW[cellIndex]*timeMult; + + } + else if ( zInf[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = 0.0; + patchV[cellIndex] = 0.0; + patchW[cellIndex] = 0.0; + } + else + { + if ( calculatedLevel[cellGroup[cellIndex]-1] + < measuredLevels[cellGroup[cellIndex]-1] ) + { + if ( (zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) ) + { + auxiliar = calculatedLevel[cellGroup[cellIndex]-1] + - zInf[cellIndex]; + auxiliarTotal = zSup[cellIndex]-zInf[cellIndex]; + auxiliarTotal = auxiliar/auxiliarTotal; + + patchU[cellIndex] = StokesIIFun :: U + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + patchHeight[cellIndex] + ); + + patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] + * sin(waveAngle)*timeMult; + patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] + * cos(waveAngle)*timeMult; + + patchW[cellIndex] = StokesIIFun :: W + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + patchHeight[cellIndex] + ); + + patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] + * timeMult; + } + } + else if ( calculatedLevel[cellGroup[cellIndex]-1] + > measuredLevels[cellGroup[cellIndex]-1] ) + { + if (zSup[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = StokesIIFun :: U + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchV[cellIndex] = patchU[cellIndex] + * sin(waveAngle)*timeMult; + patchU[cellIndex] = patchU[cellIndex] + * cos(waveAngle)*timeMult; + + patchW[cellIndex] = StokesIIFun :: W + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchW[cellIndex] = patchW[cellIndex]*timeMult; + } + else if ((zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1])) + { + auxiliar = calculatedLevel[cellGroup[cellIndex]-1] + - zInf[cellIndex]; + auxiliarTotal = zSup[cellIndex]-zInf[cellIndex]; + auxiliarTotal = auxiliar/auxiliarTotal; + + patchU[cellIndex] = StokesIIFun :: U + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] + * sin(waveAngle) * timeMult; + patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] + * cos(waveAngle) * timeMult; + + patchW[cellIndex] = StokesIIFun :: W + ( + waveHeight_, + RealwaterDepth_, + waveKx, + xGroup[cellGroup[cellIndex]-1], + waveKy, + yGroup[cellGroup[cellIndex]-1], + waveOmega, + currTime, + wavePhase_, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] + * timeMult; + } + } + } diff --git a/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileStokesV.H b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileStokesV.H new file mode 100644 index 0000000000000000000000000000000000000000..00e593aa12b61374183b0cc3583d794e1006ebd8 --- /dev/null +++ b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/profileStokesV.H @@ -0,0 +1,186 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + + if (zSup[cellIndex] < heights[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = stokesVFun :: U + ( + RealwaterDepth_, + waveKx, + waveKy, + lambdaStokesV_, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + wavePhase_, + patchHeight[cellIndex] + ); + + patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) * timeMult; + patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) * timeMult; + + patchW[cellIndex] = stokesVFun :: W + ( + RealwaterDepth_, + waveKx, + waveKy, + lambdaStokesV_, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + wavePhase_, + patchHeight[cellIndex] + ); + + patchW[cellIndex] = patchW[cellIndex]*timeMult; + + } + else if ( zInf[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = 0.0; + patchV[cellIndex] = 0.0; + patchW[cellIndex] = 0.0; + } + else + { + if ( calculatedLevel[cellGroup[cellIndex]-1] < + measuredLevels[cellGroup[cellIndex]-1] ) + { + if ( (zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) ) + { + auxiliar = calculatedLevel[cellGroup[cellIndex]-1] + - zInf[cellIndex]; + auxiliarTotal = zSup[cellIndex]-zInf[cellIndex]; + auxiliarTotal = auxiliar/auxiliarTotal; + + patchU[cellIndex] = stokesVFun :: U + ( + RealwaterDepth_, + waveKx, + waveKy, + lambdaStokesV_, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + wavePhase_, + patchHeight[cellIndex] + ); + + patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] + * sin(waveAngle)*timeMult; + patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] + * cos(waveAngle)*timeMult; + + patchW[cellIndex] = stokesVFun :: W + ( + RealwaterDepth_, + waveKx, + waveKy, + lambdaStokesV_, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + wavePhase_, + patchHeight[cellIndex] + ); + + patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] + * timeMult; + } + } + else if ( calculatedLevel[cellGroup[cellIndex]-1] + > measuredLevels[cellGroup[cellIndex]-1] ) + { + if (zSup[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) + { + patchU[cellIndex] = stokesVFun :: U + ( + RealwaterDepth_, + waveKx, + waveKy, + lambdaStokesV_, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + wavePhase_, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchV[cellIndex] = patchU[cellIndex] * sin(waveAngle) + * timeMult; + patchU[cellIndex] = patchU[cellIndex] * cos(waveAngle) + * timeMult; + + patchW[cellIndex] = stokesVFun :: W + ( + RealwaterDepth_, + waveKx, + waveKy, + lambdaStokesV_, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + wavePhase_, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchW[cellIndex] = patchW[cellIndex]*timeMult; + } + else if ( (zSup[cellIndex] > calculatedLevel[cellGroup[cellIndex]-1]) + & (zInf[cellIndex] < calculatedLevel[cellGroup[cellIndex]-1]) ) + { + auxiliar = calculatedLevel[cellGroup[cellIndex]-1] + - zInf[cellIndex]; + auxiliarTotal = zSup[cellIndex]-zInf[cellIndex]; + auxiliarTotal = auxiliar/auxiliarTotal; + + patchU[cellIndex] = stokesVFun :: U + ( + RealwaterDepth_, + waveKx, + waveKy, + lambdaStokesV_, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + wavePhase_, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchV[cellIndex] = auxiliarTotal * patchU[cellIndex] + * sin(waveAngle) * timeMult; + patchU[cellIndex] = auxiliarTotal * patchU[cellIndex] + * cos(waveAngle) * timeMult; + + patchW[cellIndex] = stokesVFun :: W + ( + RealwaterDepth_, + waveKx, + waveKy, + lambdaStokesV_, + wavePeriod_, + xGroup[cellGroup[cellIndex]-1], + yGroup[cellGroup[cellIndex]-1], + currTime, + wavePhase_, + measuredLevels[cellGroup[cellIndex]-1] + ); + + patchW[cellIndex] = auxiliarTotal * patchW[cellIndex] + * timeMult; + } + } + } diff --git a/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/velocityProfile.H b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/velocityProfile.H new file mode 100644 index 0000000000000000000000000000000000000000..1473aa8eb15cb2c58724872399358eaa60b0146a --- /dev/null +++ b/integration/genAbs/waveGeneration/IH_Waves_InletVelocity/velProfiles/velocityProfile.H @@ -0,0 +1,36 @@ +/*---------------------------------------------------------------------------*\ + IH-Cantabria 2015 (http://www.ihcantabria.com/en/) + IHFOAM 2015 (http://ihfoam.ihcantabria.com/) + + Author(s): Javier Lopez Lara (jav.lopez@unican.es) + Gabriel Barajas (barajasg@unican.es) +\*---------------------------------------------------------------------------*/ + +if ( waveType_ == "regular" ) +{ + if ( waveTheory_ == "StokesI" ) + { + #include "profileStokesI.H" + #include "profileGENAB.H" + } + else if ( waveTheory_ == "StokesII" ) + { + #include "profileStokesII.H" + #include "profileGENAB.H" + } + else if ( waveTheory_ == "StokesV" ) + { + #include "profileStokesV.H" + #include "profileGENAB.H" + } + else if ( waveTheory_ == "Cnoidal" ) + { + #include "profileCnoidal.H" + #include "profileGENAB.H" + } +} +else if ( waveType_ == "solitary" ) +{ + #include "profileBoussinesq.H" + #include "profileGENAB.H" +} diff --git a/integration/genAbs/waveGeneration/Make/files b/integration/genAbs/waveGeneration/Make/files new file mode 100644 index 0000000000000000000000000000000000000000..0d13760da933554a99a67e5ad59736902f23c6eb --- /dev/null +++ b/integration/genAbs/waveGeneration/Make/files @@ -0,0 +1,5 @@ +IH_Waves_InletAlpha/IH_Waves_InletAlphaFvPatchScalarField.C +IH_Waves_InletVelocity/IH_Waves_InletVelocityFvPatchVectorField.C +../common/waveFun.C + +LIB = $(FOAM_USER_LIBBIN)/libIHwaveGenerationV2esi diff --git a/integration/genAbs/waveGeneration/Make/options b/integration/genAbs/waveGeneration/Make/options new file mode 100644 index 0000000000000000000000000000000000000000..d09270c43b7fdb670c50a2ad82489599a351082c --- /dev/null +++ b/integration/genAbs/waveGeneration/Make/options @@ -0,0 +1,16 @@ +sinclude $(GENERAL_RULES)/mplib$(WM_MPLIB) +sinclude $(RULES)/mplib$(WM_MPLIB) + +EXE_INC = \ + -DOFVERSION=$(OF_VERSION) \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I./IH_Waves_InletVelocity/velProfiles \ + -I./IH_Waves_InletVelocityWind/velProfilesWind \ + -I../common \ + -I../common/checks \ + -I../common/calculateWaterLevel \ + $(PFLAGS) $(PINC) + +LIB_LIBS = \ + -lfiniteVolume \ + $(PLIBS) diff --git a/integration/genAbs/waveGeneration/localMake b/integration/genAbs/waveGeneration/localMake new file mode 100755 index 0000000000000000000000000000000000000000..1b24c72f88d481b5bb646ad960ec8c994fd6be71 --- /dev/null +++ b/integration/genAbs/waveGeneration/localMake @@ -0,0 +1,25 @@ +#!/bin/bash + +if [ $WM_PROJECT == "foam" ]; then + ofversion=`echo $WM_PROJECT_VERSION | sed -e 's/\.x/-9/' -e 's/\./\'$'\n/g' -e 's/-/\'$'\n/g' | grep "[0-9]" | head -2 | tr -d '\n'` +else + ofversion=`echo $WM_PROJECT_VERSION"-0" | sed -e 's/\.x/-9/' -e 's/\./\'$'\n/g' -e 's/-/\'$'\n/g' -e 's/+/\'$'\n/g' -e 's/v/\'$'\n/g' | grep "[0-9]" | head -3 | tr -d '\n'` +fi + +#IHC_dbg +echo $ofversion +#---- + +export OF_VERSION=$ofversion + +wclean + +wmake libso + +if (( $? )) ; then + echo "\n\nIH Wave generation boundary conditions (V2.0) compilation failed" + exit; else + echo -e "\n\nIH Wave generation boundary conditions (V2.0) compiled successfully for $WM_PROJECT $WM_PROJECT_VERSION\n\n\n"; +fi + +wclean