hexch_ani.c

Go to the documentation of this file.
00001 /*
00002     This file is part of magpar.
00003 
00004     Copyright (C) 2002-2009 Werner Scholz
00005 
00006     www:   http://www.magpar.net/
00007     email: magpar(at)magpar.net
00008 
00009     magpar is free software; you can redistribute it and/or modify
00010     it under the terms of the GNU General Public License as published by
00011     the Free Software Foundation; either version 2 of the License, or
00012     (at your option) any later version.
00013 
00014     magpar is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00017     GNU General Public License for more details.
00018 
00019     You should have received a copy of the GNU General Public License
00020     along with magpar; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
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   /* create gdata->VHexchani in any case
00056      otherwise we get SEGV in WriteAVS
00057      TODO: not most elegant way to solve this issue
00058   */
00059   ierr = VecDuplicate(gdata->M,&gdata->VHexchani);CHKERRQ(ierr);
00060 
00061   /* activate only if Ms>0 && (K1!=0 or A!=0) */
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   /* get ownership range */
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       /* we could also use MPI_Send/Recv to get everything just to the first processor */
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   /* get ownership range */
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   /* calculate number of nonzeroes for each line */
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   /* create matrix to calculate Ad2E_dMidMk */
00144 /* normal MPIAIJ format: */
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 /* block MPIAIJ format:
00154   ierr = MatCreateMPIBAIJ(
00155     PETSC_COMM_WORLD,
00156     ND,
00157     ND*gdata->ln_vert,ND*gdata->ln_vert,
00158     ND*gdata->n_vert,ND*gdata->n_vert,
00159     PREALLOC_DG,PETSC_NULL,PREALLOC_OD,PETSC_NULL,
00160     &Ad2E_dMidMk
00161   );CHKERRQ(ierr);
00162   ierr = MatSetFromOptions(Ad2E_dMidMk);CHKERRQ(ierr);
00163 */
00164   ierr = VecDuplicate(gdata->M,&gdata->VHexchani);CHKERRQ(ierr);
00165 
00166 #ifdef EXCH
00167   /* calculate number of nonzeroes for each line */
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   /* matrix for exchange field */
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     /* skip elements with Ms<=0 */
00213     if (gdata->propdat[NP*gdata->eleprop[i]+4]<=0.0) {
00214       continue;
00215     }
00216 
00217     /* get element volume */
00218     PetscReal t_elevol;
00219     t_elevol=ta_elevol[i];
00220 
00221     /* get vertex ids */
00222     int *t_v;
00223     t_v=gdata->elevert+i*NV;
00224 
00225     /* Get Vertex coordinates */
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     /* set matrix elements */
00236 
00237     /* prefactor for exchange field: pre_exch=2*A */
00238     PetscReal   pre_exch;
00239     pre_exch=2.0*gdata->propdat[NP*gdata->eleprop[i]+5];
00240 /* must not divide by Js to account for different materials:
00241    |M|=1 -> Js has to be considered in local fields !!!
00242       /gdata->propdat[NP*gdata->eleprop[i]+4];
00243 */
00244     for (int l=0;l<NV;l++) {
00245       for (int j=0;j<NV;j++) {
00246         PetscReal matele;
00247 
00248         /* set matrix element for jacobian */
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     /* skip elements with Ms<=0 */
00295     if (gdata->propdat[NP*gdata->eleprop[i]+4]<=0.0) {
00296       continue;
00297     }
00298     /* skip elements with K1==0 */
00299     if (PetscAbsReal(gdata->propdat[NP*gdata->eleprop[i]+2])<D_EPS) {
00300       continue;
00301     }
00302     /* skip matrix element if material has cubic anisotropy */
00303     if (gdata->propdat[NP*gdata->eleprop[i]+16]) {
00304       continue;
00305     }
00306 
00307     /* prefactor for anisotropy field: pre_ani=2*K1*V */
00308     pre_ani=2.0*ta_elevol[i]*gdata->propdat[NP*gdata->eleprop[i]+2];
00309     /* division by Js to account for different materials done below (vertMomrec3) */
00310 
00311     Eanioffset += pre_ani/2.0;
00312 
00313     /* TODO: more elegant with MatSetValuesLocal !? */
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           /* approximation for Int(phi[i]*phi[j])
00327 
00328              i==j ? Int(phi[i]*phi[j])=V/10 : Int(phi[i]*phi[j])=V/20
00329              source: http://caswww.colorado.edu/courses.d/AFEM.d/AFEM.Ch15.d/
00330                      http://smirnov.mae.wvu.edu/hedra/poisson.pdf
00331                      Wenjie Chen, D. R. Fredkin, and T. R. Koehler,
00332                      "A New Finite Element Method in Micromagnetics,"
00333                      IEEE Trans. Magn. 29 (1993) 2124.
00334                      http://ieeexplore.ieee.org/iel1/20/5774/00221033.pdf
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   /* calculate field */
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   /* calculate exchange field separately */
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   /* calculate exchange+anisotropy energy */
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   /* calculate exchange energy */
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 

magpar - Parallel Finite Element Micromagnetics Package
Copyright (C) 2002-2009 Werner Scholz