可能我描述的问题不够准确,下面附上我的部分代码:
在.H文件中,定义了一些私有变量与私有变量的接口:
// private perameters
//- freestream Ma
scalar Ma_;
//- freestream T
scalar Tinf_;
scalar gamma_;
scalar Prt_;
//- munber of modes
scalar nmodes_;
//- smallest wavenumber
scalar dxmin_;
//- ratio of ke and kmin (in wavenumber)
scalar wew1fct_;
//-
scalar amp_;
//- kinetic viscosity of the smallest eddy
scalar visc_;
//- delta Inlet
scalar deltaInlet_;
//- freestream Velocity
scalar Uinf_;
//- freestream pressure
scalar pInf_;
//- read file
IFstream IF_;
//- rows of file data
label ny_;
//- MeanField
vectorField Ubar_;
//- primeField
vectorField UPrime_;
//- ReynoldStress tensor
symmTensorField ReynoldStress_;
label curTimeIndex_;
// member function
void calculateProfiles(IFstream&, vectorField&, symmTensorField&, scalar&);
// Access
scalar nmodes() const
{
return nmodes_;
}
scalar deltaInlet() const
{
return deltaInlet_;
}
scalar& deltaInlet()
{
return deltaInlet_;
}
scalar amp() const
{
return amp_;
}
scalar dxmin() const
{
return dxmin_;
}
scalar visc() const
{
return visc_;
}
scalar wew1fct() const
{
return wew1fct_;
}
scalar Ma() const
{
return Ma_;
}
scalar Tinf() const
{
return Tinf_;
}
scalar gamma() const
{
return gamma_;
}
scalar Prt() const
{
return Prt_;
}
scalar Uinf() const
{
return Uinf_;
}
scalar pInf() const
{
return pInf_;
}
IFstream IF() const
{
return IF_;
}
IFstream& IF()
{
return IF_;
}
vectorField& Ubar()
{
return Ubar_;
}
vectorField Ubar() const
{
return Ubar_;
}
symmTensorField& ReynoldStress()
{
return ReynoldStress_;
}
symmTensorField ReynoldStress() const
{
return ReynoldStress_;
}
label ny()
{
return ny_;
}
然后在.C中,我的构造函数:
// constrctor
Foam::syntheticVelocityFixedValueFvPatchField::syntheticVelocityFixedValueFvPatchField
(
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF,
const dictionary& dict
)
:
fixedValueFvPatchVectorField(p, iF),
Ma_(readScalar(dict.lookup("Ma"))),
Tinf_(readScalar(dict.lookup("Tinf"))),
gamma_(dict.lookupOrDefault("gamma", 1.4)),
Prt_(dict.lookupOrDefault("Prt", 0.9)),
nmodes_(readScalar(dict.lookup("nmodes"))),
dxmin_(readScalar(dict.lookup("dxmin"))),
wew1fct_(dict.lookupOrDefault("wew1fct", 2.0)),
amp_(dict.lookupOrDefault("amp", 1.452762113)),
visc_(dict.lookupOrDefault("visc", 15.2e-6)),
Uinf_(readScalar(dict.lookup("Uinf"))),
pInf_(readScalar(dict.lookup("pInf"))),
IF_(dict.lookup("pathName")), // stored meanStreamVelocity and ReynoldStress profile data
ny_(dict.lookupOrDefault("ny", 535)),
curTimeIndex_(-1)
{
if (dict.found("Ubar") && dict.found("ReynoldStress") && dict.found("deltaInlet"))
{
deltaInlet_= scalar(readScalar(dict.lookup("deltaInlet")));
Ubar_ = vectorField("Ubar", dict, p.size());
ReynoldStress_ = symmTensorField("ReynoldStress", dict, p.size());
}
else
{
calculateProfiles(IF(), Ubar(), ReynoldStress(), deltaInlet());
}
if (dict.found("UPrime"))
{
UPrime_ = vectorField("UPrime", dict, p.size());
}
else
{
UPrime_ = vectorField(p.size(), vector::zero);
}
if (dict.found("value"))
{
fvPatchField<vector>::operator=
(
vectorField("value", dict, p.size())
);
}
else
{
fvPatchField<vector>::operator=(Ubar());
}
}
这里调用的成员函数calculateProperties如下:
// member function
void Foam::syntheticVelocityFixedValueFvPatchField::calculateProfiles
(
IFstream& IF,
vectorField& Ubar,
symmTensorField& ReynoldStress,
scalar& deltaInlet
)
{
if (!IF)
{
FatalErrorIn("void Foam::syntheticVelocityFixedValueFvPatchField::calculateMeanVelocity")
<< "can not find " << IF.name() << exit(FatalError);
}
else
{
scalarField UvdPlus(ny()), vPlus(ny()), wPlus(ny()), Ux(ny()), Uy(ny()), Uz(ny());
symmTensorField stressPlus(ny(), symmTensor::zero);
scalarField yPlus(ny()), y(ny());
scalar yInf, ySup;
scalar R(287.0), max_UvdPlus(26.457427), As(1.512e-06), Ts(120.0);
scalar r=pow(Prt(), 1.0/3.0);
//- calculate nuwInlet and Utao
scalar Taw = Tinf()*(1.0+(gamma()-1.0)/2.0*r*sqr(Ma()));
scalar Tw = (Taw - Tinf())/r + Tinf();
scalar rhow = pInf()/(R*Tw);
scalar muw = As*sqrt(Tw)/(1.0+Ts/Tw);
scalar nuwInlet = muw/rhow;
scalar A = sqrt((gamma()-1.0)/2.0*Prt()*sqr(Ma())*Tinf()/Tw);
scalar B = (1.0+sqrt(Prt())*(gamma()-1.0)/2.0*sqr(Ma()))*Tinf()/Tw - 1.0;
scalar Utao = Uinf()/A*(
asin((2.0*sqr(A)-B)/sqrt(sqr(B)+4.0*sqr(A)))
+ asin(B/sqrt(sqr(B)+4.0*sqr(A)))
)/max_UvdPlus; // max(Uvd+) = 25.6803733;
const scalarField T = this->patch().lookupPatchField<volScalarField, scalar>("T");
//- flag to decide whether to calculate deltaInlet
bool delta=true;
forAll(UvdPlus, patchI)
{
IF >> yPlus[patchI] >> UvdPlus[patchI] >> vPlus[patchI] >> wPlus[patchI]
>> stressPlus[patchI].xx() >> stressPlus[patchI].yy() >> stressPlus[patchI].zz() >> stressPlus[patchI].xy();
// calculate dimensional y
y[patchI] = yPlus[patchI]*nuwInlet/Utao;
// calculate dimensional Ux
scalar tmp = A*(UvdPlus[patchI]*Utao)/Uinf() - asin(B/sqrt(sqr(B)+4.0*sqr(A)));
Ux[patchI] = (sin(tmp)*sqrt(sqr(B)+4.0*sqr(A)) + B)/(2.0*sqr(A))*Uinf();
Uy[patchI] = vPlus[patchI]*Utao;
Uz[patchI] = wPlus[patchI]*Utao;
if (Ux[patchI] > 0.99*Uinf() && delta)
{
deltaInlet = y[patchI];
delta = false; // deltaInlet calculated only once;
}
}
vectorField Cf = this->patch().Cf();
Ubar = vectorField(Cf.size());
symmTensorField Reynold(Cf.size(), symmTensor::zero);
ReynoldStress = symmTensorField(Cf.size(), symmTensor::zero);
// interpolate at the cell
forAll(Cf, faceI)
{
for(label patchI=0; patchI < y.size()-1; patchI++)
{
yInf = Cf[faceI][1]-y[patchI];
ySup = Cf[faceI][1]-y[patchI+1];
if(yInf*ySup <= 0)
{
Ubar[faceI][0] = yInf/(yInf-ySup)*Ux[patchI+1] - ySup/(yInf-ySup)*Ux[patchI];
Ubar[faceI][1] = yInf/(yInf-ySup)*Uy[patchI+1] - ySup/(yInf-ySup)*Uy[patchI];
Ubar[faceI][2] = yInf/(yInf-ySup)*Uz[patchI+1] - ySup/(yInf-ySup)*Uz[patchI];
Reynold[faceI] = yInf/(yInf-ySup)*stressPlus[patchI+1] - ySup/(yInf-ySup)*stressPlus[patchI];
// calculate dimensional Reynold stress
scalar ksi = sqrt(Tw/T[faceI]);
ReynoldStress[faceI].xx() = sqr(Reynold[faceI].xx()/ksi)*Utao;
ReynoldStress[faceI].xy() = -sqr(Reynold[faceI].xy()/ksi)*Utao;
ReynoldStress[faceI].yy() = sqr(Reynold[faceI].yy()/ksi)*Utao;
ReynoldStress[faceI].zz() = sqr(Reynold[faceI].zz()/ksi)*Utao;
break;
}
else // Cf.component(1) is higher than y.last()
{
Ubar[faceI][0] = Ux.last();
Ubar[faceI][1] = Uy.last();
Ubar[faceI][2] = Uz.last();
scalar ksi = sqrt(Tw/T[faceI]);
ReynoldStress[faceI].xx() = sqr(stressPlus.last().xx()/ksi)*Utao;
ReynoldStress[faceI].xy() = -sqr(stressPlus.last().xy()/ksi)*Utao;
ReynoldStress[faceI].yy() = sqr(stressPlus.last().yy()/ksi)*Utao;
ReynoldStress[faceI].zz() = sqr(stressPlus.last().zz()/ksi)*Utao;
}
}
}
scalar ymin = min(Cf.component(1));
scalar Uxmin = min(Ubar.component(0));
scalar dUxdy = Uxmin/ymin;
scalar dUdy = sqr(Utao)/nuwInlet;
Info << "dUxdy (from case simulation) =" << dUxdy << nl
<< "dUdy (from theory calculation) =" << dUdy << nl
<< "Utao =" << Utao << nl
<< "nuwInlet =" << nuwInlet << nl
<< "deltaInlet =" << deltaInlet << endl;
}
}
我在单核下测试时时可以计算的,但是想并行计算的时候,decomposePar过不去,终端输出信息如下:
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.4.0 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 2.4.0-dcea1e13ff76
Exec : decomposePar
Date : Nov 13 2016
Time : 10:59:43
Host : "zwk"
PID : 12485
Case : /home/Ubuntu/WORKPLACE4/SyntheticCase-test
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Allowing user-supplied system call operations
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time
Decomposing mesh region0
Create mesh
Calculating distribution of cells
Selecting decompositionMethod simple
Finished decomposition in 0.04 s
Calculating original mesh data
Distributing cells to processors
Distributing faces to processors
Distributing points to processors
Constructing processor meshes
Processor 0
Number of cells = 23686
Number of faces shared with processor 1 = 1997
Number of faces shared with processor 1 = 58
Number of processor patches = 2
Number of processor faces = 2055
Number of boundary faces = 4243
Processor 1
Number of cells = 23686
Number of faces shared with processor 0 = 1997
Number of faces shared with processor 0 = 58
Number of faces shared with processor 2 = 1978
Number of faces shared with processor 2 = 52
Number of processor patches = 4
Number of processor faces = 4085
Number of boundary faces = 2189
Processor 2
Number of cells = 23686
Number of faces shared with processor 1 = 1978
Number of faces shared with processor 1 = 52
Number of faces shared with processor 3 = 1954
Number of faces shared with processor 3 = 11
Number of processor patches = 4
Number of processor faces = 3995
Number of boundary faces = 2329
Processor 3
Number of cells = 23686
Number of faces shared with processor 2 = 1954
Number of faces shared with processor 2 = 11
Number of faces shared with processor 4 = 1992
Number of faces shared with processor 4 = 57
Number of processor patches = 4
Number of processor faces = 4014
Number of boundary faces = 2276
Processor 4
Number of cells = 23686
Number of faces shared with processor 3 = 1992
Number of faces shared with processor 3 = 57
Number of faces shared with processor 5 = 1964
Number of faces shared with processor 5 = 43
Number of processor patches = 4
Number of processor faces = 4056
Number of boundary faces = 2208
Processor 5
Number of cells = 23686
Number of faces shared with processor 4 = 1964
Number of faces shared with processor 4 = 43
Number of faces shared with processor 6 = 1984
Number of faces shared with processor 6 = 41
Number of processor patches = 4
Number of processor faces = 4032
Number of boundary faces = 2330
Processor 6
Number of cells = 23685
Number of faces shared with processor 5 = 1984
Number of faces shared with processor 5 = 41
Number of faces shared with processor 7 = 1986
Number of faces shared with processor 7 = 55
Number of processor patches = 4
Number of processor faces = 4066
Number of boundary faces = 2218
Processor 7
Number of cells = 23685
Number of faces shared with processor 6 = 1986
Number of faces shared with processor 6 = 55
Number of processor patches = 2
Number of processor faces = 2041
Number of boundary faces = 4211
Number of processor faces = 14172
Max number of cells = 23686 (0.00105549% above average 23685.8)
Max number of processor patches = 4 (14.2857% above average 3.5)
Max number of faces between processors = 4085 (15.2978% above average 3543)
Time = 0
dUxdy (from case simulation) =2.46405e+06
dUdy (from theory calculation) =2.46471e+06
Utao =23.1659
nuwInlet =0.000217736
deltaInlet =0.0144036
Processor 0: field transfer
terminate called after throwing an instance of 'std::logic_error'
what(): basic_string::_S_construct null not valid
还麻烦各位可以给个意见,这个问题就是这个边条造成的,但不知道如何解决,谢谢各位了