Matlab mex memory handling - same input different result

2012 6 16
4 个回答
回答已采纳
5 次查看（30 天）

0 个投票

I ve converted a computationally expensive part of my Matlab code into a C mex file, which greatly speeded up the calculation. The problem is that in some runs the C mex code gives exactly the same result as the Matlab code (within the machine accuracy), whereas in some others there is a significant discrepancy between the two (reaching even 1e-3). Furthermore, there are times that the C mex code disagrees with itself by giving in consecutive runs different results with exactly the same inputs. I guess that the problem is probably caused by some incorrect memory handling due to my programming. This is highly likely as thats my first attempt to create a C mex file.
I ve read somewhere that differences between the mex function and Matlab are normal, especially when using trigonometric functions like acos (which I am using). However does this justify a different result even between consecutive runs with the same input???
I attach my mexFunction which calls another C function (in the same file) to perform the calculations. I would be grateful if you could give me some insight on my problem. Please let me know if you need any other part of the code.
Thanks in advance,
Panos
/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
mwIndex InterestBE;
mwIndex BladeElement;
double *WakeXYZ;
double *RQuarter;
double *BE_RControl;
double *BE_Gamma;
double *GCBladeElement;
const mwSize *dimsWakeXYZ;
const mwSize *dimsBE_RControl;
mwSize NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems;
mwSize NoOfElements, NoOfDimensions;
double p_CompressibleBiot;
double *p_M;
double *p_CompressibleWakeMachMax;
mwSize p_WakePointsNo;
double p_BoundVorticity;
double *p_VortexCoreOffset;
double *p_VortexDelta;
double *p_kin_visc;
double *p_VatistasN;
//----------------------------------------------------------
// First retrieve inputs
//----------------------------------------------------------
/* Get the value of the scalar input InterestBE */
InterestBE = (mwIndex)mxGetScalar(prhs[0]);
/* Get the value of the scalar input Filament */
BladeElement = (mwIndex)mxGetScalar(prhs[1]);
/* create a pointer to the real data in the WakeXYZ matrix */
WakeXYZ = mxGetPr(prhs[2]);
/* create a pointer to the real data in the RQuarter matrix */
RQuarter = mxGetPr(prhs[3]);
/* create a pointer to the real data in the BE_RControl matrix */
BE_RControl = mxGetPr(prhs[4]);
/* create a pointer to the real data in the BE_RControl matrix */
BE_Gamma = mxGetPr(prhs[5]);
/* create pointers to the fields of the conf structure */
p_CompressibleBiot = mxGetScalar(mxGetField(prhs[6], 0, "CompressibleBiot"));
p_M = mxGetPr(mxGetField(prhs[6], 0, "M"));
p_CompressibleWakeMachMax = mxGetPr(mxGetField(prhs[6], 0, "CompressibleWakeMachMax"));
p_VortexCoreOffset = mxGetPr(mxGetField(prhs[6], 0, "VortexCoreOffset"));
p_VortexDelta = mxGetPr(mxGetField(prhs[6], 0, "VortexDelta"));
p_kin_visc = mxGetPr(mxGetField(prhs[6], 0, "kin_visc"));
p_VatistasN = mxGetPr(mxGetField(prhs[6], 0, "VatistasN"));
p_WakePointsNo = (mwSize)mxGetScalar(mxGetField(prhs[6],0,"WakePoints"));
p_BoundVorticity = mxGetScalar(mxGetField(prhs[6],0,"BoundVorticity"));
// Get the dimensions of the WakeXYZ array
dimsWakeXYZ = mxGetDimensions(prhs[2]);
// Get the dimensions of the BE_RControl array
dimsBE_RControl = mxGetDimensions(prhs[4]);
NoOfBlades = dimsWakeXYZ[0];
NoOfFilaments = dimsWakeXYZ[1];
NoOfSegments = dimsWakeXYZ[2];
NoOfWakeItems = dimsWakeXYZ[3];
NoOfElements = dimsBE_RControl[1];
//----------------------------------------------------------
// Call the calculation subroutine
//----------------------------------------------------------
/* create the output matrix */
plhs[0] = mxCreateDoubleMatrix(1,3,mxREAL);
GCBladeElement = mxGetPr(plhs[0]);
BladeElementGC(InterestBE, BladeElement,
WakeXYZ, RQuarter, BE_RControl, BE_Gamma,
NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems,
NoOfElements,
p_CompressibleBiot, *p_M, *p_CompressibleWakeMachMax,
*p_VortexCoreOffset,
*p_VortexDelta,
*p_kin_visc,
*p_VatistasN,
p_BoundVorticity,
GCBladeElement);
}

