hexch_ani.c

Go to the documentation of this file.
00001 /*
00002     This file is part of magpar.
00003 
00004     Copyright (C) 2002-2010 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 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   /* 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 
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   /* get ownership range */
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       /* we could also use MPI_Send/Recv to get everything just to the first processor */
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   /* get ownership range */
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   /* calculate number of nonzeroes for each line */
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   /* create matrix to calculate Ad2E_dMidMk */
00150 /* normal MPIAIJ format: */
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 /* block MPIAIJ format:
00160   ierr = MatCreateMPIBAIJ(
00161     PETSC_COMM_WORLD,
00162     ND,
00163     ND*gdata->ln_vert,ND*gdata->ln_vert,
00164     ND*gdata->n_vert,ND*gdata->n_vert,
00165     PREALLOC_DG,PETSC_NULL,PREALLOC_OD,PETSC_NULL,
00166     &Ad2E_dMidMk
00167   );CHKERRQ(ierr);
00168   ierr = MatSetFromOptions(Ad2E_dMidMk);CHKERRQ(ierr);
00169 */
00170 
00171 #ifdef EXCH
00172   /* calculate number of nonzeroes for each line */
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   /* matrix for exchange field */
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     /* skip elements with Ms<=0 */
00218     if (gdata->propdat[NP*gdata->eleprop[i]+4]<=0.0) {
00219       continue;
00220     }
00221 
00222     /* get element volume */
00223     PetscReal t_elevol;
00224     t_elevol=ta_elevol[i];
00225 
00226     /* get vertex ids */
00227     int *t_v;
00228     t_v=gdata->elevert+i*NV;
00229 
00230     /* Get Vertex coordinates */
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     /* set matrix elements */
00241 
00242     /* prefactor for exchange field: pre_exch=2*A */
00243     PetscReal   pre_exch;
00244     pre_exch=2.0*gdata->propdat[NP*gdata->eleprop[i]+5];
00245 /* must not divide by Js to account for different materials:
00246    |M|=1 -> Js has to be considered in local fields !!!
00247       /gdata->propdat[NP*gdata->eleprop[i]+4];
00248 */
00249     for (int l=0;l<NV;l++) {
00250       for (int j=0;j<NV;j++) {
00251         PetscReal matele;
00252 
00253         /* set matrix element for jacobian */
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     /* skip elements with Ms<=0 */
00300     if (gdata->propdat[NP*gdata->eleprop[i]+4]<=0.0) {
00301       continue;
00302     }
00303     /* skip elements with K1==0 */
00304     if (PetscAbsReal(gdata->propdat[NP*gdata->eleprop[i]+2])<D_EPS) {
00305       continue;
00306     }
00307     /* skip matrix element if material has cubic anisotropy */
00308     if (gdata->propdat[NP*gdata->eleprop[i]+16]) {
00309       continue;
00310     }
00311 
00312     /* prefactor for anisotropy field: pre_ani=2*K1*V */
00313     pre_ani=2.0*ta_elevol[i]*gdata->propdat[NP*gdata->eleprop[i]+2];
00314     /* division by Js to account for different materials done below (vertMomrec3) */
00315 
00316     Eanioffset += pre_ani/2.0;
00317 
00318     /* TODO: more elegant with MatSetValuesLocal !? */
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           /* calculate Int(phi[i]*phi[j])
00332 
00333              i==j ? Int(phi[i]*phi[j])=V/10 : Int(phi[i]*phi[j])=V/20
00334              source: Advanced Finite Element Methods (ASEN 6367) - Spring 2009
00335                      Aerospace Engineering Sciences - University of Colorado at Boulder
00336                      Chapter 16: THE LINEAR TETRAHEDRON
00337                      page 10, Eq. (16.35)
00338                      http://www.colorado.edu/engineering/cas/courses.d/AFEM.d/AFEM.Ch16.d/AFEM.Ch16.pdf
00339 
00340                      http://smirnov.mae.wvu.edu/hedra/poisson.pdf
00341 
00342                      Wenjie Chen, D. R. Fredkin, and T. R. Koehler,
00343                      "A New Finite Element Method in Micromagnetics,"
00344                      IEEE Trans. Magn. 29 (1993) 2124.
00345                      http://ieeexplore.ieee.org/iel1/20/5774/00221033.pdf
00346 
00347                      see also hdemag.c
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   /* calculate field */
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   /* calculate exchange field separately */
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   /* calculate exchange+anisotropy energy */
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   /* calculate exchange energy */
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 

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