SimFlowCFD / RapidCFD-dev

RapidCFD is an OpenFOAM fork running fully on CUDA platform. Brought to you by
https://sim-flow.com
Other
325 stars 94 forks source link

supersonicFreestream #102

Closed chaous closed 1 year ago

chaous commented 1 year ago

For some reason I'm not able to use supersonicFreestream boundary.

--> FOAM FATAL IO ERROR: Unknown patchField type supersonicFreestream for patch type patch

Valid patchField types are :

63 ( SRFFreestreamVelocity SRFVelocity advective calculated codedFixedValue codedMixed compressible::kqRWallFunction cyclic cyclicACMI cyclicAMI cyclicSlip cylindricalInletVelocity directionMixed empty externalCoupled fixedGradient fixedInternalValue fixedJump fixedJumpAMI fixedMean fixedNormalInletOutletVelocity fixedNormalSlip fixedValue flowRateInletVelocity fluxCorrectedVelocity freestream inletOutlet interstitialInletVelocity mixed movingWallVelocity nonuniformTransformCyclic oscillatingFixedValue outletInlet outletPhaseMeanVelocity partialSlip pressureDirectedInletOutletVelocity pressureDirectedInletVelocity pressureInletOutletParSlipVelocity pressureInletOutletVelocity pressureInletUniformVelocity pressureInletVelocity pressureNormalInletOutletVelocity processor processorCyclic rotatingPressureInletOutletVelocity rotatingWallVelocity sliced slip surfaceNormalFixedValue swirlFlowRateInletVelocity symmetry symmetryPlane translatingWallVelocity turbulentInlet uniformFixedGradient uniformFixedValue uniformInletOutlet uniformJump uniformJumpAMI variableHeightFlowRateInletVelocity waveTransmissive wedge zeroGradient )

TonkomoLLC commented 1 year ago

For reasons that I do not know (maybe it does not work?), the patch you refer to is commented out here.

You can try removing the comments and see what happens...

chaous commented 1 year ago

For some reason Rapid-cfd doesn't see this bc. I replaced another bc with code for supersonicFreestream (waiting for results)

chaous commented 1 year ago

