[petsc-users] petsc error disappears when I print something in the function

[Kaushik Vijaykumar]

Thanks Dave and Mathew. I will keep looking for the bug and here is the
full piece of code.

PetscInt kel_tet10CX(PetscScalar *X, PetscScalar *Y, PetscScalar *Z,
PetscScalar *ke, PetscInt elmat, PetscScalar dmg, PetscInt el_has_orient,
PetscScalar rot[3][3], PetscScalar temperature, PetscScalar *thermforce,
PetscInt damflag, PetscScalar Cmat[6][6], PetscScalar alphavec[6]){
/////////PETSc declarations
PetscErrorCode ierr;
PetscInt ii,jj,kk,k1,ll,l1,i,j,k,l,i1,j1,jj1,ii1,iflag;
PetscInt nope=10;   //// nodes on the element
PetscScalar
alpha2[4][4],alpha2unrot[4][4],C[6][6],C4[4][4][4][4],C4unrot[4][4][4][4],eta,etas[9],GP[4],s[100][100],shp[10][4],shpj[4][11],w[4][4],weight,wts[4],xi,xis[4],xl[10][3],xsj,xsjj,zeta,zetas[4];

/////////////Prototypes
PetscInt constit(PetscInt elmat,PetscScalar C4[4][4][4][4], PetscScalar
alpha2[4][4], PetscScalar dmg, PetscScalar temperature, PetscInt damflag,
PetscScalar Cmat[6][6], PetscScalar alphavec[6]);
void shape10tet_(PetscScalar *xi, PetscScalar *eta, PetscScalar *zeta,
PetscScalar xl[4][11], PetscScalar *xsj, PetscScalar shp[5][11], PetscInt
*iflag); // you have to dimension the array arguments to the size Fortran
sees them

// printf("\n Entered SPAF element Tet10 \n");
/////////////Code
if (el_has_orient == 1){ierr =
constit(elmat,C4unrot,alpha2unrot,dmg,temperature,damflag,Cmat,alphavec);ierr
= rot4(rot,C4unrot,C4);ierr = rot2(rot,alpha2unrot,alpha2);}
else {ierr =
constit(elmat,C4,alpha2,dmg,temperature,damflag,Cmat,alphavec);}

// printf("\n Came back from spaf materials to SPAF element Tet10 \n");
// ierr = PetscPrintf(PETSC_COMM_SELF,"In tet10 \n");CHKERRQ(ierr);


// Gauss points and weights
xis[0] = 0.138196601125011;
xis[1] = 0.585410196624968;
xis[2] = 0.138196601125011;
xis[3] = 0.138196601125011;
etas[0] = 0.138196601125011;
etas[1] = 0.138196601125011;
etas[2] = 0.585410196624968;
etas[3] = 0.138196601125011;
zetas[0] = 0.138196601125011;
zetas[1] = 0.138196601125011;
zetas[2] = 0.138196601125011;
zetas[3] = 0.585410196624968;
wts[0] = 0.04166666666667;
wts[1] = 0.04166666666667;
wts[2] = 0.04166666666667;
wts[3] = 0.04166666666667;

// assemble nodal coordinates into xl
// when fortran gets it, it will be transposed and the indices increased by
1
for (ii=0; ii<nope; ii++)
{
xl[ii][0] = X[ii];//fprintf(outfile,"X[%i] = %f\n",ii,X[ii]);
xl[ii][1] = Y[ii];//fprintf(outfile,"Y[%i] = %f\n",ii,Y[ii]);
xl[ii][2] = Z[ii];//fprintf(outfile,"Z[%i] = %f\n",ii,Z[ii]);
}

// ke
memset(ke, 0, sizeof(ke[0])*nope*3*nope*3);  // why doesn't compiler know
the real size of ke?
memset(s, 0, sizeof(s));
iflag = 3; /// tells shape routine to calculate shape functions,
derivatives and det J

  for (kk=0; kk<4; kk++) // Loop 3
  {
  weight = wts[kk];
  // Integration point
  xi = xis[kk];
  eta = etas[kk];
  zeta = zetas[kk];
  shape10tet_(&xi,&eta,&zeta,xl,&xsj,shp,&iflag);
       if (xsj<0){fprintf(outfile,"Negative det J on TET10 element!!
(Integration point %i)\n",kk);}
  xsjj = sqrt(xsj);
  for (i1=1; i1<nope+1; i1++) // Loop 4a
  {
    shpj[1][i1] = shp[i1-1][0]*xsjj; // shape function derivative wrt xi
    shpj[2][i1] = shp[i1-1][1]*xsjj; // shape function derivative wrt eta
    shpj[3][i1] = shp[i1-1][2]*xsjj; // shape function derivative wrt zeta
    shpj[4][i1] = shp[i1-1][3]*xsj;  // shape function itself
    // printf("\n shp[1][il] %f %f %f %f
\n",shpj[1][i1],shpj[2][i1],shpj[3][i1],shpj[4][i1]);
  } // Loop 4a
  jj1=1;
  for (j=1; j<nope+1; j++) // Loop 4b
  {
    ii1=1;
      for (i=1; i<nope+1; i++) // Loop 5
      {
        for (i1=1; i1<4; i1++) // Loop 6a
        {
      for (j1=1; j1<4; j1++) // Loop 7a
      {
      w[i1][j1] = shpj[i1][i]*shpj[j1][j];

      } // Loop 7a
    }  // Loop 6a
    for (i1=1; i1<4; i1++)  // Loop 6b
      {
        for (j1=1; j1<4; j1++) // Loop 7b
      {
      for (k1=1; k1<4; k1++) // Loop 8
{
        for (l1=1; l1<4; l1++) // Loop 9
{
  s[ii1+i1-1][jj1+j1-1] =
s[ii1+i1-1][jj1+j1-1]+C4[i1][k1][j1][l1]*w[k1][l1]*weight;
  ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"C4 %f
\n",C4[i1][k1][j1][l1]);CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"w %f
\n",w[k1][l1]);CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"weight %f
\n",weight);CHKERRQ(ierr);
}   // Loop 9
}  // Loop 8
      } // Loop 7b
        } // Loop 6b
    ii1 = ii1 + 3;
    } // Loop 5
    jj1 = jj1 + 3;
  } // Loop 4b
  } // Loop 3
  // ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"C4 %f
\n",C4[1][1][1][1]);CHKERRQ(ierr);


for (ii=1; ii<nope*3+1; ii++)
{
for (jj=1; jj<nope*3+1; jj++)
{ ke[(ii-1)*nope*3+(jj-1)] = s[ii][jj];
// printf("\n ke[(%d-1)*nope*3+(%d-1)] = %f \n", ii,jj,
ke[(ii-1)*nope*3+(jj-1)] );
}
}

// ierr = PetscPrintf(PETSC_COMM_SELF,"exiting spaf elemnts tet10
%f\n",ke[3200]);CHKERRQ(ierr);

return 0;

}

Thanks in advance
Kaushik

On Tue, Apr 21, 2020 at 10:51 AM Dave May <dave.mayhem23 at gmail.com> wrote:

>
>
> On Tue 21. Apr 2020 at 16:47, Matthew Knepley <knepley at gmail.com> wrote:
>
>> You are overwriting memory somewhere. The prints just move it around. I
>> suggest running with valgrind.
>>
>
> Matt is right. However, judging by the code snippet I bet all the arrays
> in question are statically allocated, thus valgrind may be of somewhat
> limited use.
>
> If you send the entire function, or all the related pieces of code,
> someone in the list might spot the error.
>
> Thanks
> Dave
>
>
>
>
>>   Thanks,
>>
>>     Matt
>>
>> On Tue, Apr 21, 2020 at 10:44 AM Kaushik Vijaykumar <
>> kaushikv318 at gmail.com> wrote:
>>
>>> Hello group,
>>>
>>> I have been trying to navigate a weird error that I have found in my FEA
>>> code that I am developing using PetSc. The error occurs when, i execute a
>>> call to stiffness generation of a Tetrahederal element. The code returns an
>>> memory error if I don't print the following statements in the function,
>>> (see the ierr print statements below):
>>>
>>> for (i1=1; i1<4; i1++)  // Loop 6b
>>>       {
>>>         for (j1=1; j1<4; j1++) // Loop 7b
>>>       {
>>>       for (k1=1; k1<4; k1++) // Loop 8
>>> {
>>>         for (l1=1; l1<4; l1++) // Loop 9
>>> {
>>>   s[ii1+i1-1][jj1+j1-1] =
>>> s[ii1+i1-1][jj1+j1-1]+C4[i1][k1][j1][l1]*w[k1][l1]*weight;
>>>   ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"C4 %f
>>> \n",C4[i1][k1][j1][l1]);CHKERRQ(ierr);
>>>   ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"w %f
>>> \n",w[k1][l1]);CHKERRQ(ierr);
>>>   ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"weight %f
>>> \n",weight);CHKERRQ(ierr);
>>> }   // Loop 9
>>> }  // Loop 8
>>>       } // Loop 7b
>>>         } // Loop 6b
>>>
>>>
>>> Any help on this is really appreciated.
>>>
>>> Thanks
>>> Kaushik
>>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> <http://www.cse.buffalo.edu/~knepley/>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200421/27b24232/attachment-0001.html>