请先登录，再回答此问题。

 采纳的回答

Panagiotis
Panagiotis 2012-6-18

0 个投票

Hello again,
thank you so much for confirming that you couldn't see any obvious mistake. That made me think that probably there was another kind of bug and I ran extensive tests by comparing Matlab and C results. Finally, there was a series of small subtle bugs, like square roots of negative numbers, which in Matlab give a complex result but in C they return nan etc. I really wonder how I could get some "correct" results before....
once again thanks for your time
Panos

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Titus Edelhofer
Titus Edelhofer 2012-6-19
Hi Panos,
it's good to hear the problem is solved. You might mark the question as answered then ...
Titus

请先登录，再进行评论。

更多回答（3 个）

Titus Edelhofer
Titus Edelhofer 2012-6-16

0 个投票

Hi,
looking on the code I don't see a reason why consecutive calls don't give the same answer. So it looks as if the answer to your question lies within "BladeElemtGC"...
Titus
James Tursa
James Tursa 2012-6-16

0 个投票

Cross-posted in the newsgroup, so I will repeat my reply here:
What does the BladeElementGC interface look like? Normally I would not expect you to be dereferenceing the pointers in the argument list. I.e., why are you passing *p_kin_visc instead of just p_kin_visc? If you are accessing array elements inside the BladeElementGC routine then dereferencing the pointers in the argument list will not work. But again, I would need to see the actual interface to be sure what is going wrong.
Panagiotis
Panagiotis 2012-6-17

0 个投票