(I'm using docker from this issue https://github.com/Atizar/RapidCFD-dev/issues/97)

TonkomoLLC commented 1 year ago

Just to make sure there is no misunderstanding, you will need to remove the comments in the referenced file and then re-compile that library in the src/finiteVolume folder with wmake

chaous commented 1 year ago

ok, I will try to do it

chaous commented 1 year ago

I uncommented supersonicFreestream in Make file and got this fields/fvPatchFields/derived/supersonicFreestream/supersonicFreestreamFvPatchVectorField.C(180): error: no operator "[]" matches these operands operand types are: const Foam::fvPatchField [ int ]

chaous commented 1 year ago

I will try to add [] to fvPatchField class

TonkomoLLC commented 1 year ago

I am not sure what version of CUDA you are using or if your GPU is "old" enough to allow it, but I would consider using CUDA 11.0 or 11.1 when debugging this boundary condition. As reported in #92, there are some occasional issues with CUDA versions starting at 11.2.

I looked at line 180 in supersonicFreestreamFvPatchVectorField.C and didn't see anything obvious that would cause the compilation error that you received.

Good luck with your debugging!

chaous commented 1 year ago

Issue was that operator [] was not implemented. I implemented it and now got error with another function and truing to fix it. If it will work I will post my code with nacaAerfol example case.

chaous commented 1 year ago

I replaced valueFraction()[facei] = 0; with valueFraction().set(facei, 0); in str/finiteVolume/fields/fvPatchFields/derived/supersonicFreestream/supersonicFreestreamFvPatchVectorField.C

now RapidCFD is compiling properly however I'm getting Courant Number mean: 0.365855 max: 2.02613e+09 diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0 AINVPBiCG: Solving for Ux, Initial residual = 0.948318, Final residual = 0.927013, No Iterations 1001 AINVPBiCG: Solving for Uy, Initial residual = 0.947164, Final residual = 0.934587, No Iterations 1001 AINVPBiCG: Solving for e, Initial residual = 1, Final residual = 1.8177e+17, No Iterations 1001 AINVPBiCG: Solving for p, Initial residual = -nan, Final residual = -nan, No Iterations 1001 diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0 time step continuity errors : sum local = -nan, global = -nan, cumulative = -nan

--> FOAM FATAL ERROR: unphysical subsonic inflow has been generated on patch INLE1 of field U in file "/home/ofuser/nacaAirfoil/0/U"

From function supersonicFreestreamFvPatchVectorField::updateCoeffs()
in file fields/fvPatchFields/derived/supersonicFreestream/supersonicFreestreamFvPatchVectorField.C at line 293.

Does someone know how to hix it? (I'm ruining nacaAerfol example https://github.com/OpenFOAM/OpenFOAM-2.3.x/tree/master/tutorials/compressible/sonicFoam/ras/nacaAirfoil)

TonkomoLLC commented 1 year ago

Hi, would you be willing to place your modified source code for the supersonicFreeStream patch here? I do not know the answer to your question right now. What version of CUDA are you using? thanks -

chaous commented 1 year ago

I'm using CUDA 11.5 (in docker from that https://github.com/Atizar/RapidCFD-dev/issues/97 issue)

chaous commented 1 year ago

I added inline const T& operator[](const label n) const { return *constiterator(v+n); }

to https://github.com/Atizar/RapidCFD-dev/blob/0fc3af285e3551234ffff9f9402c4c42f088af66/src/OpenFOAM/containers/Lists/gpuList/gpuList.H

chaous commented 1 year ago

/---------------------------------------------------------------------------\ ========= | \ / F ield | OpenFOAM: The Open Source CFD Toolbox \ / O peration | \ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \/ 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/>.

*---------------------------------------------------------------------------*/

include "supersonicFreestreamFvPatchVectorField.H"

include "addToRunTimeSelectionTable.H"

include "fvPatchFieldMapper.H"

include "volFields.H"

// Constructors //

Foam::supersonicFreestreamFvPatchVectorField:: supersonicFreestreamFvPatchVectorField ( const fvPatch& p, const DimensionedField<vector, volMesh>& iF ) : mixedFvPatchVectorField(p, iF), TName("T"), pName("p"), psiName("thermo:psi"), UInf(vector::zero), pInf(0), TInf(0), gamma_(0) { refValue() = patchInternalField(); refGrad() = vector::zero; valueFraction() = 1; }

Foam::supersonicFreestreamFvPatchVectorField:: supersonicFreestreamFvPatchVectorField ( const supersonicFreestreamFvPatchVectorField& ptf, const fvPatch& p, const DimensionedField<vector, volMesh>& iF, const fvPatchFieldMapper& mapper ) : mixedFvPatchVectorField(ptf, p, iF, mapper), TName(ptf.TName), pName(ptf.pName), psiName(ptf.psiName), UInf(ptf.UInf), pInf(ptf.pInf), TInf(ptf.TInf), gamma(ptf.gamma) {}

Foam::supersonicFreestreamFvPatchVectorField:: supersonicFreestreamFvPatchVectorField ( const fvPatch& p, const DimensionedField<vector, volMesh>& iF, const dictionary& dict ) : mixedFvPatchVectorField(p, iF), TName(dict.lookupOrDefault("T", "T")), pName(dict.lookupOrDefault("p", "p")), psiName(dict.lookupOrDefault("psi", "thermo:psi")), UInf(dict.lookup("UInf")), pInf(readScalar(dict.lookup("pInf"))), TInf(readScalar(dict.lookup("TInf"))), gamma_(readScalar(dict.lookup("gamma"))) { if (dict.found("value")) { fvPatchField::operator= ( vectorField("value", dict, p.size()) ); } else { fvPatchField::operator=(patchInternalField()); }

refValue() = *this;
refGrad() = vector::zero;
valueFraction() = 1;

if (pInf_ < SMALL)
{
    FatalIOErrorIn
    (
        "supersonicFreestreamFvPatchVectorField::"
        "supersonicFreestreamFvPatchVectorField"
        "(const fvPatch&, const vectorField&, const dictionary&)",
        dict
    )   << "    unphysical pInf specified (pInf <= 0.0)"
        << "\n    on patch " << this->patch().name()
        << " of field " << this->dimensionedInternalField().name()
        << " in file " << this->dimensionedInternalField().objectPath()
        << exit(FatalIOError);
}

}

Foam::supersonicFreestreamFvPatchVectorField:: supersonicFreestreamFvPatchVectorField ( const supersonicFreestreamFvPatchVectorField& sfspvf ) : mixedFvPatchVectorField(sfspvf), TName(sfspvf.TName), pName(sfspvf.pName), psiName(sfspvf.psiName), UInf(sfspvf.UInf), pInf(sfspvf.pInf), TInf(sfspvf.TInf), gamma(sfspvf.gamma) {}

Foam::supersonicFreestreamFvPatchVectorField:: supersonicFreestreamFvPatchVectorField ( const supersonicFreestreamFvPatchVectorField& sfspvf, const DimensionedField<vector, volMesh>& iF ) : mixedFvPatchVectorField(sfspvf, iF), TName(sfspvf.TName), pName(sfspvf.pName), psiName(sfspvf.psiName), UInf(sfspvf.UInf), pInf(sfspvf.pInf), TInf(sfspvf.TInf), gamma(sfspvf.gamma) {}

// Member Functions //

void Foam::supersonicFreestreamFvPatchVectorField::updateCoeffs() { if (!size() || updated()) { return; }

const fvPatchField<scalar>& pT =
    patch().lookupPatchField<volScalarField, scalar>(TName_);

const fvPatchField<scalar>& pp =
    patch().lookupPatchField<volScalarField, scalar>(pName_);

const fvPatchField<scalar>& ppsi =
    patch().lookupPatchField<volScalarField, scalar>(psiName_);

// Need R of the free-stream flow.  Assume R is independent of location
// along patch so use face 0
scalar R = 1.0/(ppsi[0]*pT[0]);

scalar MachInf = mag(UInf_)/sqrt(gamma_*R*TInf_);

if (MachInf < 1.0)
{
    FatalErrorIn
    (
        "supersonicFreestreamFvPatchVectorField::updateCoeffs()"
    )   << "    MachInf < 1.0, free stream must be supersonic"
        << "\n    on patch " << this->patch().name()
        << " of field " << this->dimensionedInternalField().name()
        << " in file " << this->dimensionedInternalField().objectPath()
        << exit(FatalError);
}

//vectorField& Up = refValue();
Foam::gpuField<Foam::vector>& Up = refValue();
valueFraction() = 1;

// get the near patch internal cell values
const vectorField U(patchInternalField());

// Find the component of U normal to the free-stream flow and in the
// plane of the free-stream flow and the patch normal

// Direction of the free-stream flow
vector UInfHat = UInf_/mag(UInf_);

// Normal to the plane defined by the free-stream and the patch normal
//tmp<vectorField> nnInfHat = UInfHat ^ patch().nf();
auto nnInfHat = UInfHat ^ patch().nf();

// Normal to the free-stream in the plane defined by the free-stream
// and the patch normal
const vectorField nHatInf(nnInfHat ^ UInfHat);

// Component of U normal to the free-stream in the plane defined by the
// free-stream and the patch normal
const vectorField Un(nHatInf*(nHatInf & U));

// The tangential component is
const vectorField Ut(U - Un);

// Calculate the Prandtl-Meyer function of the free-stream
scalar nuMachInf =
    sqrt((gamma_ + 1)/(gamma_ - 1))
   *atan(sqrt((gamma_ - 1)/(gamma_ + 1)*(sqr(MachInf) - 1)))
  - atan(sqr(MachInf) - 1);

// Set the patch boundary condition based on the Mach number and direction
// of the flow dictated by the boundary/free-stream pressure difference

forAll(Up, facei)
{
    if (pp[facei] >= pInf_) // If outflow
    {
        // Assume supersonic outflow and calculate the boundary velocity
        // according to ???

        scalar fpp =
            sqrt(sqr(MachInf) - 1)
           /(gamma_*sqr(MachInf))*mag(Ut[facei])*log(pp[facei]/pInf_);

        //Up[facei] = Ut[facei] + fpp*nHatInf[facei];
        Up.set(facei, Ut[facei] + fpp*nHatInf[facei]);

        // Calculate the Mach number of the boundary velocity
        scalar Mach = mag(Up.get(facei))/sqrt(gamma_/ppsi[facei]);

        if (Mach <= 1) // If subsonic
        {
            // Zero-gradient subsonic outflow

            Up.set(facei, U[facei]);
            //valueFraction()[facei] = 0;
            valueFraction().set(facei, 0);
        }
    }
    else // if inflow
    {
        // Calculate the Mach number of the boundary velocity
        // from the boundary pressure assuming constant total pressure
        // expansion from the free-stream
        scalar Mach =
            sqrt
            (
                (2/(gamma_ - 1))*(1 + ((gamma_ - 1)/2)*sqr(MachInf))
               *pow(pp[facei]/pInf_, (1 - gamma_)/gamma_)
              - 2/(gamma_ - 1)
            );

        if (Mach > 1) // If supersonic
        {
            // Supersonic inflow is assumed to occur according to the
            // Prandtl-Meyer expansion process

            scalar nuMachf =
                sqrt((gamma_ + 1)/(gamma_ - 1))
               *atan(sqrt((gamma_ - 1)/(gamma_ + 1)*(sqr(Mach) - 1)))
              - atan(sqr(Mach) - 1);

            scalar fpp = (nuMachInf - nuMachf)*mag(Ut[facei]);

            Up.set(facei, Ut[facei] + fpp*nHatInf[facei]);
        }
        else // If subsonic
        {
            FatalErrorIn
            (
                "supersonicFreestreamFvPatchVectorField::updateCoeffs()"
            )   << "unphysical subsonic inflow has been generated"
                << "\n    on patch " << this->patch().name()
                << " of field " << this->dimensionedInternalField().name()
                << " in file "
                << this->dimensionedInternalField().objectPath()
                << exit(FatalError);
        }
    }
}

mixedFvPatchVectorField::updateCoeffs();

}

void Foam::supersonicFreestreamFvPatchVectorField::write(Ostream& os) const { fvPatchVectorField::write(os); writeEntryIfDifferent(os, "T", "T", TName); writeEntryIfDifferent(os, "p", "p", pName); writeEntryIfDifferent(os, "psi", "thermo:psi", psiName); os.writeKeyword("UInf") << UInf << token::ENDSTATEMENT << nl; os.writeKeyword("pInf") << pInf << token::ENDSTATEMENT << nl; os.writeKeyword("TInf") << TInf << token::ENDSTATEMENT << nl; os.writeKeyword("gamma") << gamma << token::END_STATEMENT << nl; writeEntry("value", os); }

// * //

namespace Foam { makePatchTypeField ( fvPatchVectorField, supersonicFreestreamFvPatchVectorField ); }

// ***** //

chaous commented 1 year ago

However now I getting ourant Number mean: 0.000174184 max: 0.948172 diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0

0 Foam::error::printStack(Foam::Ostream&) at ??:?

1 Foam::sigSegv::sigHandler(int) at ??:?

2 in "/lib/x86_64-linux-gnu/libc.so.6"

3 Foam::supersonicFreestreamFvPatchVectorField::updateCoeffs() at ??:?

4 Foam::fvMatrix<Foam::Vector >::fvMatrix(Foam::GeometricField<Foam::Vector, Foam::fvPatchField, Foam::volMesh> const&, Foam::dimensionSet const&) in "/opt/OpenFOAM/RapidCFD-dev/platforms/linux64NvccDPOpt/bin/sonicFoam"

5 in "/opt/OpenFOAM/RapidCFD-dev/platforms/linux64NvccDPOpt/bin/sonicFoam"

6 in "/opt/OpenFOAM/RapidCFD-dev/platforms/linux64NvccDPOpt/bin/sonicFoam"

7 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"

8 in "/opt/OpenFOAM/RapidCFD-dev/platforms/linux64NvccDPOpt/bin/sonicFoam"

Segmentation fault (core dumped)

TonkomoLLC commented 1 year ago

Hi, Presently I do not know the answer to your question.

I think you can drag text files into the text box. That will be easier for me to implement your changes than copying text from your reply. Alternatively you can fork RapidCFD from GitHub and make your changes on your own fork. Then I can git clone your fork and compile it.

chaous commented 1 year ago

I forked repository and added changes https://github.com/chaous/RapidCFD-sonic/blob/master/src/finiteVolume/fields/fvPatchFields/derived/supersonicFreestream/supersonicFreestreamFvPatchVectorField.C

chaous commented 1 year ago

Courant Number mean: 0.000174184 max: 0.948172 diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0

0 Foam::error::printStack(Foam::Ostream&) at ??:?

1 Foam::sigSegv::sigHandler(int) at ??:?

2 in "/lib/x86_64-linux-gnu/libc.so.6"

3 Foam::supersonicFreestreamFvPatchVectorField::updateCoeffs() at ??:?

4 Foam::fvMatrix<Foam::Vector >::fvMatrix(Foam::GeometricField<Foam::Vector, Foam::fvPatchField, Foam::volMesh> const&, Foam::dimensionSet const&) in "/opt/OpenFOAM/RapidCFD-dev/platforms/linux64NvccDPOpt/bin/sonicFoam"

5 in "/opt/OpenFOAM/RapidCFD-dev/platforms/linux64NvccDPOpt/bin/sonicFoam"

6 in "/opt/OpenFOAM/RapidCFD-dev/platforms/linux64NvccDPOpt/bin/sonicFoam"

7 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"

8 in "/opt/OpenFOAM/RapidCFD-dev/platforms/linux64NvccDPOpt/bin/sonicFoam"

Segmentation fault (core dumped)

chaous commented 1 year ago

Time = 6e-07

Courant Number mean: 0.000173943 max: 1.55758 diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0 AINVPBiCG: Solving for Ux, Initial residual = 0.901529, Final residual = 0.653648, No Iterations 1001 AINVPBiCG: Solving for Uy, Initial residual = 0.984108, Final residual = 0.792287, No Iterations 1001 AINVPBiCG: Solving for e, Initial residual = 0.99692, Final residual = 0.961589, No Iterations 1001 AINVPBiCG: Solving for p, Initial residual = -nan, Final residual = -nan, No Iterations 1001 diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0 time step continuity errors : sum local = -nan, global = -nan, cumulative = -nan

--> FOAM FATAL ERROR: unphysical subsonic inflow has been generated on patch INLE1 of field U in file "/home/ofuser/nacaAirfoil/0/U"

From function supersonicFreestreamFvPatchVectorField::updateCoeffs()
in file fields/fvPatchFields/derived/supersonicFreestream/supersonicFreestreamFvPatchVectorField.C at line 301.

FOAM exiting

now getting this error again (replaced all [] with get and set comand)

https://github.com/chaous/RapidCFD-sonic/blob/master/src/finiteVolume/fields/fvPatchFields/derived/supersonicFreestream/supersonicFreestreamFvPatchVectorField.C

chaous commented 1 year ago

I looks like I fixed it by making dt = 4e-09 instead of 4e-08

TonkomoLLC commented 1 year ago

Thank you for making the fork. It was easier for me to see what you changed. I see that the changes you made mainly appear to handle the reference to field face IDs compatible with the compiler. I think your point was that these changes were innocuous. If so, I agree. Also, the time step you need to use is very small. I wonder if this eats into the GPU advantage... ?

Upon review of the updated supersonicFreestream patch, I wonder if calculations are setup to utilize the GPU? I refer to, as an example, turbulentInletFvPatchField. You will notice many thrust operations in this boundary condition.

What do you think?

chaous commented 1 year ago

anim 0000 anim 0001 anim 0002 I think my results are not very physical. I will try to understand why because the issue might be in the solver too. Also, I will look at your example too. (I will try to rewrite supersonicFreestream, so it will utilize the GPU) Thank you for your help

chaous commented 1 year ago

(It's velocity visualization)

chaous commented 1 year ago

@TonkomoLLC

TonkomoLLC commented 1 year ago

I agree your results do not look physical ; good luck with the troubleshooting and the check of the GPU calculation in the supersonicFreestream BC.

chaous commented 1 year ago

anim 0001 I fixed my results. Now it works with only dt = 3e-08 instead of 4e-8 as in original example