Matlab mex memory handling - same input different result

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);

   }

Answer 1

[Panagiotis]

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

[Titus Edelhofer]

Hi Panos,

it's good to hear the problem is solved. You might mark the question as answered then ...

Titus

Answer 2

[Titus Edelhofer]

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

Answer 3

[James Tursa]

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.

Answer 4

[Panagiotis]

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];

}

[Titus Edelhofer]

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]

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