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: hexch_ani.c 3023 2010-03-27 23:42:04Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "field.h"
00028 #include "init/init.h"
00029 #include "util/util.h"
00030
00031 static int doinit=1;
00032 static int fieldon=0;
00033
00034 Mat Ad2E_dMidMk;
00035 static PetscReal Eanioffset;
00037 #ifdef EXCH
00038 static Mat AHexchonly;
00039 #endif
00040
00041 #define PREALLOC_DG 15
00042 #define PREALLOC_OD 5
00043
00044 int Hexchani_Init(GridData *gdata)
00045 {
00046 MagparFunctionLogBegin;
00047
00048 int rank,size;
00049 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00050 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00051
00052 fieldon=0;
00053
00054 Ad2E_dMidMk=PETSC_NULL;
00055
00056
00057
00058
00059 ierr = VecDuplicate(gdata->M,&gdata->VHexchani);CHKERRQ(ierr);
00060
00061
00062 for (int i=0;i<gdata->n_prop;i++) {
00063 if (gdata->propdat[NP*i+4]>0.0 && (gdata->propdat[NP*i+2]!=0.0 || gdata->propdat[NP*i+5]!=0.0)) {
00064 fieldon=1;
00065 break;
00066 }
00067 }
00068 if (!fieldon) {
00069 PetscPrintf(PETSC_COMM_WORLD,"Hexchani off\n");
00070 ierr = PetscOptionsSetValue("-precon","0");CHKERRQ(ierr);
00071 PetscPrintf(PETSC_COMM_WORLD,"Disabling preconditioning\n");
00072
00073 gdata->VHexchani=PETSC_NULL;
00074 #ifdef EXCH
00075 gdata->VHexchonly=PETSC_NULL;
00076 #endif
00077
00078 MagparFunctionLogReturn(0);
00079 }
00080 PetscPrintf(PETSC_COMM_WORLD,"Hexchani on\n");
00081
00082 PetscPrintf(PETSC_COMM_WORLD,"Calculating sparsity pattern for matrix preallocation\n");
00083
00084
00085 int low,high;
00086 ierr = VecGetOwnershipRange(gdata->elevol,&low,&high);CHKERRQ(ierr);
00087
00088 int *elevertall;
00089
00090 #ifdef HAVE_ELEVERTALL
00091 elevertall=gdata->elevertall;
00092 #else
00093 if (size>1) {
00094 ierr = PetscMalloc(NV*gdata->n_ele*sizeof(int),&elevertall);CHKERRQ(ierr);
00095 ierr = PetscMemcpy(elevertall+NV*low,gdata->elevert,NV*gdata->ln_ele*sizeof(int));CHKERRQ(ierr);
00096
00097 int *p;
00098 p=elevertall;
00099
00100 int recvcount[size];
00101 recvcount[rank]=NV*gdata->ln_ele;
00102 for (int i=0;i<size;i++) {
00103 ierr = MPI_Bcast(recvcount+i,1,MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
00104
00105 ierr = MPI_Bcast(p,recvcount[i],MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
00106 p+=recvcount[i];
00107 }
00108 assert(p==elevertall+NV*gdata->n_ele);
00109 }
00110 else {
00111 elevertall=gdata->elevert;
00112 }
00113 #endif
00114
00115 int *ia,*ja;
00116 #ifdef HAVE_M2IJ
00117 ia=gdata->mesh2nodal_ia;
00118 ja=gdata->mesh2nodal_ja;
00119 #else
00120 ierr = Mesh2Nodal(gdata->n_ele,gdata->n_vert,elevertall,&ia,&ja);CHKERRQ(ierr);
00121 #endif
00122
00123
00124 ierr = VecGetOwnershipRange(gdata->M,&low,&high);CHKERRQ(ierr);
00125 low/=ND;
00126 high/=ND;
00127
00128 int *d_nnz;
00129 int *o_nnz;
00130 ierr = PetscMalloc(ND*(high-low)*sizeof(int),&d_nnz);CHKERRQ(ierr);
00131 ierr = PetscMalloc(ND*(high-low)*sizeof(int),&o_nnz);CHKERRQ(ierr);
00132
00133 for (int i=low;i<high;i++) {
00134 for (int j=0;j<ND;j++) {
00135 d_nnz[ND*(i-low)+j]=ND;
00136 o_nnz[ND*(i-low)+j]=0;
00137 for (int l=ia[i];l<ia[i+1];l++) {
00138 if (ja[l]>=low && ja[l]<high) {
00139 d_nnz[ND*(i-low)+j]+=ND;
00140 }
00141 else {
00142 o_nnz[ND*(i-low)+j]+=ND;
00143 }
00144 }
00145 }
00146 }
00147
00148 PetscPrintf(PETSC_COMM_WORLD,"Creating matrices\n");
00149
00150
00151 ierr = MatCreateMPIAIJ(
00152 PETSC_COMM_WORLD,
00153 ND*gdata->ln_vert,ND*gdata->ln_vert,
00154 ND*gdata->n_vert,ND*gdata->n_vert,
00155 0,d_nnz,0,o_nnz,
00156 &Ad2E_dMidMk
00157 );CHKERRQ(ierr);
00158 ierr = MatSetFromOptions(Ad2E_dMidMk);CHKERRQ(ierr);
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171 #ifdef EXCH
00172
00173 for (int i=low;i<high;i++) {
00174 int k;
00175 k=ia[i+1]-ia[i]+1;
00176 for (int j=0;j<ND;j++) {
00177 if (i>=low || i<high) {
00178 d_nnz[ND*(i-low)+j]=k;
00179 }
00180 else {
00181 o_nnz[ND*(i-low)+j]=k;
00182 }
00183 }
00184 }
00185
00186
00187 ierr = MatCreateMPIAIJ(
00188 PETSC_COMM_WORLD,
00189 ND*gdata->ln_vert,ND*gdata->ln_vert,
00190 ND*gdata->n_vert,ND*gdata->n_vert,
00191 0,d_nnz,0,o_nnz,
00192 &AHexchonly
00193 );CHKERRQ(ierr);
00194 ierr = MatSetFromOptions(AHexchonly);CHKERRQ(ierr);
00195
00196 ierr = VecDuplicate(gdata->M,&gdata->VHexchonly);CHKERRQ(ierr);
00197 #endif
00198
00199 #ifndef HAVE_M2IJ
00200 ierr = PetscFree(ia);CHKERRQ(ierr);
00201 ierr = PetscFree(ja);CHKERRQ(ierr);
00202 #endif
00203 ierr = PetscFree(d_nnz);CHKERRQ(ierr);
00204 ierr = PetscFree(o_nnz);CHKERRQ(ierr);
00205 #ifndef HAVE_ELEVERTALL
00206 if (size>1) {
00207 ierr = PetscFree(elevertall);CHKERRQ(ierr);
00208 }
00209 #endif
00210
00211 PetscPrintf(PETSC_COMM_WORLD,"Calculating matrix elements\n");
00212
00213 PetscReal *ta_elevol;
00214 ierr = VecGetArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);
00215
00216 for (int i=0; i<gdata->ln_ele; i++) {
00217
00218 if (gdata->propdat[NP*gdata->eleprop[i]+4]<=0.0) {
00219 continue;
00220 }
00221
00222
00223 PetscReal t_elevol;
00224 t_elevol=ta_elevol[i];
00225
00226
00227 int *t_v;
00228 t_v=gdata->elevert+i*NV;
00229
00230
00231 PetscReal *x1,*x2,*x3,*x4;
00232 x1=gdata->vertxyz+(t_v[0])*ND;
00233 x2=gdata->vertxyz+(t_v[1])*ND;
00234 x3=gdata->vertxyz+(t_v[2])*ND;
00235 x4=gdata->vertxyz+(t_v[3])*ND;
00236
00237 PetscReal D_etaj[NV][ND];
00238 tetgrad(x1,x2,x3,x4,D_etaj);
00239
00240
00241
00242
00243 PetscReal pre_exch;
00244 pre_exch=2.0*gdata->propdat[NP*gdata->eleprop[i]+5];
00245
00246
00247
00248
00249 for (int l=0;l<NV;l++) {
00250 for (int j=0;j<NV;j++) {
00251 PetscReal matele;
00252
00253
00254 matele=-t_elevol*my_ddot(ND,D_etaj[j],1,D_etaj[l],1)*pre_exch;
00255 if (PetscAbsReal(matele)>D_EPS) {
00256 for (int k=0;k<ND;k++) {
00257 ierr = MatSetValue(
00258 Ad2E_dMidMk,
00259 ND*t_v[j]+k,
00260 ND*t_v[l]+k,
00261 matele,
00262 ADD_VALUES
00263 );CHKERRQ(ierr);
00264 #ifdef EXCH
00265 ierr = MatSetValue(
00266 AHexchonly,
00267 ND*t_v[j]+k,
00268 ND*t_v[l]+k,
00269 matele,
00270 ADD_VALUES
00271 );CHKERRQ(ierr);
00272 #endif
00273 }
00274 }
00275 }
00276 }
00277 }
00278
00279 ierr = MatAssemblyBegin(Ad2E_dMidMk,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00280 ierr = MatAssemblyEnd(Ad2E_dMidMk,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00281
00282 #ifdef EXCH
00283 ierr = MatAssemblyBegin(AHexchonly,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00284 ierr = MatAssemblyEnd(AHexchonly,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00285 #endif
00286
00287 ierr = PetscGetTime(&t_t2);CHKERRQ(ierr);
00288 PetscPrintf(PETSC_COMM_WORLD,
00289 "exchange matrix calculation took %g s = %g min\n",
00290 t_t2-t_t1,
00291 (t_t2-t_t1)/60.0
00292 );
00293
00294
00295 Eanioffset = 0.0;
00296
00297 PetscReal pre_ani;
00298 for (int i=0;i<gdata->ln_ele;i++) {
00299
00300 if (gdata->propdat[NP*gdata->eleprop[i]+4]<=0.0) {
00301 continue;
00302 }
00303
00304 if (PetscAbsReal(gdata->propdat[NP*gdata->eleprop[i]+2])<D_EPS) {
00305 continue;
00306 }
00307
00308 if (gdata->propdat[NP*gdata->eleprop[i]+16]) {
00309 continue;
00310 }
00311
00312
00313 pre_ani=2.0*ta_elevol[i]*gdata->propdat[NP*gdata->eleprop[i]+2];
00314
00315
00316 Eanioffset += pre_ani/2.0;
00317
00318
00319 for (int k=0;k<ND;k++) {
00320 for (int l=0;l<ND;l++) {
00321 PetscReal matele;
00322
00323 matele=
00324 pre_ani*
00325 gdata->propdat[NP*gdata->eleprop[i]+6+k]*
00326 gdata->propdat[NP*gdata->eleprop[i]+6+l];
00327
00328 if (PetscAbsReal(matele)<D_EPS) continue;
00329
00330 for (int m=0;m<NV;m++) {
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349 for (int j=0;j<NV;j++) {
00350 PetscReal matele2;
00351
00352 if (j==m) {
00353 matele2=matele/10.0;
00354 }
00355 else {
00356 matele2=matele/20.0;
00357 }
00358
00359 ierr = MatSetValue(
00360 Ad2E_dMidMk,
00361 gdata->elevert[NV*i+j]*ND+k,
00362 gdata->elevert[NV*i+m]*ND+l,
00363 matele2,
00364 ADD_VALUES
00365 );CHKERRQ(ierr);
00366 }
00367 }
00368 }
00369 }
00370 }
00371
00372 ierr = MatAssemblyBegin(Ad2E_dMidMk,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00373 ierr = MatAssemblyEnd(Ad2E_dMidMk,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00374
00375 PetscLogDouble t_t3;
00376 ierr = PetscGetTime(&t_t3);CHKERRQ(ierr);
00377 PetscPrintf(PETSC_COMM_WORLD,
00378 "anisotropy matrix calculation took %g s = %g min\n",
00379 t_t3-t_t2,
00380 (t_t3-t_t2)/60.0
00381 );
00382
00383 ierr = VecRestoreArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);
00384
00385 PetscReal tmp;
00386 ierr = PetscGlobalSum(&Eanioffset,&tmp,PETSC_COMM_WORLD);CHKERRQ(ierr);
00387 Eanioffset=tmp;
00388
00389 PetscPrintf(PETSC_COMM_WORLD,"Eanioffset: %g = %g J = %g J/m^3\n",
00390 Eanioffset,
00391 Eanioffset*gdata->escale,
00392 Eanioffset*gdata->escale/gdata->totvol
00393 );
00394
00395 Vec vertMomrec3;
00396 ierr = VecDuplicate(gdata->M,&vertMomrec3);CHKERRQ(ierr);
00397 ierr = VecCopy(gdata->VMs3,vertMomrec3);CHKERRQ(ierr);
00398 ierr = VecReciprocal(vertMomrec3);CHKERRQ(ierr);
00399
00400 PetscInfo(0,"Ad2E_dMidMk matrix assembly complete\n");
00401 ierr = MatAssemblyBegin(Ad2E_dMidMk,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00402 ierr = MatAssemblyEnd(Ad2E_dMidMk,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00403
00404 ierr = MatDiagonalScale(Ad2E_dMidMk,vertMomrec3,PETSC_NULL);CHKERRQ(ierr);
00405 ierr = PrintMatInfoAll(Ad2E_dMidMk);CHKERRQ(ierr);
00406
00407 #ifdef EXCH
00408 PetscInfo(0,"AHexchonly matrix assembly complete\n");
00409 ierr = MatAssemblyBegin(AHexchonly,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00410 ierr = MatAssemblyEnd(AHexchonly,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00411 ierr = MatDiagonalScale(AHexchonly,vertMomrec3,PETSC_NULL);CHKERRQ(ierr);
00412 ierr = PrintMatInfoAll(AHexchonly);CHKERRQ(ierr);
00413 #endif
00414
00415 ierr = VecDestroy(vertMomrec3);CHKERRQ(ierr);
00416
00417 MagparFunctionLogReturn(0);
00418 }
00419
00420
00421 int Hexchani(GridData *gdata,Vec VHtotsum)
00422 {
00423 MagparFunctionInfoBegin;
00424
00425 if (doinit) {
00426 ierr = Hexchani_Init(gdata);CHKERRQ(ierr);
00427 doinit=0;
00428 }
00429 if (!fieldon) {
00430 MagparFunctionInfoReturn(0);
00431 }
00432
00433
00434 ierr = MatMult(Ad2E_dMidMk,gdata->M,gdata->VHexchani);CHKERRQ(ierr);
00435
00436 ierr = VecAXPY(VHtotsum,1.0,gdata->VHexchani);CHKERRQ(ierr);
00437
00438 #ifdef EXCH
00439
00440 ierr = MatMult(AHexchonly,gdata->M,gdata->VHexchonly);CHKERRQ(ierr);
00441 #endif
00442
00443 MagparFunctionProfReturn(0);
00444 }
00445
00446
00447 int Hexchani_Energy(GridData *gdata,Vec VMom,PetscReal *energy)
00448 {
00449 PetscReal e;
00450
00451 MagparFunctionInfoBegin;
00452
00453 if (!fieldon) {
00454 *energy=0.0;
00455 MagparFunctionInfoReturn(0);
00456 }
00457
00458
00459 ierr = VecDot(VMom,gdata->VHexchani,&e);CHKERRQ(ierr);
00460 e /= -2.0;
00461 e += Eanioffset;
00462
00463 *energy=e;
00464
00465 MagparFunctionProfReturn(0);
00466 }
00467
00468
00469 #ifdef EXCH
00470 int Hexchonly_Energy(GridData *gdata,Vec VMom,PetscReal *energy)
00471 {
00472 PetscReal e;
00473
00474 MagparFunctionInfoBegin;
00475
00476 if (!fieldon) {
00477 *energy=0.0;
00478 MagparFunctionInfoReturn(0);
00479 }
00480
00481
00482 ierr = VecDot(VMom,gdata->VHexchonly,&e);CHKERRQ(ierr);
00483 e /= -2.0;
00484
00485 *energy=e;
00486
00487 MagparFunctionProfReturn(0);
00488 }
00489 #endif
00490