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