SolverLoop_DKE_c

C/C++ source

0001 static char help[] = "Solver loop";
0002 #include <stdio.h>
0003 #include <stdlib.h>
0004 #include <math.h>
0005 #include <string.h>
0006 #include "Fokker_Planck_DKE_loop.h"
0007 #include "petscksp.h"
0008 
0009 int main(int argc, char **args){
0010 PetscErrorCode ierr;
0011  ierr=PetscInitialize(&argc,&args,(char *)0,help);CHKERRQ(ierr);
0012 PetscMPIInt rank;
0013 MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
0014 DOUBLE Pi=4.0*atan(1.0);
0015 int it_f,ir;
0016 int i,j,k;
0017 int sizeZmask;
0018 /*                                     */
0019 
0020 /*Intermediate variables*/
0021 struct Tensor A;
0022 int *sdim;
0023 struct Tensor normf0,normf1,curr0,curr1,prec_curr,prec_normf,residu_f,prec_residu_f;
0024 struct Tensor XI;
0025 struct Tensor SXXf0,SXXf02,SXXsinksource2,SXXfM,SXXRlambda_b_p1p1,SXXfinit;
0026 struct Tensor MMXq_t,MMXR_t;
0027 struct Tensor BT,BP,BR;
0028 struct Tensor B,C,D;
0029 double sum1=0.0,sum4=0.0;
0030 int sum2=0,sum3=0;
0031 struct Tensor timescheme;
0032 struct IntTensor *pnr,*psm_f,*pnmhu,*pns_f,*pnpn;
0033 double *xvals;
0034 struct Tensor MMXS_t;
0035 int i0,j0,k0;
0036 /*INPUT FOR LOOP=============================*/
0037 struct Tensor XXsinksource,Xmask_r_edge;
0038 struct Tensor Xmask_r_edge_t; 
0039 struct Tensor Xmhu;
0040 struct Tensor Xpn2;
0041 struct SparseMatrix MMXT_f_t,MMXR_f_t,MMXP_f_t,MMXPC_f_t,MMXPDI_f_t;
0042 struct Tensor MMXf1_t;
0043 struct Tensor Xf1_t;
0044 struct Tensor XXf0,XXI,XXfM,XXRlambda_b_p1p1,XXfinit;
0045 struct IntTensor *Zmask_f0_t;
0046 struct Tensor* Zpn2_t;
0047 struct Tensor dmhu;
0048 struct Tensor dtn;
0049 struct Tensor mhu;
0050 struct IntTensor *pnit_f;
0051 struct IntTensor* pboundary_mode_f;
0052 struct IntTensor* pbounce_mode;
0053 struct IntTensor* pcoll_mode;
0054 struct IntTensor* pnorm_mode_f;
0055 struct Tensor pprec0_f;
0056 struct IntTensor* pdisplay_mode;
0057 char dir_in[1000],dir_out[1000];
0058 char buff[1000];
0059 
0060 
0061 
0062 
0063 /* printf("argc:%i\n",argc);
0064 printf("args:%s\n",args[0]);
0065 printf("args:%s\n",args[1]);
0066 printf("args:%s\n",args[2]);*/
0067 
0068 
0069 strcpy(dir_in,args[1]);
0070 strcpy(dir_out,args[2]);
0071 
0072 /*printf("dir_in = %s\n",dir_in);
0073 printf("dir_out = %s\n",dir_out);*/
0074 
0075 
0076 
0077 
0078 strcpy(buff,dir_in);strcat(buff,"");
0079 
0080 /*printf("READING INPUT VARIABLES\n");*/
0081 
0082 /*fprintf(stdout,"Read nmhu\n");*/strcpy(buff,dir_in);strcat(buff,"nmhu");ReadBinaryIntArray_DKE(buff,&pnmhu);
0083 int nmhu=(*pnmhu).vals[0];
0084 /*fprintf(stdout,"Read npn\n");*/strcpy(buff,dir_in);strcat(buff,"npn");ReadBinaryIntArray_DKE(buff,&pnpn);
0085 int npn=(*pnpn).vals[0];
0086 /*fprintf(stdout,"Read ns_f\n");*/strcpy(buff,dir_in);strcat(buff,"ns_f");ReadBinaryIntArray_DKE(buff,&(pns_f));
0087 int* ns_f=(*pns_f).vals;
0088 /*fprintf(stdout,"Read nr\n");*/strcpy(buff,dir_in);strcat(buff,"nr");ReadBinaryIntArray_DKE(buff,&(pnr));
0089 int nr=(*pnr).vals[0];
0090 /*fprintf(stdout,"Read sm_f\n");*/strcpy(buff,dir_in);strcat(buff,"sm_f");ReadBinaryIntArray_DKE(buff,&(psm_f));
0091 int* sm_f=(*psm_f).vals;
0092 /*fprintf(stdout,"Read time_scheme\n");*/strcpy(buff,dir_in);strcat(buff,"time_scheme");ReadBinary_DKE(buff,&timescheme);
0093 double time_scheme=timescheme.vals[0];
0094 /*fprintf(stdout,"Read dtn\n");*/strcpy(buff,dir_in);strcat(buff,"dtn");ReadBinary_DKE(buff,&(dtn));
0095 /*fprintf(stdout,"Read mhu\n");*/strcpy(buff,dir_in);strcat(buff,"mhu");ReadBinary_DKE(buff,&mhu);
0096 
0097 /*fprintf(stdout,"Read nit_f\n");*/strcpy(buff,dir_in);strcat(buff,"nit_f");ReadBinaryIntArray_DKE(buff,&(pnit_f));
0098 int nit_f=(*pnit_f).vals[0];
0099 /*fprintf(stdout,"Read boundary_mode_f\n");*/strcpy(buff,dir_in);strcat(buff,"boundary_mode_f");ReadBinaryIntArray_DKE(buff,&pboundary_mode_f);
0100 int boundary_mode_f=(*pboundary_mode_f).vals[0];
0101 /*fprintf(stdout,"Read bounce_mode\n");*/strcpy(buff,dir_in);strcat(buff,"bounce_mode");ReadBinaryIntArray_DKE(buff,&pbounce_mode);
0102 int bounce_mode=(*pbounce_mode).vals[0];
0103 /*fprintf(stdout,"Read coll_mode\n");*/strcpy(buff,dir_in);strcat(buff,"coll_mode");ReadBinaryIntArray_DKE(buff,&pcoll_mode);
0104 int coll_mode=(*pcoll_mode).vals[0];
0105 /*fprintf(stdout,"Read norm_mode_f\n");*/strcpy(buff,dir_in);strcat(buff,"norm_mode_f");ReadBinaryIntArray_DKE(buff,&pnorm_mode_f);
0106 int norm_mode_f=(*pnorm_mode_f).vals[0];
0107 /*fprintf(stdout,"Read prec0_f\n");*/strcpy(buff,dir_in);strcat(buff,"prec0_f");ReadBinary_DKE(buff,&pprec0_f);
0108 /*printf("prec0=%i\n",(pprec0_f).ndims);*/
0109 double prec0_f=(pprec0_f).vals[0];
0110 /*fprintf(stdout,"Read display_mode\n");*/strcpy(buff,dir_in);strcat(buff,"display_mode");ReadBinaryIntArray_DKE(buff,&pdisplay_mode);
0111 int display_mode=(*pdisplay_mode).vals[0];
0112 
0113 
0114 /*fprintf(stdout,"Read Xmhu\n");*/strcpy(buff,dir_in);strcat(buff,"Xmhu");ReadBinary_DKE(buff,&Xmhu);
0115 /*fprintf(stdout,"Read XXsinksource\n");*/strcpy(buff,dir_in);strcat(buff,"XXsinksource");ReadBinary_DKE(buff,&XXsinksource);
0116 /*fprintf(stdout,"Read Xmask_r_edge\n");*/strcpy(buff,dir_in);strcat(buff,"Xmask_r_edge");ReadBinary_DKE(buff,&Xmask_r_edge);
0117 /*fprintf(stdout,"Read Xpn2\n");*/strcpy(buff,dir_in);strcat(buff,"Xpn2");ReadBinary_DKE(buff,&Xpn2);
0118 /*fprintf(stdout,"Read XXI\n");*/strcpy(buff,dir_in);strcat(buff,"XXI");ReadBinary_DKE(buff,&XXI);
0119 /*fprintf(stdout,"Read XXfM\n");*/strcpy(buff,dir_in);strcat(buff,"XXfM");ReadBinary_DKE(buff,&XXfM);
0120 /*fprintf(stdout,"Read XXRlambda_b_p1p1\n");*/strcpy(buff,dir_in);strcat(buff,"XXRlambda_b_p1p1");ReadBinary_DKE(buff,&XXRlambda_b_p1p1);
0121 /*fprintf(stdout,"Read XXfinit\n");*/strcpy(buff,dir_in);strcat(buff,"XXfinit");ReadBinary_DKE(buff,&XXfinit);
0122  /*MMXT_f_t,MMXR_f_t,MMXP_f_t  must be read from sparse format*/
0123 /*fprintf(stdout,"Read MMXT_f_t\n");*/strcpy(buff,dir_in);strcat(buff,"MMXT_f_t");ReadBinarySparseMatrix_DKE(buff,&MMXT_f_t);
0124 /*fprintf(stdout,"Read MMXR_f_t\n");*/strcpy(buff,dir_in);strcat(buff,"MMXR_f_t");ReadBinarySparseMatrix_DKE(buff,&MMXR_f_t);
0125 /*fprintf(stdout,"Read MMXP_f_t\n");*/strcpy(buff,dir_in);strcat(buff,"MMXP_f_t");ReadBinarySparseMatrix_DKE(buff,&MMXP_f_t);
0126 /*fprintf(stdout,"ReadMMXPC_f_t\n");*/strcpy(buff,dir_in);strcat(buff,"MMXPC_f_t");ReadBinarySparseMatrix_DKE(buff,&MMXPC_f_t);
0127 /*fprintf(stdout,"ReadMMXPDI_f_t\n");*/strcpy(buff,dir_in);strcat(buff,"MMXPDI_f_t");ReadBinarySparseMatrix_DKE(buff,&MMXPDI_f_t);
0128 
0129 /*fprintf(stdout,"Read Zmask_f0_t\n");*/strcpy(buff,dir_in);strcat(buff,"Zmask_f0_t");ReadBinaryIntArray_DKE(buff,&Zmask_f0_t);
0130 
0131 /*fprintf(stdout,"ReadZpn2_t\n");*/strcpy(buff,dir_in);strcat(buff,"Zpn2_t");ReadBinaryArray_DKE(buff,&Zpn2_t);
0132 
0133 /*fprintf(stdout,"Read dmhu\n");*/strcpy(buff,dir_in);strcat(buff,"dmhu");ReadBinary_DKE(buff,&dmhu);
0134 
0135 /*Need to initialize the following variables*/
0136 /*fprintf(stdout,"Read XXf0\n");*/strcpy(buff,dir_in);strcat(buff,"XXf0");ReadBinary_DKE(buff,&XXf0);
0137 
0138 /*===========================================*/
0139 
0140 /*INPUT FOR firstordercoolop_dke_ft============================================*/
0141 
0142 struct Tensor f1;
0143 struct Tensor betath_ref;
0144 struct Tensor masku,maskl; 
0145 struct Tensor dpn; 
0146 struct Tensor xpn, xpn2,xpn3,xsigma,xgamma;
0147 struct Tensor xz, xz2;
0148 struct Tensor pn,pn2 ,pn3; 
0149 struct Tensor sigma,gamma;
0150 struct Tensor z,z2;
0151 struct Tensor xJ1,xJ2,xJ3,xJ4;
0152 struct Tensor J1,J2,J3,J4;
0153 struct Tensor xTe_norm,xTe_normir,xne_norm;
0154 
0155 /*OUTPUT FOR firstordercollop_dke_ft==============================================*/
0156 
0157 struct Tensor I1,I2;
0158 
0159 /*=============================================================================*/
0160 
0161 /*READING VARIABLES FOR firstordercollop_dke_ft*/
0162 
0163 /*fprintf(stdout,"Readf1\n");*//*ReadBinary_DKE("/Users/Yves/CEA/Projets/LUKE/Packages/petsc-2.3.2-p10/DKE/DKESolver_variables/f1",&f1);*/
0164 /*fprintf(stdout,"Read dpn\n");*/strcpy(buff,dir_in);strcat(buff,"dpn");ReadBinary_DKE(buff,&dpn);
0165 /*fprintf(stdout,"Read xpn\n");*/strcpy(buff,dir_in);strcat(buff,"xpn");ReadBinary_DKE(buff,&xpn);
0166 /*fprintf(stdout,"Read xpn2\n");*/strcpy(buff,dir_in);strcat(buff,"xpn2");ReadBinary_DKE(buff,&xpn2);
0167 /*fprintf(stdout,"Read xpn3\n");*/strcpy(buff,dir_in);strcat(buff,"xpn3");ReadBinary_DKE(buff,&xpn3);
0168 /*fprintf(stdout,"Read masku\n");*/strcpy(buff,dir_in);strcat(buff,"masku");ReadBinary_DKE(buff,&masku);
0169 /*fprintf(stdout,"Read maskl\n");*/strcpy(buff,dir_in);strcat(buff,"maskl");ReadBinary_DKE(buff,&maskl);
0170 /*fprintf(stdout,"Read xgamma\n");*/strcpy(buff,dir_in);strcat(buff,"xgamma");ReadBinary_DKE(buff,&xgamma);
0171 /*fprintf(stdout,"Read xsigma\n");*/strcpy(buff,dir_in);strcat(buff,"xsigma");ReadBinary_DKE(buff,&xsigma);
0172 /*fprintf(stdout,"Read xz\n");*/strcpy(buff,dir_in);strcat(buff,"xz");ReadBinary_DKE(buff,&xz);
0173 /*fprintf(stdout,"Read xz2\n");*/strcpy(buff,dir_in);strcat(buff,"xz2");ReadBinary_DKE(buff,&xz2);
0174 /*fprintf(stdout,"Read pn\n");*/strcpy(buff,dir_in);strcat(buff,"pn");ReadBinary_DKE(buff,&pn);
0175 /*fprintf(stdout,"Read pn2\n");*/strcpy(buff,dir_in);strcat(buff,"pn2");ReadBinary_DKE(buff,&pn2);
0176 /*fprintf(stdout,"Read pn3\n");*/strcpy(buff,dir_in);strcat(buff,"pn3");ReadBinary_DKE(buff,&pn3);
0177 /*fprintf(stdout,"Read gamma\n");*/strcpy(buff,dir_in);strcat(buff,"gamma");ReadBinary_DKE(buff,&gamma);
0178 /*fprintf(stdout,"Read sigma\n");*/strcpy(buff,dir_in);strcat(buff,"sigma");ReadBinary_DKE(buff,&sigma);
0179 /*fprintf(stdout,"Read z\n");*/strcpy(buff,dir_in);strcat(buff,"z");ReadBinary_DKE(buff,&z);
0180 /*fprintf(stdout,"Read z2\n");*/strcpy(buff,dir_in);strcat(buff,"z2");ReadBinary_DKE(buff,&z2);
0181 /*fprintf(stdout,"Read xJ1\n");*/strcpy(buff,dir_in);strcat(buff,"xJ1");ReadBinary_DKE(buff,&xJ1);
0182 /*fprintf(stdout,"Read xJ2\n");*/strcpy(buff,dir_in);strcat(buff,"xJ2");ReadBinary_DKE(buff,&xJ2);
0183 /*fprintf(stdout,"Read xJ3\n");*/strcpy(buff,dir_in);strcat(buff,"xJ3");ReadBinary_DKE(buff,&xJ3);
0184 /*fprintf(stdout,"Read xJ4\n");*/strcpy(buff,dir_in);strcat(buff,"xJ4");ReadBinary_DKE(buff,&xJ4);
0185 /*fprintf(stdout,"Read J1\n");*/strcpy(buff,dir_in);strcat(buff,"J1");ReadBinary_DKE(buff,&J1);
0186 /*fprintf(stdout,"Read J2\n");*/strcpy(buff,dir_in);strcat(buff,"J2");ReadBinary_DKE(buff,&J2);
0187 /*fprintf(stdout,"Read J3\n");*/strcpy(buff,dir_in);strcat(buff,"J3");ReadBinary_DKE(buff,&J3);
0188 /*fprintf(stdout,"Read J4\n");*/strcpy(buff,dir_in);strcat(buff,"J4");ReadBinary_DKE(buff,&J4);
0189 /*fprintf(stdout,"Read xne_norm\n");*/strcpy(buff,dir_in);strcat(buff,"xne_norm");ReadBinary_DKE(buff,&xne_norm);
0190 /*fprintf(stdout,"Read xTe_norm\n");*/strcpy(buff,dir_in);strcat(buff,"xTe_norm");ReadBinary_DKE(buff,&xTe_norm);
0191 /*fprintf(stdout,"Read betath_ref\n");*/strcpy(buff,dir_in);strcat(buff,"betath_ref");ReadBinary_DKE(buff,&betath_ref);
0192 if (rank==0){
0193 printf("END READ\n");
0194 }
0195 /*END READ*/
0196 
0197 sdim=malloc(2*sizeof(int));
0198 /*Intermediate variables*/
0199 struct Tensor XI_t,MMXf0_t,MMXsource_t,MMXIpn2_t,Xf1;
0200 int sizeMMX;
0201 /*                      */
0202 if (rank==0){
0203 printf("INITIALIZING VARIABLES\n");
0204 }
0205 /*Begin Fokker_Planck solver loop*/
0206 
0207 sdim[0]=nr;
0208 sdim[1]=nit_f+1;
0209 Tensor_Create(2,sdim,&normf0);
0210 Tensor_Create(2,sdim,&normf1);
0211 Tensor_Create(2,sdim,&curr0);
0212 Tensor_Create(2,sdim,&curr1);
0213 Tensor_Create(2,sdim,&prec_curr);
0214 Tensor_Create(2,sdim,&prec_normf);
0215 Tensor_Create(2,sdim,&residu_f);
0216 Tensor_Create(2,sdim,&prec_residu_f);
0217 Tensor_Fill(&normf0,0.0);
0218 Tensor_Fill(&normf1,0.0);
0219 Tensor_Fill(&curr0,0.0);
0220 Tensor_Fill(&curr1,0.0);
0221 Tensor_Fill(&prec_curr,0.0);
0222 Tensor_Fill(&prec_normf,0.0);
0223 Tensor_Fill(&residu_f,1.0);
0224 Tensor_Fill(&prec_residu_f,1.0);
0225 
0226 for(ir=0;ir/*printf("XXfinitsize=%i\n",XXfinit.size);*/
0231     for(i=0;i/*printf("dmhudims=%i %i\n",dmhu.sizedim[0],dmhu.sizedim[1]);    */
0235     integral_dke_jd(dmhu,SXXfinit,2,&B);
0236     /*printf("SXXfinit=%i %i\n",SXXfinit.sizedim[0],SXXf0.sizedim[1]);
0237     printf("Bdims=%i %i\n",B.sizedim[0],B.sizedim[1]);*/
0238     /*fprintf(stdout,"Read pn2\n");ReadBinary_DKE("/cgc_data/zineb/data/LUKE/Packages/petsc-2.3.2-p8/DKE/DKESolver_variables/pn2",&pn2);*/
0239     /*printf("pn2dims=%i %i\n",pn2.sizedim[0],pn2.sizedim[1]);*/
0240     product(pn2,B,&B);
0241     /*printf("Bdims=%i %i\n",B.sizedim[0],B.sizedim[1]);*/
0242     Tensor_Transpose(&B);
0243     integral_dke_jd(dpn,B,1,&C);
0244         normf0.vals[ir+0*normf0.sizedim[0]]=2*Pi*C.vals[0];
0245     product(Xmhu,SXXfinit,&B);
0246     integral_dke_jd(dmhu,B,2,&C);
0247     product(C,pn,&C);
0248     product(C,pn2,&C);
0249     division(C,gamma,&C);
0250     Tensor_Transpose(&C);
0251     integral_dke_jd(dpn,C,1,&D);           
0252         curr0.vals[ir+0*curr0.sizedim[0]]=2.0*Pi*D.vals[0]*xne_norm.vals[ir]/normf0.vals[ir+normf0.sizedim[0]*0];
0253 
0254 }
0255 if (rank==0){printf("END INITIALIZATION\n");}
0256 
0257 /*PETSC variables*/
0258 Vec x,b;
0259 PetscInt low,high;
0260 PetscInt* ix;
0261 DOUBLE * y;
0262 Mat PMMXPC_f_t;
0263 KSP ksp;
0264 PC pc;
0265 PetscTruth flg;
0266 PetscViewer fd;
0267 VecScatter vecscat;
0268 Vec xseq;
0269 
0270 /*               */
0271 
0272 /*Initialize PETSc environment*/
0273 /*                            */
0274 
0275 if (rank==0){printf("SETTING UP PETSc OBJECTS\n");}
0276 /*LOAD SYSTEM MATRIX AS PETSC OBJECT MAT*/
0277 strcpy(buff,dir_in);strcat(buff,"PMMXPC_f_t");
0278 ierr=PetscViewerBinaryOpen(PETSC_COMM_WORLD,buff,FILE_MODE_READ,&fd);CHKERRQ(ierr);
0279 ierr = MatLoad(fd,MATAIJ,&PMMXPC_f_t);CHKERRQ(ierr);
0280 if (flg){ierr = PetscViewerDestroy(fd);CHKERRQ(ierr);}
0281 ierr=KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr); 
0282 ierr=KSPSetOperators(ksp,PMMXPC_f_t,PMMXPC_f_t,SAME_PRECONDITIONER);CHKERRQ(ierr);
0283   ierr=KSPSetFromOptions(ksp);CHKERRQ(ierr);
0284   ierr=KSPGetPC(ksp,&pc);CHKERRQ(ierr);
0285   ierr=PCSetFromOptions(pc);CHKERRQ(ierr);
0286 if (rank==0){printf("END SETUP\n");}
0287 
0288 if (rank==0){printf("BEGIN LOOP\n");}
0289 for(it_f=1;it_f<=nit_f;it_f++){
0290   
0291   sum3=0;
0292   sum4=0.0; 
0293   for(i=0;i=prec0_f);
0295   sum4+=(prec_residu_f.vals[i+(it_f-1)*prec_residu_f.sizedim[0]]);
0296   }
0297   /*printf("prec=%g\n",prec0_f);
0298   printf("sum3=%i\n",sum3);printf("sum4=%g\n",sum4);*/
0299     if(sum3>=1)/*convergence condition)*/
0300 {
0301 /*for(i=0;i<residu_f.sizedim[0];i++){
0302  printf("residu=%g ",residu_f.vals[i+(it_f-1)*residu_f.sizedim[0]]);
0303  }*/
0304        if (rank==0){printf("RHS COMPUTATION\n");}
0305        /*fprintf(stdout,"nrvals=%i\n",nr);*/
0306        /*Allocation for MMXf0_t,MMXsource_t,MMXIpn2_t*/
0307        sizeMMX=0;
0308        for(ir=0;ir/*   
0318                                             */
0319                         
0320        i0=0;
0321        j0=0;    
0322        k0=0;                
0323        for(ir=0;irfor(i=0;i/*Size of submatrices SXXf0 and SXXsinksource*/
0327       }
0328       
0329 
0330           /*Computation of XI_t*/
0331       sdim[0]=XXI.sizedim[0];
0332           sdim[1]=sizeZmask;
0333       /*sdim[1]=(Zmask_f0_t[ir]).sizedim[1];*/
0334       Tensor_Create(2,sdim,&XI_t);      
0335       
0336       j=0;
0337       for(k=0;kif(Zmask_f0_t[ir].vals[k]==1){
0339             for(i=0;i/*printf("XI_t done\n");*/
0347        if(boundary_mode_f>0){
0348               for(i=0;ifor(j=0;j/*Computation of MMXf0_t, MMXq_t,MMXIpn2_t)*/
0355            /*printf("SXXf02\n");*/
0356        sdim[0]=XXf0.sizedim[0];
0357        sdim[1]=sizeZmask; /*Zmask_f0_t[ir] are vectors*/
0358        Tensor_Create(2,sdim,&SXXf02);
0359        j=0;
0360        for(k=0;kif(Zmask_f0_t[ir].vals[k]==1){
0362             for(i=0;i/*printf("SXXf02 done\n");*/
0371       
0372 /*Tensor_Save("/cgc_data/zineb/data/LUKE/Packages/petsc-2.3.2-p8/DKE/DKESolver_output/SXXf02",SXXf02);*/
0373        
0374            Tensor_Transpose(&SXXf02);
0375       
0376       for(i=0;i/*Zmask_f0_t[ir] are vectors*/
0386        Tensor_Create(2,sdim,&SXXsinksource2);
0387        j=0;
0388        for(k=0;kif(Zmask_f0_t[ir].vals[k]==1){
0390             for(i=0;ifor(i=0;i/*printf("i=%i,j=%i\n",i,j0);
0400       printf("MMXsize=%i\n", MMXsource_t.size);*/
0401               MMXsource_t.vals[i+j0]=SXXsinksource2.vals[i];
0402       } 
0403       
0404       free(SXXsinksource2.vals);free(SXXsinksource2.sizedim);
0405        
0406        j0+=i;  
0407        
0408       product(XI_t,Zpn2_t[ir],&A);
0409       Tensor_Transpose(&A);    
0410       /*fprintf(stdout,"MMXIpn2_t\n");  */
0411       for(j=0;j/*MMXIpn2_t.vals[ir*A.size+j]=A.vals[j];*/
0413           MMXIpn2_t.vals[j+k0]=A.vals[j];
0414       }
0415       
0416       k0+=j;
0417       
0418       free(A.vals);free(A.sizedim);
0419       free(XI_t.sizedim);free(XI_t.vals);
0420       
0421       
0422        }/*for(ir=0,....)*/    
0423           
0424        /*Scalar_Tensor_product(-1.0,XI_t,&XI_t);*/    /*Weird sign problem for XI_t*/
0425                    
0426        if(1){/*isempty(Xmask_r_edge),No anomalous radial transport*/
0427           /*fprintf(stdout,"================================\n");*/
0428           sizeZmask=0;
0429           for(i=0;i/*Size of submatrices SXXf0 and SXXsinksource*/
0431       }
0432       
0433        /*Computation of XI_t*/
0434       sdim[0]=XXI.sizedim[0];
0435       sdim[1]=sizeZmask;
0436       Tensor_Create(2,sdim,&XI_t);
0437        
0438             j=0;
0439       for(k=0;kif(Zmask_f0_t[nr-1].vals[k]==1){
0441             for(i=0;iif(boundary_mode_f >0){
0450                for(i=0;ifor(j=0;j/*Zmask_f0_t[nr-1] are vectors*/
0461             Tensor_Create(2,sdim,&SXXf02);
0462             j=0;
0463        for(k=0;kif(Zmask_f0_t[nr-1].vals[k]==1){
0465             for(i=0;i/*Tensor_Save("/cgc_data/zineb/data/LUKE/Packages/petsc-2.3.2-p8/DKE/DKESolver_output/SXXf02test3",SXXf02);*/
0476                     
0477             for(i=0;i/*Zmask_f0_t[(nr-1)] are vectors*/
0487            Tensor_Create(2,sdim,&SXXsinksource2);
0488        
0489          j=0;
0490        for(k=0;kif(Zmask_f0_t[nr-1].vals[k]==1){
0492             for(i=0;i/*Tensor_Save("SXXsinksource2test",SXXsinksource2);*/
0500                Tensor_Transpose(&SXXsinksource2);
0501       
0502            for(i=0;i/*Tensor_Save("Zpn2",Zpn2_t[(nr-1)]);*/           
0510            product(XI_t,Zpn2_t[(nr-1)],&A);
0511            Tensor_Transpose(&A);    
0512            for(j=0;j/*MMXIpn2_t.vals[(nr-1)*A.size+j]=A.vals[j];*/
0514           MMXIpn2_t.vals[j+k0]=A.vals[j];
0515            } 
0516            k0+=j;
0517            free(A.vals);free(A.sizedim);
0518 
0519         
0520       /*printf("done");*/
0521 /* printf("MMXf0.size1=%i =========================================\n",MMXf0_t.size);      
0522 printf("MMXsource.size=%i=========================================\n",MMXsource_t.size);
0523   */  
0524 
0525   }/*end if(1)*/
0526        
0527        else {
0528           fprintf(stdout,"else\n");
0529        
0530           /*Declaration and allocation of Xmask_r_edge_t*/      
0531           sdim[0]=Xmask_r_edge.sizedim[0];
0532           sdim[1]==Zmask_f0_t[nr-1].size;
0533           Tensor_Create(2,sdim,&Xmask_r_edge_t);
0534        
0535            sdim[0]=XXf0.sizedim[0];
0536        sdim[1]=Zmask_f0_t[nr-1].size; /*Zmask_f0_t[nr-1] are vectors*/
0537        Tensor_Create(2,sdim,&SXXf02);
0538        for(i=0;ifor(j=0;jfor(i=0;i/*Zmask_f0_t[(nr-1)] are vectors*/
0555        Tensor_Create(2,sdim,&SXXsinksource2);
0556        
0557        for(i=0;ifor(j=0;jfor(i=0;i/*fprintf(stdout,"MMXIpn2_t\n");*/  
0574        for(j=0;jfor(i=0;ifor(j=0;j/*SAVES===================================================================*/
0589 
0590 /*Tensor_Save("/cgc_data/zineb/data/LUKE/Packages/petsc-2.3.2-p8/DKE/DKESolver_output/XI_t_test",XI_t);
0591 Tensor_Save("/cgc_data/zineb/data/LUKE/Packages/petsc-2.3.2-p8/DKE/DKESolver_output/MMXf0_t_test",MMXf0_t);
0592 Tensor_Save("/cgc_data/zineb/data/LUKE/Packages/petsc-2.3.2-p8/DKE/DKESolver_output/MMXsource_t_test",MMXsource_t);
0593 Tensor_Save("/cgc_data/zineb/data/LUKE/Packages/petsc-2.3.2-p8/DKE/DKESolver_output/MMXIpn2_t_test",MMXIpn2_t);*/
0594 
0595  /*Tensor_Save("XI_t_test",XI_t);
0596  Tensor_Save("MMXf0_t_test",MMXf0_t);
0597  Tensor_Save("MMXsource_t_test",MMXsource_t);
0598  Tensor_Save("MMXIpn2_t_test",MMXIpn2_t);*/
0599  /*=========================================================================*/    
0600      
0601      /*printf("Multiplication");*/
0602     
0603   /*Computation of MMXq_t*/
0604     SparseMat_Tensorproduct(MMXT_f_t,MMXf0_t,&BT);
0605     SparseMat_Tensorproduct(MMXP_f_t,MMXf0_t,&BP);
0606     SparseMat_Tensorproduct(MMXR_f_t,MMXf0_t,&BR);
0607     Tensor_Sum(BT,BP,-(time_scheme-1.0)/time_scheme,&MMXq_t);
0608     Tensor_Sum(MMXq_t,BR,-(time_scheme-1.0)/time_scheme,&MMXq_t);
0609     Tensor_Sum(MMXIpn2_t,MMXq_t,1.0,&MMXq_t);
0610     Tensor_Sum(MMXsource_t,MMXq_t,1.0,&MMXq_t);
0611  /*Computation of MMXR_t*/    
0612         SparseMat_Tensorproduct(MMXPDI_f_t,MMXq_t,&BT); 
0613     SparseMat_Tensorproduct(MMXPC_f_t,MMXf0_t,&MMXR_t); 
0614     Tensor_Sum(BT,MMXR_t,-1.0,&MMXR_t);  
0615     
0616         free(BT.sizedim);free(BT.vals);
0617     free(BP.sizedim);free(BP.vals);
0618     free(BR.sizedim);free(BR.vals);
0619     
0620     if (rank==0){printf("RHS COMPUTED\n");}
0621 
0622 /*Tensor_Save("/cgc_data/zineb/data/LUKE/Packages/petsc-2.3.2-p8/DKE/DKESolver_output/MMXq_tc",MMXq_t);
0623 Tensor_Save("/cgc_data/zineb/data/LUKE/Packages/petsc-2.3.2-p8/DKE/DKESolver_output/MMXR_tc",MMXR_t);*/
0624     
0625   /*Solving step.Computation of MMS_t*/
0626   if (rank==0){printf("CREATING PETSc PARALLEL OBJECTS\n");}
0627   ierr=VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
0628   ierr=VecSetSizes(b,PETSC_DECIDE,MMXR_t.size);CHKERRQ(ierr);
0629   ierr=VecSetType(b,VECMPI);CHKERRQ(ierr);
0630   ierr=VecGetOwnershipRange(b,&low,&high);CHKERRQ(ierr);
0631   ix=malloc((high-low)*sizeof(int));
0632   y=malloc((high-low)*sizeof(DOUBLE));
0633   for(i=low;iif (rank==0){printf("SOLVING\n");}
0642   ierr=KSPSolve(ksp,b,x);CHKERRQ(ierr);
0643   
0644   if (rank==0){printf("RETRIEVING SOLUTION\n");}
0645   VecScatterCreateToAll(x,&vecscat,&xseq);
0646   VecScatterBegin(x,xseq,INSERT_VALUES,SCATTER_FORWARD,vecscat);
0647   VecScatterEnd(x,xseq,INSERT_VALUES,SCATTER_FORWARD,vecscat);
0648   VecScatterDestroy(vecscat);
0649   VecGetArray(xseq, &xvals);
0650   sdim[0]=1;
0651   sdim[1]=MMXR_t.size;
0652   Tensor_Create(2,sdim,&MMXS_t);
0653   for(i=0;iif(it_f==1){
0661 strcpy(buff,dir_out);strcat(buff,"MMXS_tLU");Tensor_Save(buff,MMXS_t);
0662 /*Tensor_Save("/cgc_data/zineb/data/LUKE/Packages/petsc-2.3.2-p8/DKE/DKESolver_output/MMXS_tLU",MMXS_t);*/
0663 } 
0664   
0665   
0666   /*Computation of MMXf1_t=MMXf0_t+MMXS_t*/
0667          Tensor_Sum(MMXf0_t,MMXS_t,1.0,&MMXf1_t);
0668      free(MMXS_t.vals);free(MMXS_t.sizedim);
0669      
0670      /*printf("COMPUTING Xf1_t\n");*/
0671          for(ir=0;ir/*Computation of Xf1_t=F(MMXf1_t,sm_f,nmhu,ns_f_npn)*/
0673        sdim[0]=nmhu-ns_f[ir];
0674        sdim[1]=npn;
0675        Tensor_Create(2,sdim,&Xf1_t);
0676        
0677        for(i=0;i/*Computation of Xf1=F(Xf1_t, nmhu, ns_f, ir)*/
0683      /*printf("COMPUTING Xf1\n");*/
0684 
0685        if(bounce_mode==0){
0686           Tensor_Create(Xf1_t.ndims,Xf1_t.sizedim,&Xf1);
0687           for(i=0;ielse if(bounce_mode==1){
0692           sdim[0]=npn;
0693           sdim[1]=nmhu;
0694           Tensor_Create(2,sdim,&Xf1);
0695           for(i=0;ifor(i=0;ifor(j=0;jfor(j=nmhu/2;j/*for(j=nmhu/2-ns_f[ir];j<nmhu/2;j++){    */                         
0706                   for(j=0;j/*printf("COMPUTING NORMS,CURRENTS, ACCURACIES\n");*/
0716 
0717    /*Computation of normf0,normf1,curr0, curr1, prec_curr, prec_norm,
0718    residu_f_prec_residu_f*/
0719     sdim[0]=XXf0.sizedim[0];
0720     sdim[1]=XXf0.sizedim[1];
0721     Tensor_Create(2,sdim,&SXXf0);
0722     for(i=0;i/*printf("SXXF0=%i %i\n",SXXf0.sizedim[0],SXXf0.sizedim[1]);
0729     printf("Bdims=%i %i\n",B.sizedim[0],B.sizedim[1]);*/
0730     /*fprintf(stdout,"Read pn2\n");ReadBinary_DKE("/cgc_data/zineb/data/LUKE/Packages/petsc-2.3.2-p8/DKE/DKESolver_variables/pn2",&pn2);*/
0731     product(pn2,B,&B);
0732     Tensor_Transpose(&B);
0733     integral_dke_jd(dpn,B,1,&C);    
0734     normf0.vals[ir+normf0.sizedim[0]*it_f]=2.0*Pi*C.vals[0];
0735     free(B.vals);free(B.sizedim);
0736     free(C.vals);free(C.sizedim);
0737     
0738     integral_dke_jd(dmhu,Xf1,2,&B);    
0739     product(pn2,B,&B);
0740     Tensor_Transpose(&B);
0741     integral_dke_jd(dpn,B,1,&C);    
0742     normf1.vals[ir+normf1.sizedim[0]*it_f]=2.0*Pi*C.vals[0];    
0743     free(B.vals);free(B.sizedim);
0744     free(C.vals);free(C.sizedim);
0745 
0746 
0747     product(Xmhu,SXXf0,&B);
0748     integral_dke_jd(dmhu,B,2,&C);
0749     product(C,pn,&C);
0750     product(C,pn2,&C);
0751     division(C,gamma,&C);
0752     Tensor_Transpose(&C);
0753     integral_dke_jd(dpn,C,1,&D);          
0754        curr0.vals[ir+curr0.sizedim[0]*it_f]=2.0*Pi*D.vals[0]*xne_norm.vals[ir]/normf0.vals[ir+normf1.sizedim[0]*it_f];
0755 
0756     free(B.vals);free(B.sizedim);
0757     free(C.vals);free(C.sizedim);
0758     free(D.vals);free(D.sizedim);
0759     
0760     product(Xmhu,Xf1,&B);
0761     integral_dke_jd(dmhu,B,2,&C);
0762     product(C,pn,&C);
0763     product(C,pn2,&C);
0764     division(C,gamma,&C);
0765     Tensor_Transpose(&C);
0766     integral_dke_jd(dpn,C,1,&D);    
0767         curr1.vals[ir+curr1.sizedim[0]*it_f]=2.0*Pi*D.vals[0]*xne_norm.vals[ir]/normf1.vals[ir+normf1.sizedim[0]*it_f];
0768 
0769     free(B.vals);free(B.sizedim);
0770     free(C.vals);free(C.sizedim);
0771     free(D.vals);free(D.sizedim);
0772            
0773         prec_curr.vals[ir+prec_curr.sizedim[0]*it_f]=(curr1.vals[ir+curr1.sizedim[0]*it_f]-curr0.vals[ir+curr0.sizedim[0]*it_f])/curr0.vals[ir+curr0.sizedim[0]*it_f];
0774 
0775 
0776         prec_normf.vals[ir+prec_normf.sizedim[0]*it_f]=(normf1.vals[ir+normf1.sizedim[0]*it_f]-normf0.vals[ir+normf0.sizedim[0]*1])/normf0.vals[ir+0*normf0.sizedim[0]];
0777     
0778     Tensor_Sum(Xf1,SXXf0,-1.0,&B);
0779     product(B,B,&B);
0780     product(B,SXXf0,&B);
0781     integral_dke_jd(dmhu,B,2,&C);    
0782     product(pn2,C,&C);
0783     Tensor_Transpose(&C);
0784     integral_dke_jd(dpn,C,1,&D);      
0785         residu_f.vals[ir+residu_f.sizedim[0]*it_f]=sqrt(2.0*Pi*D.vals[0]/normf1.vals[ir+normf1.sizedim[0]*it_f]/dtn.vals[0]/dtn.vals[0]);
0786         
0787     free(B.vals);free(B.sizedim);
0788     free(C.vals);free(C.sizedim);
0789     free(D.vals);free(D.sizedim);
0790       
0791         prec_residu_f.vals[ir+prec_residu_f.sizedim[0]*it_f]=absolute((residu_f.vals[ir+residu_f.sizedim[0]*it_f]-residu_f.vals[ir+residu_f.sizedim[0]*(it_f-1)])/residu_f.vals[ir+residu_f.sizedim[0]*it_f]);
0792 
0793         
0794 
0795         if(coll_mode==2){
0796        /*Computation of f1,I1,I2,XI=F(I1,I2,f1,ir,XXFM,XXRlambda_b_p1p1,gamma,
0797    mhu,pn,pn2),XXI=F(X1)*/
0798    /*printf("COMPUTING f1,I1,I2,XI,XXI\n");*/
0799            sdim[0]=npn;
0800            sdim[1]=mhu.sizedim[1];
0801            Tensor_Create(2,sdim,&B);
0802        for(i=0;ifor(j=0;jfor(i=0;ifor(i=0;ifor(i=0;ifor(i=0;ifor(j=0;jfor(i=0;ifor(j=0;jfor(i=0;ifor(j=0;j/*printf("%i!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n",it_f);*/
0866           free(f1.vals);free(f1.sizedim);      
0867           free(B.vals);free(B.sizedim);
0868           free(I1.vals);free(I1.sizedim);
0869           free(I2.vals);free(I2.sizedim);
0870           free(XI.vals);free(XI.sizedim);
0871           free(SXXfM.vals);free(SXXfM.sizedim);
0872           free(SXXRlambda_b_p1p1.vals);free(SXXRlambda_b_p1p1.sizedim);
0873           /*free(C.vals);free(C.sizedim);*/
0874           
0875     }
0876        
0877     else{
0878        for(i=0;ifor(j=0;j/*printf("COMPUTING XXf0\n");*/
0885     if(norm_mode_f==1){
0886        /*Computation of XXfo=F(X1,xne_norm(ir),normf1*/
0887        
0888       for(i=0;ifor(j=0;jelse{
0897        
0898        for(i=0;ifor(j=0;j/*printf("%g\n",XXf0.vals[i+j*XXf0.sizedim[0]+ir*XXf0.sizedim[0]*XXf0.sizedim[1]]);*/  }    
0904 /*printf("%g\n",Xf1.vals[i+j*XXf0.sizedim[0]]);*/
0905        }                 
0906     }    
0907 
0908    /*Computation of XXfo=F(X1,xne_norm(ir),normf1*/
0909    /*Computation of normf0,curr0*/
0910    
0911       normf0.vals[ir+normf0.sizedim[0]*it_f]=normf1.vals[ir+normf1.sizedim[0]*it_f];
0912       curr0.vals[ir+curr0.sizedim[0]*it_f]=curr1.vals[ir+curr1.sizedim[0]*it_f];
0913        
0914       free(Xf1_t.sizedim);free(Xf1_t.vals);free(Xf1.sizedim);free(Xf1.vals);    
0915         }
0916     
0917   free(MMXf0_t.vals); free(MMXf0_t.sizedim);    
0918   free(MMXsource_t.vals); free(MMXsource_t.sizedim);      
0919   free(MMXIpn2_t.vals);free(MMXIpn2_t.sizedim);
0920   
0921  }
0922     else{
0923     sum1=0.0;
0924        for(i=0;iif(sum1if (rank==0){printf("WARNING: stationnary residu but convergence not achieved for f0 after %i iterations\n",it_f);}
0928        }
0929        sum2=0;
0930        for(i=0;i=prec0_f);
0932        }
0933         if(sum2==0 && display_mode>0){
0934       if (rank==0){printf("Successfull convergence achieved for f0 after %i iterations\n",it_f-2);}
0935     }  
0936     break;     
0937     }
0938     
0939     if(it_f==nit_f){
0940       if (rank==0){printf("WARNING: convergence not achieved for f0 after %i iterations.Time step too small",nit_f);}
0941     }
0942    
0943 strcpy(buff,dir_out);strcat(buff,"XXf0");Tensor_Save(buff,XXf0);
0944 strcpy(buff,dir_out);strcat(buff,"residu_f1");Tensor_Save(buff,residu_f);
0945 
0946 
0947 
0948 
0949  }    
0950  
0951 /*free(normf0.sizedim);free(normf0.vals);
0952 free(normf1.sizedim);free(normf1.vals);
0953 free(curr0.sizedim);free(curr0.vals);
0954 free(curr1.sizedim);free(curr1.vals);
0955 free(prec_curr.sizedim);free(prec_curr.vals);
0956 free(prec_normf.sizedim);free(prec_normf.vals);
0957 free(residu_f.sizedim);free(residu_f.vals);
0958 free(prec_residu_f.sizedim);free(prec_residu_f.vals);*/
0959 
0960 /*FREE AND SAVES*/
0961 printf("FREEING MEMORY\n");
0962 KSPDestroy(ksp);
0963 
0964 /*free(XXI.vals);free(XXf0.vals);
0965 free(XXI.sizedim);free(XXf0.sizedim);
0966 free(XXsinksource.vals);
0967 free(dpn.vals);
0968  free(xpn.vals);free(xpn2.vals);free(xpn3.vals);
0969  free(masku.vals);free(maskl.vals);
0970  free(xsigma.vals);free(xgamma.vals);
0971  free(xz.vals);free(xz2.vals);
0972  free(z.vals);free(z2.vals);
0973  free(pn.vals);free(pn2.vals);free(pn3.vals);
0974  free(sigma.vals);free(gamma.vals);
0975  free(xJ1.vals);free(xJ2.vals);free(xJ3.vals);free(xJ4.vals);
0976  free(J1.vals);free(J2.vals);free(J3.vals);free(J4.vals);
0977  free(xTe_norm.vals);free(betath_ref.vals);
0978 free(XXsinksource.sizedim);
0979 free(dpn.sizedim);
0980  free(xpn.sizedim);free(xpn2.sizedim);free(xpn3.sizedim);
0981  free(masku.sizedim);free(maskl.sizedim);
0982  free(xsigma.sizedim);free(xgamma.sizedim);
0983  free(xz.sizedim);free(xz2.sizedim);
0984  free(z.sizedim);free(z2.sizedim);
0985  free(pn.sizedim);free(pn2.sizedim);free(pn3.sizedim);
0986  free(sigma.sizedim);free(gamma.sizedim);
0987  free(xJ1.sizedim);free(xJ2.sizedim);free(xJ3.sizedim);free(xJ4.sizedim);
0988  free(J1.sizedim);free(J2.sizedim);free(J3.sizedim);free(J4.sizedim);
0989  free(xTe_norm.sizedim);free(betath_ref.sizedim);*/
0990  
0991 
0992 /*free(index.vals);free(index.sizedim);*/
0993 /*Tensor_Save("curr0",curr0);
0994 Tensor_Save("curr1",curr1);
0995 Tensor_Save("normf0",normf0);
0996 Tensor_Save("normf1",normf1);
0997 Tensor_Save("prec_normf",prec_normf);
0998 Tensor_Save("prec_curr",prec_curr);
0999 Tensor_Save("residu_f",residu_f);
1000 Tensor_Save("prec_residu_f",prec_residu_f);*/
1001 
1002 ierr=PetscFinalize();CHKERRQ(ierr);
1003 
1004 }
1005