00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 static char const Id[] = "$Id: hdemag.c 3013 2010-03-26 16:12:31Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #include "field/field.h"
00036 #include "init/init.h"
00037 #include "util/util.h"
00038
00039 static int doinit=1;
00040 static int fieldon=0;
00041
00042 static Mat Bmat;
00043 static Mat AdivM;
00044 static Mat Agradu;
00045 static Mat ABu1;
00046 static Mat ABu2;
00047 static Vec u2bnd;
00048 static Vec u1;
00049 static Vec u1bnd;
00050 static Vec u2;
00051 static Vec bu1;
00052 static Vec bu2;
00053 static Vec utot;
00055 KSP ksp_Au1;
00056 KSP ksp_Au2;
00057
00058 #ifdef HDEMAGUEPOL
00059 static Vec u1bak;
00060 static Vec u2bak;
00061 static Vec du1dt;
00062 static Vec du2dt;
00063 #endif
00064
00065 #include "bmatrix.c"
00066
00067 int Hdemag_Init(GridData *gdata)
00068 {
00069 MagparFunctionLogBegin;
00070
00071 int rank,size;
00072 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00073 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00074
00075 gdata->VHdem=PETSC_NULL;
00076
00077 PetscTruth flg;
00078 ierr = PetscOptionsGetInt(PETSC_NULL,"-demag",(PetscInt*)&fieldon,&flg);CHKERRQ(ierr);
00079 if (flg!=PETSC_TRUE) {
00080 fieldon=1;
00081 PetscPrintf(PETSC_COMM_WORLD,
00082 "Option -demag not found, using default value: %i\n",
00083 fieldon
00084 );
00085 }
00086
00087 if (!fieldon) {
00088 PetscPrintf(PETSC_COMM_WORLD,"Hdemag off\n");
00089 MagparFunctionLogReturn(0);
00090 }
00091 PetscPrintf(PETSC_COMM_WORLD,"Hdemag on\n");
00092
00093
00094
00095
00096 #if 0
00097 #define U1DIRICHLET 0
00098 PetscPrintf(PETSC_COMM_WORLD,"Using 1 Dirichlet BC for u1 on node %i\n",U1DIRICHLET);
00099 #else
00100 PetscPrintf(PETSC_COMM_WORLD,"No Dirichlet BC for u1 (floating)\n");
00101 #endif
00102
00103
00104
00105
00106
00107
00108 int ln_bmatrixrows;
00109 ln_bmatrixrows=gdata->n_vert_bnd/size;
00110 if ((gdata->n_vert_bnd%size)/(rank+1)>=1) ln_bmatrixrows++;
00111
00112 int ln_vert;
00113 ln_vert=gdata->ln_vert;
00114 int n_vert;
00115 n_vert=gdata->n_vert;
00116
00117 PetscPrintf(PETSC_COMM_WORLD,"Calculating sparsity pattern for matrix preallocation\n");
00118
00119
00120 int low,high;
00121 ierr = VecGetOwnershipRange(gdata->elevol,&low,&high);CHKERRQ(ierr);
00122
00123 PetscSynchronizedPrintf(PETSC_COMM_WORLD,
00124 "<%i>element ownership: %i - %i\n",
00125 rank,low,high-1
00126 );
00127 PetscSynchronizedFlush(PETSC_COMM_WORLD);
00128
00129 int *elevertall;
00130
00131 #ifdef HAVE_ELEVERTALL
00132 elevertall=gdata->elevertall;
00133 #else
00134 if (size>1) {
00135 ierr = PetscMalloc(NV*gdata->n_ele*sizeof(int),&elevertall);CHKERRQ(ierr);
00136 ierr = PetscMemcpy(elevertall+NV*low,gdata->elevert,NV*gdata->ln_ele*sizeof(int));CHKERRQ(ierr);
00137
00138 int *p;
00139 p=elevertall;
00140
00141 int recvcount[size];
00142 recvcount[rank]=NV*gdata->ln_ele;
00143 for (int i=0;i<size;i++) {
00144 ierr = MPI_Bcast(recvcount+i,1,MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
00145
00146 ierr = MPI_Bcast(p,recvcount[i],MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
00147 p+=recvcount[i];
00148 }
00149 assert(p==elevertall+NV*gdata->n_ele);
00150 }
00151 else {
00152 elevertall=gdata->elevert;
00153 }
00154 #endif
00155
00156 int *ia,*ja;
00157 #ifdef HAVE_M2IJ
00158 ia=gdata->mesh2nodal_ia;
00159 ja=gdata->mesh2nodal_ja;
00160 #else
00161 ierr = Mesh2Nodal(gdata->n_ele,gdata->n_vert,elevertall,&ia,&ja);CHKERRQ(ierr);
00162 #endif
00163
00164
00165 ierr = VecGetOwnershipRange(gdata->M,&low,&high);CHKERRQ(ierr);
00166 low/=ND;
00167 high/=ND;
00168
00169
00170 ierr = VecCreate(PETSC_COMM_WORLD,&utot);CHKERRQ(ierr);
00171 ierr = VecSetSizes(utot,ln_vert,n_vert);CHKERRQ(ierr);
00172 ierr = VecSetFromOptions(utot);CHKERRQ(ierr);
00173
00174 #define XND 1
00175 int *d_nnz;
00176 int *o_nnz;
00177 ierr = PetscMalloc(ND*gdata->n_vert*sizeof(int),&d_nnz);CHKERRQ(ierr);
00178 ierr = PetscMalloc(ND*gdata->n_vert*sizeof(int),&o_nnz);CHKERRQ(ierr);
00179
00180 for (int i=0;i<gdata->n_vert;i++) {
00181 for (int j=0;j<XND;j++) {
00182 d_nnz[XND*i+j]=XND;
00183 o_nnz[XND*i+j]=0;
00184 for (int l=ia[i];l<ia[i+1];l++) {
00185 if (ja[l]>=low && ja[l]<high) {
00186 d_nnz[XND*i+j]+=1;
00187 }
00188 else {
00189 o_nnz[XND*i+j]+=1;
00190 }
00191 }
00192 }
00193 }
00194
00195 int *d_nnzm;
00196 int *o_nnzm;
00197 ierr = PetscMalloc(ND*n_vert*sizeof(int),&d_nnzm);CHKERRQ(ierr);
00198 ierr = PetscMalloc(ND*n_vert*sizeof(int),&o_nnzm);CHKERRQ(ierr);
00199
00200
00201 int lowm,highm;
00202
00203 lowm=low;
00204 highm=high;
00205
00206 ierr = PetscMemcpy(d_nnzm,d_nnz,ND*n_vert*sizeof(int));CHKERRQ(ierr);
00207 ierr = PetscMemcpy(o_nnzm,o_nnz,ND*n_vert*sizeof(int));CHKERRQ(ierr);
00208
00209
00210 Mat Au1;
00211 ierr = MatCreateMPIAIJ(
00212 PETSC_COMM_WORLD,
00213 ln_vert,ln_vert,
00214 n_vert,n_vert,
00215 0,d_nnzm+lowm,0,o_nnzm+lowm,
00216 &Au1
00217 );CHKERRQ(ierr);
00218 ierr = MatSetFromOptions(Au1);CHKERRQ(ierr);
00219
00220 ierr = MatCreateMPIAIJ(
00221 PETSC_COMM_WORLD,
00222 ln_bmatrixrows,ln_vert,
00223 gdata->n_vert_bnd,n_vert,
00224 1,PETSC_NULL,1,PETSC_NULL,
00225 &ABu1
00226 );CHKERRQ(ierr);
00227 ierr = MatSetFromOptions(ABu1);CHKERRQ(ierr);
00228
00229 Mat Au2;
00230 ierr = MatCreateMPIAIJ(
00231 PETSC_COMM_WORLD,
00232 ln_vert,ln_vert,
00233 n_vert,n_vert,
00234 0,d_nnzm+lowm,0,o_nnzm+lowm,
00235 &Au2
00236 );CHKERRQ(ierr);
00237 ierr = MatSetFromOptions(Au2);CHKERRQ(ierr);
00238
00239
00240
00241 for (int i=lowm; i<highm; i++) {
00242 d_nnzm[i] = PetscMin(d_nnzm[i],ln_bmatrixrows);
00243 o_nnzm[i] = PetscMin(d_nnzm[i],gdata->n_vert_bnd-ln_bmatrixrows);
00244 assert(d_nnzm[i]>=0);
00245 assert(o_nnzm[i]>=0);
00246 }
00247
00248 ierr = MatCreateMPIAIJ(
00249 PETSC_COMM_WORLD,
00250 ln_vert,ln_bmatrixrows,
00251 n_vert,gdata->n_vert_bnd,
00252 0,d_nnzm+lowm,0,o_nnzm+lowm,
00253 &ABu2
00254 );CHKERRQ(ierr);
00255 ierr = MatSetFromOptions(ABu2);CHKERRQ(ierr);
00256
00257
00258 for (int i=0; i<n_vert; i++) {
00259 d_nnzm[i] = 0;
00260 o_nnzm[i] = 0;
00261 }
00262 for (int i=0; i<gdata->n_vert; i++) {
00263 int i2;
00264 i2=i;
00265 if (i<lowm || i>=highm) continue;
00266 assert(d_nnz[i]>=0);
00267 assert(o_nnz[i]>=0);
00268 d_nnzm[i2] += ND*d_nnz[i]+ND*o_nnz[i];
00269 o_nnzm[i2] += ND*d_nnz[i]+ND*o_nnz[i];
00270 assert(i2<n_vert);
00271 assert(d_nnzm[i2]>=0);
00272 assert(o_nnzm[i2]>=0);
00273 }
00274
00275
00276 ierr = MatCreateMPIAIJ(
00277 PETSC_COMM_WORLD,
00278 ln_vert,ND*gdata->ln_vert,
00279 n_vert,ND*gdata->n_vert,
00280 0,d_nnzm+lowm,0,o_nnzm+lowm,
00281 &AdivM
00282 );CHKERRQ(ierr);
00283 ierr = MatSetFromOptions(AdivM);CHKERRQ(ierr);
00284
00285 for (int i=ND*gdata->n_vert-1; i>=0; i--) {
00286 d_nnz[i] = d_nnz[i/ND]+o_nnz[i/ND];
00287 o_nnz[i] = d_nnz[i/ND]+o_nnz[i/ND];
00288 assert(d_nnz[i]>=0);
00289 assert(o_nnz[i]>=0);
00290 }
00291
00292
00293 ierr = MatCreateMPIAIJ(
00294 PETSC_COMM_WORLD,
00295 ND*gdata->ln_vert,ln_vert,
00296 ND*gdata->n_vert,n_vert,
00297 0,d_nnz+ND*low,0,o_nnz+ND*low,
00298 &Agradu
00299 );CHKERRQ(ierr);
00300 ierr = MatSetFromOptions(Agradu);CHKERRQ(ierr);
00301
00302 #ifndef HAVE_M2IJ
00303 ierr = PetscFree(ia);CHKERRQ(ierr);
00304 ierr = PetscFree(ja);CHKERRQ(ierr);
00305 #endif
00306 ierr = PetscFree(d_nnz);CHKERRQ(ierr);
00307 ierr = PetscFree(o_nnz);CHKERRQ(ierr);
00308 #ifndef HAVE_ELEVERTALL
00309 if (size>1) {
00310 ierr = PetscFree(elevertall);CHKERRQ(ierr);
00311 }
00312 #endif
00313 ierr = PetscFree(d_nnzm);CHKERRQ(ierr);
00314 ierr = PetscFree(o_nnzm);CHKERRQ(ierr);
00315
00316
00317 PetscReal *ta_elevol;
00318 ierr = VecGetArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);
00319
00320 PetscPrintf(PETSC_COMM_WORLD,"Calculating stiffness matrices...\n");
00321 for (int i=0; i<gdata->ln_ele; i++) {
00322 ierr = ProgressBar(i,gdata->ln_ele,10);
00323
00324
00325 int t_v[NV];
00326 int t_w[NV];
00327
00328 for (int j=0; j<NV; j++) {
00329 t_v[j]=gdata->elevert[i*NV+j];
00330 t_w[j]=t_v[j];
00331 }
00332
00333
00334 PetscReal D_etaj[NV][ND];
00335 tetgrad(
00336 gdata->vertxyz+t_v[0]*ND,
00337 gdata->vertxyz+t_v[1]*ND,
00338 gdata->vertxyz+t_v[2]*ND,
00339 gdata->vertxyz+t_v[3]*ND,
00340 D_etaj
00341 );
00342
00343
00344
00345 for (int l=0;l<NV;l++) {
00346 for (int j=0;j<NV;j++) {
00347
00348 PetscReal matele;
00349 matele=-ta_elevol[i]*my_ddot(ND,D_etaj[j],1,D_etaj[l],1);
00350
00351
00352
00353
00354 if (PetscAbsReal(matele)>D_EPS) {
00355 #ifdef U1DIRICHLET
00356
00357
00358
00359 if (t_w[j] != U1DIRICHLET &&
00360 t_w[l] != U1DIRICHLET) {
00361 #endif
00362 ierr = MatSetValue(
00363 Au1,
00364 t_w[j],
00365 t_w[l],
00366 matele,
00367 ADD_VALUES
00368 );CHKERRQ(ierr);
00369
00370 assert(t_w[j]<n_vert);
00371 assert(t_w[l]<n_vert);
00372 #ifdef U1DIRICHLET
00373 }
00374 #endif
00375
00376 int tvj,tvl;
00377 tvj=t_v[j];
00378 tvl=t_v[l];
00379 if (gdata->vertbndg2bnd[tvj]==C_INT) {
00380
00381 if (gdata->vertbndg2bnd[tvl]==C_INT) {
00382
00383
00384
00385 ierr = MatSetValue(
00386 Au2,
00387 t_w[j],
00388 t_w[l],
00389 matele,
00390 ADD_VALUES
00391 );CHKERRQ(ierr);
00392 assert(t_w[j]>=0);
00393 assert(t_w[l]>=0);
00394 assert(t_w[j]<n_vert);
00395 assert(t_w[l]<n_vert);
00396 }
00397 else {
00398
00399 ierr = MatSetValue(
00400 ABu2,
00401 t_w[j],
00402 gdata->vertbndg2bnd[tvl],
00403 -matele,
00404 ADD_VALUES
00405 );CHKERRQ(ierr);
00406 assert(t_w[j]<n_vert);
00407 assert(gdata->vertbndg2bnd[tvl]<gdata->n_vert_bnd);
00408 }
00409 }
00410 }
00411 }
00412 }
00413
00414 for (int l=0;l<NV;l++) {
00415 for (int k=0;k<ND;k++) {
00416 PetscReal matele;
00417 matele=ta_elevol[i]/NV*D_etaj[l][k];
00418 if (PetscAbsReal(matele)>D_EPS) {
00419 for (int j=0;j<NV;j++) {
00420 ierr = MatSetValue(
00421 Agradu,
00422 ND*t_v[j]+k,
00423 t_w[l],
00424 matele,
00425 ADD_VALUES
00426 );CHKERRQ(ierr);
00427 assert(t_v[j]<gdata->n_vert);
00428 assert(t_w[l]<n_vert);
00429
00430 #ifdef U1DIRICHLET
00431
00432 if (t_w[l] != U1DIRICHLET) {
00433 #endif
00434
00435 if (gdata->propdat[NP*gdata->eleprop[i]+4]>0.0) {
00436 ierr = MatSetValue(
00437 AdivM,
00438 t_w[l],
00439 ND*t_v[j]+k,
00440 matele*gdata->propdat[NP*gdata->eleprop[i]+4],
00441 ADD_VALUES
00442 );CHKERRQ(ierr);
00443 assert(t_w[l]<n_vert);
00444 assert(t_v[j]<gdata->n_vert);
00445 }
00446 #ifdef U1DIRICHLET
00447 }
00448 #endif
00449 }
00450 }
00451 }
00452 }
00453 }
00454 ierr = ProgressBar(1,1,10);
00455
00456 #undef MATFLUSH
00457 #ifdef MATFLUSH
00458 ierr = MatAssemblyBegin(Au1,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00459 ierr = MatAssemblyEnd(Au1,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00460 ierr = MatAssemblyBegin(Au2,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00461 ierr = MatAssemblyEnd(Au2,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00462 ierr = MatAssemblyBegin(ABu2,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00463 ierr = MatAssemblyEnd(ABu2,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00464 ierr = MatAssemblyBegin(Agradu,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00465 ierr = MatAssemblyEnd(Agradu,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00466 ierr = MatAssemblyBegin(AdivM,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00467 ierr = MatAssemblyEnd(AdivM,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00468 #endif
00469
00470 ierr = PetscGetTime(&t_t2);CHKERRQ(ierr);
00471 PetscPrintf(PETSC_COMM_WORLD,
00472 "stiffness matrix calculation took %g s = %g min\n",
00473 t_t2-t_t1,
00474 (t_t2-t_t1)/60.0
00475 );
00476
00477 ierr = VecRestoreArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);
00478
00479
00480 if (!rank) {
00481 #ifdef U1DIRICHLET
00482 ierr = MatSetValue(Au1,U1DIRICHLET,U1DIRICHLET,-1.0,ADD_VALUES);CHKERRQ(ierr);
00483 #endif
00484
00485 for (int i=0; i<gdata->n_vert; i++) {
00486 int j,k;
00487 j=i;
00488 k=j;
00489 if (gdata->vertbndg2bnd[j] >= 0) {
00490 assert(gdata->vertbndg2bnd[j]<gdata->n_vert_bnd);
00491
00492 ierr = MatSetValue(Au2,k,k,-1.0,ADD_VALUES);CHKERRQ(ierr);
00493 ierr = MatSetValue(ABu2,k,gdata->vertbndg2bnd[j],-1.0,ADD_VALUES);CHKERRQ(ierr);
00494 ierr = MatSetValue(ABu1,gdata->vertbndg2bnd[j],k,1.0,INSERT_VALUES);CHKERRQ(ierr);
00495 }
00496 }
00497 }
00498
00499 #undef DOMATCOMPRESS
00500 PetscPrintf(PETSC_COMM_WORLD,"Au1 matrix assembly complete\n");
00501 ierr = MatAssemblyBegin(Au1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00502 ierr = MatAssemblyEnd(Au1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00503
00504 #ifdef DOMATCOMPRESS
00505 ierr = MatCompress(Au1);CHKERRQ(ierr);
00506 #endif
00507 ierr = PrintMatInfoAll(Au1);CHKERRQ(ierr);
00508
00509 PetscPrintf(PETSC_COMM_WORLD,"Au2 matrix assembly complete\n");
00510 ierr = MatAssemblyBegin(Au2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00511 ierr = MatAssemblyEnd(Au2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00512 #ifdef DOMATCOMPRESS
00513 ierr = MatCompress(Au2);CHKERRQ(ierr);
00514 #endif
00515 ierr = PrintMatInfoAll(Au2);CHKERRQ(ierr);
00516
00517 PetscPrintf(PETSC_COMM_WORLD,"AdivM matrix assembly complete\n");
00518 ierr = MatAssemblyBegin(AdivM,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00519 ierr = MatAssemblyEnd(AdivM,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00520 #ifdef DOMATCOMPRESS
00521 ierr = MatCompress(AdivM);CHKERRQ(ierr);
00522 #endif
00523 ierr = PrintMatInfoAll(AdivM);CHKERRQ(ierr);
00524
00525 PetscPrintf(PETSC_COMM_WORLD,"Agradu matrix assembly complete\n");
00526 ierr = MatAssemblyBegin(Agradu,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00527 ierr = MatAssemblyEnd(Agradu,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00528 #ifdef DOMATCOMPRESS
00529 ierr = MatCompress(Agradu);CHKERRQ(ierr);
00530 #endif
00531 ierr = PrintMatInfoAll(Agradu);CHKERRQ(ierr);
00532
00533 PetscPrintf(PETSC_COMM_WORLD,"ABu1 matrix assembly complete\n");
00534 ierr = MatAssemblyBegin(ABu1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00535 ierr = MatAssemblyEnd(ABu1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00536 #ifdef DOMATCOMPRESS
00537 ierr = MatCompress(ABu1);CHKERRQ(ierr);
00538 #endif
00539 ierr = PrintMatInfoAll(ABu1);CHKERRQ(ierr);
00540
00541 PetscPrintf(PETSC_COMM_WORLD,"ABu2 matrix assembly complete\n");
00542 ierr = MatAssemblyBegin(ABu2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00543 ierr = MatAssemblyEnd(ABu2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00544 #ifdef DOMATCOMPRESS
00545 ierr = MatCompress(ABu2);CHKERRQ(ierr);
00546 #endif
00547 ierr = PrintMatInfoAll(ABu2);CHKERRQ(ierr);
00548
00549 Vec vertvolrec3;
00550 ierr = VecDuplicate(gdata->M,&vertvolrec3);CHKERRQ(ierr);
00551 for (int i=0;i<ND;i++) {
00552 ierr = VecStrideScatter(gdata->vertvol,i,vertvolrec3,INSERT_VALUES);CHKERRQ(ierr);
00553 }
00554 ierr = VecReciprocal(vertvolrec3);CHKERRQ(ierr);
00555 ierr = MatDiagonalScale(Agradu,vertvolrec3,PETSC_NULL);CHKERRQ(ierr);
00556
00557 ierr = VecDestroy(vertvolrec3);CHKERRQ(ierr);
00558
00559 ierr = MatScale(Au1,-1.0);CHKERRQ(ierr);
00560 ierr = MatScale(AdivM,-1.0);CHKERRQ(ierr);
00561
00562 ierr = MatScale(Au2,-1.0);CHKERRQ(ierr);
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580 #ifdef SUPERLU
00581 Mat Au1tmp;
00582 ierr = MatConvert(Au1,MATSUPERLU_DIST,MAT_INITIAL_MATRIX,&Au1tmp);CHKERRQ(ierr);
00583 ierr = MatDestroy(Au1);CHKERRQ(ierr);
00584 Au1=Au1tmp;
00585 #endif
00586
00587 #ifdef SUPERLU
00588 Mat Au2tmp;
00589 ierr = MatConvert(Au2,MATSUPERLU_DIST,MAT_INITIAL_MATRIX,&Au2tmp);CHKERRQ(ierr);
00590 ierr = MatDestroy(Au2);CHKERRQ(ierr);
00591 Au2=Au2tmp;
00592 #endif
00593
00594
00595 #if PETSC_VERSION >= 300
00596 ierr = MatSetOption(Au1,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
00597 ierr = MatSetOption(Au2,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
00598 ierr = MatSetOption(Au1,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr);
00599 ierr = MatSetOption(Au2,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr);
00600 #else
00601 ierr = MatSetOption(Au1,MAT_SYMMETRIC);CHKERRQ(ierr);
00602 ierr = MatSetOption(Au2,MAT_SYMMETRIC);CHKERRQ(ierr);
00603 ierr = MatSetOption(Au1,MAT_SYMMETRY_ETERNAL);CHKERRQ(ierr);
00604 ierr = MatSetOption(Au2,MAT_SYMMETRY_ETERNAL);CHKERRQ(ierr);
00605 #endif
00606
00607 ierr = VecDuplicate(utot,&u1);CHKERRQ(ierr);
00608 ierr = VecDuplicate(utot,&u2);CHKERRQ(ierr);
00609
00610 ierr = VecDuplicate(utot,&bu1);CHKERRQ(ierr);
00611 ierr = VecDuplicate(utot,&bu2);CHKERRQ(ierr);
00612
00613 ierr = VecCreate(PETSC_COMM_WORLD,&u1bnd);CHKERRQ(ierr);
00614 ierr = VecSetSizes(u1bnd,ln_bmatrixrows,gdata->n_vert_bnd);CHKERRQ(ierr);
00615 ierr = VecSetFromOptions(u1bnd);CHKERRQ(ierr);
00616
00617 ierr = VecDuplicate(u1bnd,&u2bnd);CHKERRQ(ierr);
00618
00619
00620
00621
00622 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_Au1);CHKERRQ(ierr);
00623 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_Au2);CHKERRQ(ierr);
00624 ierr = KSPSetOperators(ksp_Au1,Au1,Au1,SAME_PRECONDITIONER);CHKERRQ(ierr);
00625 ierr = KSPSetOperators(ksp_Au2,Au2,Au2,SAME_PRECONDITIONER);CHKERRQ(ierr);
00626
00627
00628
00629
00630 ierr = PetscOptionsHasName(PETSC_NULL,"-ksp_type",&flg);CHKERRQ(ierr);
00631 if (flg==PETSC_TRUE) {
00632 PetscPrintf(PETSC_COMM_WORLD,"Warning: Options '-ksp_*' are deprecated! Please use -hdemag_u{1,2}_ksp_* instead!\n");
00633 ierr = KSPSetFromOptions(ksp_Au1);CHKERRQ(ierr);
00634 ierr = KSPSetFromOptions(ksp_Au2);CHKERRQ(ierr);
00635 }
00636 else {
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647 for (int i=1;i<=2;i++) {
00648 char str[256],ostr[256];
00649
00650
00651 sprintf(ostr,"-hdemag_u%i_ksp_type",i);
00652 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00653 if (flg!=PETSC_TRUE) {
00654 sprintf(str,"cg");
00655 PetscOptionsSetValue(ostr,str);
00656 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00657 }
00658
00659 if (size>1) {
00660
00661 sprintf(ostr,"-hdemag_u%i_pc_type",i);
00662 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00663 if (flg!=PETSC_TRUE) {
00664 sprintf(str,"bjacobi");
00665 PetscOptionsSetValue(ostr,str);
00666 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00667 }
00668
00669
00670 sprintf(ostr,"-hdemag_u%i_sub_pc_type",i);
00671 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00672 if (flg!=PETSC_TRUE) {
00673 sprintf(str,"icc");
00674 PetscOptionsSetValue(ostr,str);
00675 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00676 }
00677
00678
00679 sprintf(ostr,"-hdemag_u%i_sub_pc_factor_shift_positive_definite",i);
00680 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00681 if (flg!=PETSC_TRUE) {
00682 sprintf(str,"1");
00683 PetscOptionsSetValue(ostr,str);
00684 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00685 }
00686 }
00687 else {
00688
00689 sprintf(ostr,"-hdemag_u%i_pc_type",i);
00690 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00691 if (flg!=PETSC_TRUE) {
00692 sprintf(str,"icc");
00693 PetscOptionsSetValue(ostr,str);
00694 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00695 }
00696
00697
00698 sprintf(ostr,"-hdemag_u%i_pc_factor_shift_positive_definite",i);
00699 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00700 if (flg!=PETSC_TRUE) {
00701 sprintf(str,"1");
00702 PetscOptionsSetValue(ostr,str);
00703 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00704 }
00705 }
00706 }
00707 }
00708
00709
00710 ierr = KSPSetOptionsPrefix(ksp_Au1,"hdemag_u1_");CHKERRQ(ierr);
00711
00712 ierr = KSPSetFromOptions(ksp_Au1);CHKERRQ(ierr);
00713
00714
00715 ierr = KSPSetOptionsPrefix(ksp_Au2,"hdemag_u2_");CHKERRQ(ierr);
00716
00717 ierr = KSPSetFromOptions(ksp_Au2);CHKERRQ(ierr);
00718
00719
00720 ierr = PetscOptionsHasName(PETSC_NULL,"-hdemag_u2_ksp_rtol",&flg);CHKERRQ(ierr);
00721
00722 if (flg!=PETSC_TRUE) {
00723 PetscReal rtol,atol,dtol;
00724 int maxits;
00725 ierr = KSPGetTolerances(ksp_Au1,&rtol,&atol,&dtol,&maxits);CHKERRQ(ierr);
00726
00727 PetscPrintf(PETSC_COMM_WORLD,"Adjusting tolerances for ksp_Au2\n");
00728 ierr = KSPSetTolerances(ksp_Au2,rtol*1e-3,atol,dtol,maxits);CHKERRQ(ierr);
00729 }
00730
00731
00732 ierr = KSPSetInitialGuessNonzero(ksp_Au1,PETSC_TRUE);CHKERRQ(ierr);
00733 ierr = KSPSetInitialGuessNonzero(ksp_Au2,PETSC_TRUE);CHKERRQ(ierr);
00734
00735
00736 ierr = KSPSetUp(ksp_Au1);CHKERRQ(ierr);
00737 ierr = KSPSetUp(ksp_Au2);CHKERRQ(ierr);
00738
00739 #if defined(SUPERLU)
00740 ierr = KSPGetTolerances(ksp_Au1,&rtol,&atol,&dtol,&maxits);CHKERRQ(ierr);
00741 PetscPrintf(PETSC_COMM_WORLD,"KSP for A*u1=divM:\n");
00742 PetscPrintf(PETSC_COMM_WORLD,"rtol: %g atol: %g dtol: %g maxits: %i\n",
00743 rtol,atol,dtol,maxits
00744 );
00745 ierr = KSPGetTolerances(ksp_Au2,&rtol,&atol,&dtol,&maxits);CHKERRQ(ierr);
00746 PetscPrintf(PETSC_COMM_WORLD,"KSP for Au2*u2=u2bnd:\n");
00747 PetscPrintf(PETSC_COMM_WORLD,"rtol: %g atol: %g dtol: %g maxits: %i\n",
00748 rtol,atol,dtol,maxits
00749 );
00750 #else
00751 PetscPrintf(PETSC_COMM_WORLD,"KSP for A*u1=divM:\n");
00752 ierr = KSPView(ksp_Au1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
00753 PetscPrintf(PETSC_COMM_WORLD,"KSP for Au2*u2=u2bnd:\n");
00754 ierr = KSPView(ksp_Au2,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
00755 #endif
00756
00757 ierr = VecDuplicate(gdata->M,&gdata->VHdem);CHKERRQ(ierr);
00758
00759 #ifdef HDEMAGUEPOL
00760 ierr = VecDuplicate(utot,&u1bak);CHKERRQ(ierr);
00761 ierr = VecDuplicate(utot,&u2bak);CHKERRQ(ierr);
00762 ierr = VecDuplicate(utot,&du1dt);CHKERRQ(ierr);
00763 ierr = VecDuplicate(utot,&du2dt);CHKERRQ(ierr);
00764 ierr = VecZeroEntries(du1dt);CHKERRQ(ierr);
00765 ierr = VecZeroEntries(du2dt);CHKERRQ(ierr);
00766 #endif
00767
00768 ierr = BMatrix(gdata,&Bmat);CHKERRQ(ierr);
00769
00770 MagparFunctionLogReturn(0);
00771 }
00772
00773
00774 int Hdemag(GridData *gdata,Vec VHtotsum)
00775 {
00776 #ifdef HDEMAGUEPOL
00777 static PetscReal tubak=0.0;
00778 PetscReal dt;
00779 #endif
00780
00781 int its1,its2;
00782
00783 MagparFunctionInfoBegin;
00784
00785 int rank,size;
00786 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00787 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00788
00789 if (doinit) {
00790 ierr = Hdemag_Init(gdata);CHKERRQ(ierr);
00791 doinit=0;
00792
00793 ierr = PetscGetTime(&t_t1);CHKERRQ(ierr);
00794 }
00795 if (!fieldon) {
00796 MagparFunctionInfoReturn(0);
00797 }
00798
00799 PetscLogDouble t_t3;
00800 t_t3=t_t1;
00801
00802
00803
00804
00805
00806
00807
00808 ierr = MatMult(AdivM,gdata->M,bu1);CHKERRQ(ierr);
00809
00810 #ifdef HDEMAGUEPOL
00811 dt=gdata->time-tubak;
00812 if (gdata->mode==0 && PetscAbsReal(dt) > 1e-8){
00813
00814 PetscReal dt2;
00815 dt2=dt*0.5;
00816 ierr = VecWAXPY(u1,dt2,du1dt,u1bak);CHKERRQ(ierr);
00817 }
00818 else {
00819 dt=0.0;
00820 }
00821 #endif
00822
00823 t_t2=t_t3;
00824 ierr = PetscGetTime(&t_t3);CHKERRQ(ierr);
00825
00826 PetscTruth oprofile;
00827 PetscOptionsHasName(PETSC_NULL,"-profile",&oprofile);
00828 if (oprofile) {
00829 PetscPrintf(PETSC_COMM_WORLD," timing: %s - div(M) - took %g s\n",__FUNCT__,t_t3-t_t2);
00830 } else {
00831 PetscInfo2(0," timing: %s - div(M) - took %g s\n",__FUNCT__,t_t3-t_t2);
00832 }
00833
00834
00835
00836
00837
00838 #ifndef U1DIRICHLET
00839
00840
00841
00842
00843 PetscReal u1a,*ta_u1;
00844 ierr = VecGetArray(u1,&ta_u1);CHKERRQ(ierr);
00845 u1a = -ta_u1[0];
00846 ierr = VecRestoreArray(u1,&ta_u1);CHKERRQ(ierr);
00847 if (fabs(u1a)>1e6) {
00848 ierr = VecShift(u1,u1a);CHKERRQ(ierr);
00849 PetscPrintf(PETSC_COMM_WORLD,"u1 shifted by %g\n",u1a);
00850 }
00851 #endif
00852
00853
00854
00855
00856 ierr = KSPSetInitialGuessNonzero(ksp_Au1,PETSC_TRUE);CHKERRQ(ierr);
00857
00858
00859
00860
00861
00862 ierr = KSPSolve(ksp_Au1,bu1,u1);CHKERRQ(ierr);
00863 ierr = KSPGetIterationNumber(ksp_Au1,(PetscInt*)&its1);CHKERRQ(ierr);
00864
00865 KSPConvergedReason reason;
00866 ierr = KSPGetConvergedReason(ksp_Au1,&reason);CHKERRQ(ierr);
00867
00868 PetscInfo2(0,"Au1*u1=divM: %i its, reason: %i\n",its1, reason);
00869 if (reason<0) {
00870 PetscReal rnorm;
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889 ierr = KSPGetResidualNorm(ksp_Au1,&rnorm);CHKERRQ(ierr);
00890 SETERRQ3(PETSC_ERR_CONV_FAILED,
00891 "Warning: Au1*u1=divM: KSP did not converge at t=%g ns - reason: %i, rnorm: %g\n",
00892 gdata->time*gdata->tscale*1e9,reason,rnorm
00893 );
00894 }
00895
00896 t_t2=t_t3;
00897 ierr = PetscGetTime(&t_t3);CHKERRQ(ierr);
00898
00899 if (oprofile) {
00900 PetscPrintf(PETSC_COMM_WORLD," timing: %s - u1 - took %g s\n",__FUNCT__,t_t3-t_t2);
00901 } else {
00902 PetscInfo2(0," timing: %s - u1 - took %g s\n",__FUNCT__,t_t3-t_t2);
00903 }
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913 if (its1>0) {
00914
00915
00916
00917
00918
00919 ierr = MatMult(ABu1,u1,u1bnd);CHKERRQ(ierr);
00920
00921 t_t2=t_t3;
00922 ierr = PetscGetTime(&t_t3);CHKERRQ(ierr);
00923
00924 if (oprofile) {
00925 PetscPrintf(PETSC_COMM_WORLD," timing: %s - u1 scatter - took %g s\n",__FUNCT__,t_t3-t_t2);
00926 } else {
00927 PetscInfo2(0," timing: %s - u1 scatter - took %g s\n",__FUNCT__,t_t3-t_t2);
00928 }
00929
00930
00931 ierr = MatMult(Bmat,u1bnd,u2bnd);CHKERRQ(ierr);
00932
00933 t_t2=t_t3;
00934 ierr = PetscGetTime(&t_t3);CHKERRQ(ierr);
00935
00936 if (oprofile) {
00937 PetscPrintf(PETSC_COMM_WORLD," timing: %s - u2=B*u1 - took %g s\n",__FUNCT__,t_t3-t_t2);
00938 } else {
00939 PetscInfo2(0," timing: %s - u2=B*u1 - took %g s\n",__FUNCT__,t_t3-t_t2);
00940 }
00941
00942
00943 ierr = MatMult(ABu2,u2bnd,bu2);CHKERRQ(ierr);
00944
00945 t_t2=t_t3;
00946 ierr = PetscGetTime(&t_t3);CHKERRQ(ierr);
00947
00948 if (oprofile) {
00949 PetscPrintf(PETSC_COMM_WORLD," timing: %s - u2 scatter - took %g s\n",__FUNCT__,t_t3-t_t2);
00950 } else {
00951 PetscInfo2(0," timing: %s - u2 scatter - took %g s\n",__FUNCT__,t_t3-t_t2);
00952 }
00953
00954 #ifdef HDEMAGUEPOL
00955
00956
00957
00958
00959
00960
00961
00962
00963 #endif
00964
00965
00966
00967
00968
00969 ierr = KSPSetInitialGuessNonzero(ksp_Au2,PETSC_TRUE);CHKERRQ(ierr);
00970
00971
00972
00973
00974
00975 ierr = KSPSolve(ksp_Au2,bu2,u2);CHKERRQ(ierr);
00976 ierr = KSPGetIterationNumber(ksp_Au2,(PetscInt*)&its2);CHKERRQ(ierr);
00977 ierr = KSPGetConvergedReason(ksp_Au2,&reason);CHKERRQ(ierr);
00978
00979 PetscInfo2(0,"Au2*u2=bu2: %i its, reason: %i\n",its2, reason);
00980 if (reason<0) {
00981 PetscReal rnorm;
00982
00983 ierr = KSPGetResidualNorm(ksp_Au2,&rnorm);CHKERRQ(ierr);
00984 PetscPrintf(PETSC_COMM_WORLD,
00985 "Warning: Au2*u2=bu2: KSP did not converge at t=%g ns - reason: %i, rnorm: %g\n",
00986 gdata->time*gdata->tscale*1e9,reason,rnorm
00987 );
00988 }
00989
00990 t_t2=t_t3;
00991 ierr = PetscGetTime(&t_t3);CHKERRQ(ierr);
00992
00993 if (oprofile) {
00994 PetscPrintf(PETSC_COMM_WORLD," timing: %s - u2 - took %g s\n",__FUNCT__,t_t3-t_t2);
00995 } else {
00996 PetscInfo2(0," timing: %s - u2 - took %g s\n",__FUNCT__,t_t3-t_t2);
00997 }
00998
00999
01000
01001
01002
01003 ierr = VecWAXPY(utot,1.0,u1,u2);CHKERRQ(ierr);
01004 ierr = MatMult(Agradu,utot,gdata->VHdem);CHKERRQ(ierr);
01005
01006 #ifdef HDEMAGUEPOL
01007 if (gdata->mode==0 && PetscAbsReal(dt) > 1e-8){
01008 dt=1.0/dt;
01009
01010
01011 ierr = VecCopy(u1,du1dt);CHKERRQ(ierr);
01012 ierr = VecAXPY(du1dt,c_mone,u1bak);CHKERRQ(ierr);
01013 ierr = VecScale(du1dt,dt);CHKERRQ(ierr);
01014 ierr = VecCopy(u1,u1bak);CHKERRQ(ierr);
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024 tubak=gdata->time;
01025 }
01026 #endif
01027
01028
01029 }
01030 else {
01031 PetscInfo2(0,
01032 "Info (t=%g ns): Skipped calculation of u2: %i\n",
01033 gdata->time*gdata->tscale*1e9,
01034 its1
01035 );
01036 }
01037
01038 ierr = VecAXPY(VHtotsum,1.0,gdata->VHdem);CHKERRQ(ierr);
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059 MagparFunctionProfReturn(0);
01060 }
01061
01062
01063 int Hdemag_Energy(GridData *gdata,Vec VMom,PetscReal *energy)
01064 {
01065 PetscReal e;
01066
01067 MagparFunctionInfoBegin;
01068
01069 if (!fieldon) {
01070 *energy=0.0;
01071 MagparFunctionInfoReturn(0);
01072 }
01073
01074 ierr = VecDot(VMom,gdata->VHdem,&e);CHKERRQ(ierr);
01075 e /= -2.0;
01076
01077 *energy=e;
01078
01079 MagparFunctionProfReturn(0);
01080 }
01081