Блог

Я пытаюсь запустить код, содержащий эту функцию, и получаю сообщение «*** обнаружено разбиение стека ***: завершено».

 void Lapenta_Markidis( long double v[3], long double E[3], long double B[3], long double c2,  long double upart[3] ){
   
   long double   upartk[3], vbar[3];
   long double  tmp[3], Fk[3], C1[3], C2[3] ;
   long double  dupartk[3];
   long double gL = gamma_v(v, c2);
   

    for(int i=0;i<3;i  ){       
        upart[i] = v[i]*gL;       // momentum at start of time step
        upartk[i] = upart[i];
    }


    /* Start of the nonlinear cycle */

    long double     abserr = 1.0;
    long double     tol=1e-14;
    int             nkmax=10;
    int             nk = 0; //


    do{
        long double  J11, J12, J13,J21, J22,  J23, J31, J32, J33, Det;
        long double  iJ11, iJ12, iJ13,iJ21, iJ22,  iJ23, iJ31, iJ32, iJ33;
        long double  gL_new;

        nk     =  nk 1;
        gL_new = gamma_u(upartk, c2);

        for(int i=0;i<3;i  ){  
            vbar[i] = (upart[i]   upartk[i])/(gL_new   gL);
        }
        
        crossP(vbar,B,tmp);

        // Compute residual vector
        for(int i=0;i<3;i  ){ 
            Fk[i] = upartk[i] - upart[i] - q*dt/mp * (E[i]   tmp[i]);
        }

         // Compute auxiliary coefficients
        for(int i=0;i<3;i  ){
            C1[i] = (gL_new   gL - (upartk[i]*(upartk[i]   upart[i])) / (gL_new*c2) )/ pow((gL   gL_new),2) ;
            C2[i] = -( upartk[i] / (gL_new*c2) )/ ((gL_new   gL),2) ;
        }

    
        // Compute Jacobian
          J11 = 1. - (q*dt/mp) * (C2[1] * (upartk[2]   upart[2]) * B[3] - C2[1] * (upartk[3]   upart[3]) * B[2]) ;
          J12 = - (q*dt/mp)*(C1[2] * B[3] - C2[2] * (upartk[3]   upart[3]) * B[2]) ;
          J13 = - (q*dt/mp) * (C2[3] * (upartk[2]   upart[2]) * B[3] - C1[3] * B[2]) ;
          J21 = - q*dt/mp * (- C1[1] * B[3]   C2[1] * (upartk[3]   upart[3]) * B[1]) ;
          J22 = 1. - q*dt/mp * (- C2[2] * (upartk[1]   upart[1]) * B[3]   C2[2] * (upartk[3]   upart[3]) * B[1]) ;
          J23 = - q*dt/mp * (- C2[3] * (upartk[1]   upart[1]) * B[3]   C1[3] * B[1]) ;
          J31 = - q*dt/mp * (C1[1] * B[2] - C2[1] * (upartk[2]   upart[2]) * B[1]) ;
          J32 = - q*dt/mp * (C2[2] * (upartk[1]   upart[1]) * B[2] - C1[2] * B[1]) ;
          J33 = 1. - q*dt/mp * (C2[3] * (upartk[1]   upart[1]) * B[2] - C2[3] * (upartk[2]   upart[2]) * B[1]) ;

          // Compute inverse Jacobian
           Det = J11*J22*J33   J21*J32*J13   J31*J12*J23 - J11*J32*J23 - J31*J22*J13 - J21*J12*J33;
           iJ11 = (J22*J33 - J23*J32) / Det ;
           iJ12 = (J13*J32 - J12*J33) / Det ;
           iJ13 = (J12*J23 - J13*J22) / Det ;
           iJ21 = (J23*J31 - J21*J33) / Det ;
           iJ22 = (J11*J33 - J13*J31) / Det ;
           iJ23 = (J13*J21 - J11*J23) / Det ;
           iJ31 = (J21*J32 - J22*J31) / Det ;
           iJ32 = (J12*J31 - J11*J32) / Det ;
           iJ33 = (J11*J22 - J12*J21) / Det ;

        // Compute new upartk = upartk - J^(-1) * F(upartk)
          dupartk[1] = - (iJ11 * Fk[1]   iJ12 * Fk[2]   iJ13 * Fk[3]);
          dupartk[2] = - (iJ21 * Fk[1]   iJ22 * Fk[2]   iJ23 * Fk[3]);
          dupartk[3] = - (iJ31 * Fk[1]   iJ32 * Fk[2]   iJ33 * Fk[3]);

        // Check convergence
        for(int i=0;i<3;i  ){
          upartk[i]  =  dupartk[i] ;
        }
          abserr = sqrt(dotP(dupartk, dupartk));
    } while(abserr > tol amp;amp; nk < nkmax); // End of while -> end of cycle

    // Update velocity 
    for(int i=0;i<3;i  ){
        upart[i] = upartk[i];
    }
   
    
}
 

Я пытаюсь запустить код, содержащий эту функцию, и получаю сообщение «*** обнаружено разбиение стека ***: завершено».

Есть какие-нибудь предположения о том, что я делаю не так? Я не так хорошо знаком с синтаксисом C, я неправильно объявил переменную, матрицу?

Вы смешиваете индексацию на основе 1 и индексацию на основе 0. Но массивы C используют индексирование на основе 0.

На нескольких позициях , которые вы используете variable[3] , только там, где variable[2] разрешено:

       dupartk[3] = - (iJ31 * Fk[1]   iJ32 * Fk[2]   iJ33 * Fk[3]);
//            ^                                               ^
 

Переместите все эти индексы по одному, и эти доступы будут в порядке:

       dupartk[0] = - (iJ11 * Fk[0]   iJ12 * Fk[1]   iJ13 * Fk[2]);
      dupartk[1] = - (iJ21 * Fk[0]   iJ22 * Fk[1]   iJ23 * Fk[2]);
      dupartk[2] = - (iJ31 * Fk[0]   iJ32 * Fk[1]   iJ33 * Fk[2]);
 

Но имейте в виду, что есть несколько других неверных показателей, например B[3] , C[3] И. Проверьте каждый индекс.

Индексы массива взяты из zero не one .

Вы пишете и читаете за пределами всех своих массивов в своем коде.

Вам просто нужно уменьшить все индексы на 1 .