Titus and James,
thank you very much for your prompt replies and your willingness to help.
James I am mainly a Fortran programmer, so I am not very experienced with C. I have dereferenced the specific arguments as they are simple scalars and I feel more comfortable working with "non-pointer" variables due to my Fortran habits. Do you think that this might be a source of problems?
thank you once more Panos
Sorry it's quite long but, here is the rest of the code. The BladeElementGC routine is at the end.
#include <mex.h>
#include <math.h>
#include <stdio.h>
#include <matrix.h>
#define wakeXYZ(i,j,k,l) WakeXYZ[(i-1)+((j-1)+((k-1)+(l-1)*NoOfSegments)*NoOfFilaments)*NoOfBlades]
#define be_RControl(i,j,k) BE_RControl[(i-1)+((j-1)+(k-1)*NoOfElements)*NoOfBlades]
#define rquarter(i,j,k) RQuarter[(i-1)+((j-1)+(k-1)*NoOfElements)*NoOfBlades]
#define be_Gamma(i) BE_Gamma[(i-1)]
/* The computational routine for FarTrailingSingleGC */
void FarTrailingSingleGC(const mwIndex InterestBE, const mwIndex Segment, const mwIndex Filament, const mwIndex BladeNo,
const double *WakeXYZ, const double *BE_RControl,
const mwSize NoOfBlades, const mwSize NoOfFilaments, const mwSize NoOfSegments, const mwSize NoOfWakeItems,
const mwSize NoOfElements,
const double conf_CompressibleBiot, const double conf_M, const double conf_CompressibleWakeMachMax,
double *GC)
{
.......
Similar to TrailingSingleGC
.......
}
/* The computational routine for TrailingSingleGC */
void TrailingSingleGC(const mwIndex InterestBE, const mwIndex Segment, const mwIndex Filament, const mwIndex BladeNo,
const double *WakeXYZ, const double *BE_RControl, const double *BE_Gamma,
const mwSize NoOfBlades, const mwSize NoOfFilaments, const mwSize NoOfSegments, const mwSize NoOfWakeItems,
const mwSize NoOfElements,
const double conf_CompressibleBiot, const double conf_M, const double conf_CompressibleWakeMachMax,
const double conf_VortexCoreOffset,
const double conf_VortexDelta,
const double conf_kin_visc,
const double conf_VatistasN,
double *GC)
{
const mwSize NoOfDimensions = 3;
double rA[3];
double rB[3];
double rp[3];
double Mach;
double CompreBeta;
double CompreMach;
double rm[3];
double rmp[3];
double lAB[3];
double r1[3];
double r2[3];
double lAB_dot_r1;
double lAB_dot_r2;
double lABxr1[3];
double cosTheta1;
double cosTheta2;
double height;
double tbase;
double BioConstant;
const double pi = 3.1415926535897932384626433832795;
//-----------------------------------------------------------
// Vortex core calculations
//-----------------------------------------------------------
double Omegam, AzimuthStart, vDelta, TrailGamma, Core0, MidAge, CoreRadius;
const double CoreAlpha = 1.25643;
// Calculate the segment center point velocity vector
Omegam = (wakeXYZ(BladeNo,Filament,Segment,6) + wakeXYZ(BladeNo,Filament,Segment,6))/2.;
// Calculate the vortex core radius
AzimuthStart = -conf_VortexCoreOffset;
if (be_Gamma(1)<-800.)
{
vDelta = -conf_VortexDelta;
}
else
{
if (Filament ==1)
{
TrailGamma = fabs(-be_Gamma(1));
}
else if (Filament ==NoOfElements+1)
{
TrailGamma = fabs(be_Gamma(NoOfElements));
}
else
{
TrailGamma = fabs(be_Gamma(Filament-1)-be_Gamma(Filament));
}
vDelta = 1. + 2.*pow(10.,-4.) * TrailGamma / conf_kin_visc;
}
Core0 = sqrt(4.*CoreAlpha*vDelta*conf_kin_visc*(-AzimuthStart/Omegam));
MidAge = (wakeXYZ(BladeNo,Filament,Segment,4) + wakeXYZ(BladeNo,Filament,Segment,4))/2.;
CoreRadius = sqrt(pow(Core0,2.) + 4.*CoreAlpha*conf_VortexDelta*conf_kin_visc*MidAge);
//-----------------------------------------------------------
// Calculate the compressibility correction factor
//-----------------------------------------------------------
if (conf_CompressibleBiot>=2)
{
if (conf_M > conf_CompressibleWakeMachMax)
{
CompreMach = conf_CompressibleWakeMachMax;
}
else
{
CompreMach = conf_M;
}
}
else
{
CompreMach = 0.;
}
CompreBeta = 1./sqrt(1.-pow(CompreMach,2.));
/* Get the rA position vector */
rA[0] = wakeXYZ(BladeNo,Filament,Segment,1);
rA[1] = wakeXYZ(BladeNo,Filament,Segment,2);
rA[2] = CompreBeta * wakeXYZ(BladeNo,Filament,Segment,3);
/* Get the rB position vector */
rB[0] = wakeXYZ(BladeNo,Filament,Segment+1,1);
rB[1] = wakeXYZ(BladeNo,Filament,Segment+1,2);
rB[2] = CompreBeta * wakeXYZ(BladeNo,Filament,Segment+1,3);
/* Get the rp position vector */
rp[0] = be_RControl(1,InterestBE,1);
rp[1] = be_RControl(1,InterestBE,2);
rp[2] = CompreBeta * be_RControl(1,InterestBE,3);
// Calculate the AB mean point
rm[0] = (rA[0] + rB[0]) / 2.;
rm[1] = (rA[1] + rB[1]) / 2.;
rm[2] = (rA[2] + rB[2]) / 2.;
// Calculate the distance from the center point
rmp[0] = rp[0] - rm[0];
rmp[1] = rp[1] - rm[1];
rmp[2] = rp[2] - rm[2];
// Calculate the segment length vector
lAB[0] = rB[0] - rA[0];
lAB[1] = rB[1] - rA[1];
lAB[2] = rB[2] - rA[2];
// Calculate the r1=rAP vector
r1[0] = rp[0] - rA[0];
r1[1] = rp[1] - rA[1];
r1[2] = rp[2] - rA[2];
// Calculate the r2=rBP vector
r2[0] = rp[0] - rB[0];
r2[1] = rp[1] - rB[1];
r2[2] = rp[2] - rB[2];
// Calculate the lAB r1 dot product
lAB_dot_r1 = lAB[0]*r1[0] + lAB[1]*r1[1] + lAB[2]*r1[2];
// Calculate the lAB r2 dot product
lAB_dot_r2 = lAB[0]*r2[0] + lAB[1]*r2[1] + lAB[2]*r2[2];
// Now calculate the lAB x r1 cross product vector
lABxr1[0] = lAB[1]*r1[2] - lAB[2]*r1[1];
lABxr1[1] = lAB[2]*r1[0] - lAB[0]*r1[2];
lABxr1[2] = lAB[0]*r1[1] - lAB[1]*r1[0];
// Calculate the cosTheta1
cosTheta1 = lAB_dot_r1 / ( sqrt(pow(lAB[0],2.) + pow(lAB[1],2.) + pow(lAB[2],2.)) * sqrt(pow(r1[0],2.) + pow(r1[1],2.) + pow(r1[2],2.)) );
// Calculate the cosTheta2
cosTheta2 = lAB_dot_r2 / ( sqrt(pow(lAB[0],2.) + pow(lAB[1],2.) + pow(lAB[2],2.)) * sqrt(pow(r2[0],2.) + pow(r2[1],2.) + pow(r2[2],2.)) );
// Calculate the triangle height
tbase = sqrt(pow(r1[0],2.) + pow(r1[1],2.) + pow(r1[2],2.)) * cosTheta1;
height = sqrt(pow(r1[0],2.) + pow(r1[1],2.) + pow(r1[2],2.) - pow(tbase,2.));
// Finaly calculate the geometry coefficient by the Biot-Savart law
if (conf_VortexCoreOffset<0)
{
// No vortex core is used
//BioConstant = (cosTheta1-cosTheta2)/(4.*pi*height) / sqrt(pow(lABxr1[0],2.)+pow(lABxr1[1],2.)+pow(lABxr1[2],2.));
BioConstant = cosTheta1/(4.*pi*height) / sqrt(pow(lABxr1[0],2.)+pow(lABxr1[1],2.)+pow(lABxr1[2],2.))
-cosTheta2/(4.*pi*height) / sqrt(pow(lABxr1[0],2.)+pow(lABxr1[1],2.)+pow(lABxr1[2],2.));
GC[0] = BioConstant * lABxr1[0];
GC[1] = BioConstant * lABxr1[1];
GC[2] = BioConstant * lABxr1[2];
}
else
{
// Calculation includes vortex core
//BioConstant = height/pow((pow(CoreRadius,(2.*conf_VatistasN)) + pow(height,(2.*conf_VatistasN))),(1./conf_VatistasN))*
// (cosTheta1-cosTheta2)/(4.*pi) / sqrt(pow(lABxr1[0],2.)+pow(lABxr1[1],2.)+pow(lABxr1[2],2.));
BioConstant = height/pow((pow(CoreRadius,(2.*conf_VatistasN)) + pow(height,(2.*conf_VatistasN))),(1./conf_VatistasN))*
cosTheta1/(4.*pi) / sqrt(pow(lABxr1[0],2.)+pow(lABxr1[1],2.)+pow(lABxr1[2],2.))
- height/pow((pow(CoreRadius,(2.*conf_VatistasN)) + pow(height,(2.*conf_VatistasN))),(1./conf_VatistasN))*
cosTheta2/(4.*pi) / sqrt(pow(lABxr1[0],2.)+pow(lABxr1[1],2.)+pow(lABxr1[2],2.));
GC[0] = BioConstant * lABxr1[0];
GC[1] = BioConstant * lABxr1[1];
GC[2] = BioConstant * lABxr1[2];
}
// A simple desingularisation
if ( isinf(GC[0]) || isnan(GC[0]) ) {
GC[0] = 0.;
}
if ( isinf(GC[1]) || isnan(GC[1]) ) {
GC[1] = 0.;
}
if ( isinf(GC[2]) || isnan(GC[2]) ) {
GC[2] = 0.;
}
}
/* The computational routine for BoundSingleGC */
void BoundSingleGC(const mwIndex InterestBE, const mwIndex BladeElement, const mwIndex BladeNo,
const double *WakeXYZ, const double *RQuarter, const double *BE_RControl, const double *BE_Gamma,
const mwSize NoOfBlades, const mwSize NoOfFilaments, const mwSize NoOfSegments, const mwSize NoOfWakeItems,
const mwSize NoOfElements,
const double conf_CompressibleBiot, const double conf_M, const double conf_CompressibleWakeMachMax,
const double conf_VortexCoreOffset,
const double conf_VortexDelta,
const double conf_kin_visc,
const double conf_VatistasN,
double *GC)
{
const mwSize NoOfDimensions = 3;
double rA[3];
double rB[3];
double rp[3];
double Mach;
double CompreBeta;
double CompreMach;
double rm[3];
double rmp[3];
double lAB[3];
double r1[3];
double r2[3];
double lAB_dot_r1;
double lAB_dot_r2;
double lABxr1[3];
double cosTheta1;
double cosTheta2;
double height;
double tbase;
double BioConstant;
const double pi = 3.1415926535897932384626433832795;
//-----------------------------------------------------------
// Vortex core calculations
//-----------------------------------------------------------
double Omegam, AzimuthStart, vDelta, TrailGamma, Core0, MidAge, CoreRadius;
const double CoreAlpha = 1.25643;
// Calculate the corresponding trailing vortex omega
Omegam = (wakeXYZ(BladeNo,NoOfElements+1,1,6) + wakeXYZ(BladeNo,NoOfElements+1,1,6))/2.;
// Calculate the vortex core radius
AzimuthStart = -conf_VortexCoreOffset;
// Calculate the bound vortex radius as equal to the starting radius
// of the corresponding trailing vortex
if (be_Gamma(1)<-800.)
{
vDelta = -conf_VortexDelta;
}
else
{
TrailGamma = fabs(be_Gamma(BladeElement));
vDelta = 1. + 2.*pow(10.,-4.) * TrailGamma / conf_kin_visc;
}
Core0 = sqrt(4.*CoreAlpha*vDelta*conf_kin_visc*(-AzimuthStart/Omegam));
CoreRadius = Core0;
//-----------------------------------------------------------
// Calculate the compressibility correction factor
//-----------------------------------------------------------
if (conf_CompressibleBiot>=2)
{
if (conf_M > conf_CompressibleWakeMachMax)
{
CompreMach = conf_CompressibleWakeMachMax;
}
else
{
CompreMach = conf_M;
}
}
else
{
CompreMach = 0.;
}
CompreBeta = 1./sqrt(1.-pow(CompreMach,2.));
/* Get the rA position vector */
rA[0] = rquarter(BladeNo,BladeElement,1);
rA[1] = rquarter(BladeNo,BladeElement,2);
rA[2] = CompreBeta * rquarter(BladeNo,BladeElement,3);
/* Get the rB position vector */
rB[0] = rquarter(BladeNo,BladeElement+1,1);
rB[1] = rquarter(BladeNo,BladeElement+1,2);
rB[2] = CompreBeta * rquarter(BladeNo,BladeElement+1,3);
/* Get the rp position vector */
rp[0] = be_RControl(1,InterestBE,1);
rp[1] = be_RControl(1,InterestBE,2);
rp[2] = CompreBeta * be_RControl(1,InterestBE,3);
// Calculate the AB mean point
rm[0] = (rA[0] + rB[0]) / 2.;
rm[1] = (rA[1] + rB[1]) / 2.;
rm[2] = (rA[2] + rB[2]) / 2.;
// Calculate the distance from the center point
rmp[0] = rp[0] - rm[0];
rmp[1] = rp[1] - rm[1];
rmp[2] = rp[2] - rm[2];
// Calculate the segment length vector
lAB[0] = rB[0] - rA[0];
lAB[1] = rB[1] - rA[1];
lAB[2] = rB[2] - rA[2];
// Calculate the r1=rAP vector
r1[0] = rp[0] - rA[0];
r1[1] = rp[1] - rA[1];
r1[2] = rp[2] - rA[2];
// Calculate the r2=rBP vector
r2[0] = rp[0] - rB[0];
r2[1] = rp[1] - rB[1];
r2[2] = rp[2] - rB[2];
// Calculate the lAB r1 dot product
lAB_dot_r1 = lAB[0]*r1[0] + lAB[1]*r1[1] + lAB[2]*r1[2];
// Calculate the lAB r2 dot product
lAB_dot_r2 = lAB[0]*r2[0] + lAB[1]*r2[1] + lAB[2]*r2[2];
// Now calculate the lAB x r1 cross product vector
lABxr1[0] = lAB[1]*r1[2] - lAB[2]*r1[1];
lABxr1[1] = lAB[2]*r1[0] - lAB[0]*r1[2];
lABxr1[2] = lAB[0]*r1[1] - lAB[1]*r1[0];
// Calculate the cosTheta1
cosTheta1 = lAB_dot_r1 / ( sqrt(pow(lAB[0],2.) + pow(lAB[1],2.) + pow(lAB[2],2.)) * sqrt(pow(r1[0],2.) + pow(r1[1],2.) + pow(r1[2],2.)) );
// Calculate the cosTheta2
cosTheta2 = lAB_dot_r2 / ( sqrt(pow(lAB[0],2.) + pow(lAB[1],2.) + pow(lAB[2],2.)) * sqrt(pow(r2[0],2.) + pow(r2[1],2.) + pow(r2[2],2.)) );
// Calculate the triangle height
tbase = sqrt(pow(r1[0],2.) + pow(r1[1],2.) + pow(r1[2],2.)) * cosTheta1;
height = sqrt(pow(r1[0],2.) + pow(r1[1],2.) + pow(r1[2],2.) - pow(tbase,2.));
// Finaly calculate the geometry coefficient by the Biot-Savart law
if (conf_VortexCoreOffset<0)
{
// No vortex core is used
//BioConstant = (cosTheta1-cosTheta2)/(4.*pi*height) / sqrt(pow(lABxr1[0],2.)+pow(lABxr1[1],2.)+pow(lABxr1[2],2.));
BioConstant = cosTheta1/(4.*pi*height) / sqrt(pow(lABxr1[0],2.)+pow(lABxr1[1],2.)+pow(lABxr1[2],2.))
-cosTheta2/(4.*pi*height) / sqrt(pow(lABxr1[0],2.)+pow(lABxr1[1],2.)+pow(lABxr1[2],2.));
GC[0] = BioConstant * lABxr1[0];
GC[1] = BioConstant * lABxr1[1];
GC[2] = BioConstant * lABxr1[2];
}
else
{
// Calculation includes vortex core
//BioConstant = height/pow((pow(CoreRadius,(2.*conf_VatistasN)) + pow(height,(2.*conf_VatistasN))),(1./conf_VatistasN))*
// (cosTheta1-cosTheta2)/(4.*pi) / sqrt(pow(lABxr1[0],2.)+pow(lABxr1[1],2.)+pow(lABxr1[2],2.));
BioConstant = height/pow((pow(CoreRadius,(2.*conf_VatistasN)) + pow(height,(2.*conf_VatistasN))),(1./conf_VatistasN))*
cosTheta1/(4.*pi) / sqrt(pow(lABxr1[0],2.)+pow(lABxr1[1],2.)+pow(lABxr1[2],2.))
- height/pow((pow(CoreRadius,(2.*conf_VatistasN)) + pow(height,(2.*conf_VatistasN))),(1./conf_VatistasN))*
cosTheta2/(4.*pi) / sqrt(pow(lABxr1[0],2.)+pow(lABxr1[1],2.)+pow(lABxr1[2],2.));
GC[0] = BioConstant * lABxr1[0];
GC[1] = BioConstant * lABxr1[1];
GC[2] = BioConstant * lABxr1[2];
}
if ((InterestBE==BladeElement) && (BladeNo==1))
{
GC[0] = 0.;
GC[1] = 0.;
GC[2] = 0.;
}
// A simple desingularisation
if ( isinf(GC[0]) || isnan(GC[0]) ) {
GC[0] = 0.;
}
if ( isinf(GC[1]) || isnan(GC[1]) ) {
GC[1] = 0.;
}
if ( isinf(GC[2]) || isnan(GC[2]) ) {
GC[2] = 0.;
}
}
/* The computational routine for TrailingFilamentGC */
void TrailingFilamentGC(const mwIndex InterestBE, const mwIndex Filament,
const double *WakeXYZ, const double *BE_RControl, const double *BE_Gamma,
const mwSize NoOfBlades, const mwSize NoOfFilaments, const mwSize NoOfSegments, const mwSize NoOfWakeItems,
const mwSize NoOfElements,
const double conf_CompressibleBiot, const double conf_M, const double conf_CompressibleWakeMachMax,
const double conf_VortexCoreOffset,
const double conf_VortexDelta,
const double conf_kin_visc,
const double conf_VatistasN,
double *GCTRfilament)
{
double GCsingle[3];
mwIndex Segment;
mwIndex BladeNo;
// Initialise the filament Geometric Coefficient
GCTRfilament[0] = 0.;
GCTRfilament[1] = 0.;
GCTRfilament[2] = 0.;
// Add influence for all blades and segments
for (BladeNo = 1; BladeNo <= NoOfBlades; BladeNo++)
{
for (Segment = 1; Segment <= NoOfSegments-1; Segment++)
{
GCsingle[0] = 0.;
GCsingle[1] = 0.;
GCsingle[2] = 0.;
/* call the single trailing segment routine */
if (Segment<=20)
{
TrailingSingleGC(InterestBE, Segment, Filament, BladeNo,
WakeXYZ, BE_RControl, BE_Gamma,
NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems,
NoOfElements,
conf_CompressibleBiot, conf_M, conf_CompressibleWakeMachMax,
conf_VortexCoreOffset,
conf_VortexDelta,
conf_kin_visc,
conf_VatistasN,
GCsingle);
}
else
{
FarTrailingSingleGC(InterestBE, Segment, Filament, BladeNo,
WakeXYZ, BE_RControl,
NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems,
NoOfElements,
conf_CompressibleBiot, conf_M, conf_CompressibleWakeMachMax,
GCsingle);
}
// Add the segment influence
GCTRfilament[0] += GCsingle[0];
GCTRfilament[1] += GCsingle[1];
GCTRfilament[2] += GCsingle[2];
}
}
}
/* The computational routine for BoundElementGC */
void BoundElementGC(const mwIndex InterestBE, const mwIndex BladeElement,
const double *WakeXYZ, const double *RQuarter, const double *BE_RControl, const double *BE_Gamma,
const mwSize NoOfBlades, const mwSize NoOfFilaments, const mwSize NoOfSegments, const mwSize NoOfWakeItems,
const mwSize NoOfElements,
const double conf_CompressibleBiot, const double conf_M, const double conf_CompressibleWakeMachMax,
const double conf_VortexCoreOffset,
const double conf_VortexDelta,
const double conf_kin_visc,
const double conf_VatistasN,
double *GCBoundElement)
{
double GCsingle[3];
mwIndex BladeNo;
// Initialise the bound element Geometric Coefficient
GCBoundElement[0] = 0.;
GCBoundElement[1] = 0.;
GCBoundElement[2] = 0.;
// Add influence for all blades
for (BladeNo = 1; BladeNo <= NoOfBlades; BladeNo++)
{
GCsingle[0] = 0.;
GCsingle[1] = 0.;
GCsingle[2] = 0.;
/* call the single bound segment routine */
BoundSingleGC(InterestBE, BladeElement, BladeNo,
WakeXYZ, RQuarter, BE_RControl, BE_Gamma,
NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems,
NoOfElements,
conf_CompressibleBiot, conf_M, conf_CompressibleWakeMachMax,
conf_VortexCoreOffset,
conf_VortexDelta,
conf_kin_visc,
conf_VatistasN,
GCsingle);
// Add the segment influence
GCBoundElement[0] += GCsingle[0];
GCBoundElement[1] += GCsingle[1];
GCBoundElement[2] += GCsingle[2];
}
}
/* The computational routine for BoundElementGC */
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Calculate the Geometric Coefficient of a blade element on an
// Interest blade element
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
void BladeElementGC(const mwIndex InterestBE, const mwIndex BladeElement,
const double *WakeXYZ, const double *RQuarter, const double *BE_RControl, const double *BE_Gamma,
const mwSize NoOfBlades, const mwSize NoOfFilaments, const mwSize NoOfSegments, const mwSize NoOfWakeItems,
const mwSize NoOfElements,
const double conf_CompressibleBiot, const double conf_M, const double conf_CompressibleWakeMachMax,
const double conf_VortexCoreOffset,
const double conf_VortexDelta,
const double conf_kin_visc,
const double conf_VatistasN,
const double conf_BoundVorticity,
double *GCBladeElement)
{
mwIndex Filament;
double Temp1[3], Temp2[3], Temp3[3];
GCBladeElement[0] = 0.;
GCBladeElement[1] = 0.;
GCBladeElement[2] = 0.;
Filament = BladeElement;
TrailingFilamentGC(InterestBE, Filament,
WakeXYZ, BE_RControl, BE_Gamma,
NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems,
NoOfElements,
conf_CompressibleBiot, conf_M, conf_CompressibleWakeMachMax,
conf_VortexCoreOffset,
conf_VortexDelta,
conf_kin_visc,
conf_VatistasN,
Temp1);
TrailingFilamentGC(InterestBE, Filament+1,
WakeXYZ, BE_RControl, BE_Gamma,
NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems,
NoOfElements,
conf_CompressibleBiot, conf_M, conf_CompressibleWakeMachMax,
conf_VortexCoreOffset,
conf_VortexDelta,
conf_kin_visc,
conf_VatistasN,
Temp2);
if (conf_BoundVorticity>0)
{
BoundElementGC(InterestBE, BladeElement,
WakeXYZ, RQuarter, BE_RControl, BE_Gamma,
NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems,
NoOfElements,
conf_CompressibleBiot, conf_M, conf_CompressibleWakeMachMax,
conf_VortexCoreOffset,
conf_VortexDelta,
conf_kin_visc,
conf_VatistasN,
Temp3);
}
else
{
Temp3[0] = 0.;
Temp3[1] = 0.;
Temp3[2] = 0.;
}
GCBladeElement[0] = Temp2[0] - Temp1[0] + Temp3[0];
GCBladeElement[1] = Temp2[1] - Temp1[1] + Temp3[1];
GCBladeElement[2] = Temp2[2] - Temp1[2] + Temp3[2];
}

2 个评论
显示 无 隐藏 无

Titus Edelhofer
Titus Edelhofer 2012-6-18
Hmm, this will not be too easy. I guess, using the debugger/print statements to find out where the results start to differ will be (the only?) option. Looking at the code it looks "o.k." to me. One thing I observed: what do the lines like
Omegam = (wakeXYZ(BladeNo,Filament,Segment,6) + wakeXYZ(BladeNo,Filament,Segment,6))/2.;
do? It looks as if the 2 numbers you add are exactly the same? In this case Omegam=wakeXYZ(BladeNo,Filament,Segment,6); ? Happens later with MidAgm again.
Titus
Panagiotis
Panagiotis 2012-6-18
Well-spotted Titus thanks a lot
I did as you said and sorted it out finally. I was initially reluctant to do so first because I didn't know whether my C mex coding was right and secondly because the subroutine is called thousands of times and the errors only occurred in some of the runs.
Panos

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Write C Functions Callable from MATLAB (MEX Files) 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by