hcubic.c

Go to the documentation of this file.
00001 /*
00002     This file is part of magpar.
00003 
00004     Copyright (C) 2005 Greg Parker
00005     Copyright (C) 2005-2010 Werner Scholz
00006 
00007     www:   http://www.magpar.net/
00008     email: magpar(at)magpar.net
00009 
00010     magpar is free software; you can redistribute it and/or modify
00011     it under the terms of the GNU General Public License as published by
00012     the Free Software Foundation; either version 2 of the License, or
00013     (at your option) any later version.
00014 
00015     magpar is distributed in the hope that it will be useful,
00016     but WITHOUT ANY WARRANTY; without even the implied warranty of
00017     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018     GNU General Public License for more details.
00019 
00020     You should have received a copy of the GNU General Public License
00021     along with magpar; if not, write to the Free Software
00022     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00023 */
00024 
00025 static char const Id[] = "$Id: hcubic.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
00026 static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";
00027 
00028 #include "field.h"
00029 
00030 static int       doinit=1;
00031 static int       fieldon=0;
00032 static Vec       VHcubic;
00033 #ifdef EXCH
00034 static PetscReal Ecubic;
00035 #endif
00036 
00037 int Hcubic_Init(GridData *gdata)
00038 {
00039   MagparFunctionLogBegin;
00040 
00041   for (int i=0;i<gdata->n_prop;i++) {
00042     if (gdata->propdat[NP*i+16]==1 ||
00043         PetscAbsReal(gdata->propdat[NP*i+ 3]) > D_EPS) {
00044       /* FIXME: latter case not really cubic, but hcubic.c calculates the K_2 term
00045          anisotropy field calculation then off!? Hani calculated twice?
00046       */
00047       fieldon=1;
00048       break;
00049     }
00050   }
00051 
00052   if (fieldon) {
00053     PetscPrintf(PETSC_COMM_WORLD,"Hcubic on\n");
00054     ierr = VecDuplicate(gdata->M,&VHcubic);CHKERRQ(ierr);
00055   }
00056   else {
00057     PetscPrintf(PETSC_COMM_WORLD,"Hcubic off\n");
00058   }
00059 
00060   MagparFunctionLogReturn(0);
00061 }
00062 
00063 
00064 int Hcubic(GridData *gdata,Vec VHtotsum)
00065 {
00066   MagparFunctionInfoBegin;
00067 
00068   /* initialize if necessary */
00069   if (doinit) {
00070     ierr = Hcubic_Init(gdata);CHKERRQ(ierr);
00071     doinit=0;
00072   }
00073   if (!fieldon) {
00074     MagparFunctionInfoReturn(0);
00075   }
00076 
00077   PetscReal *ta_Hcubic,*ta_M;
00078   ierr = VecGetArray(VHcubic,&ta_Hcubic);CHKERRQ(ierr);
00079   ierr = VecGetArray(gdata->M,&ta_M);CHKERRQ(ierr);
00080 
00081 #ifdef EXCH
00082   PetscReal *ta_VMs3;
00083   ierr = VecGetArray(gdata->VMs3,&ta_VMs3);CHKERRQ(ierr);
00084   Ecubic=0.0;
00085 #endif
00086 
00087   /* TODO: for nodes at interfaces, the assignment is "at will" */
00088   for (int j=0; j<gdata->ln_vert; j++) {
00089     if (gdata->propdat[NP*gdata->vertprop[j]+16]) {
00090 
00091       /* only collect cubic Hk if 'grain' has cubic Hk */
00092       PetscReal k1,k2,jsi;
00093       jsi=gdata->propdat[NP*gdata->vertprop[j]+4];
00094       /* FIXME: why 1e-20? */
00095       jsi /= (jsi*jsi+1.e-20);
00096       k1=gdata->propdat[NP*gdata->vertprop[j]+2]*jsi;
00097       k2=gdata->propdat[NP*gdata->vertprop[j]+3]*jsi;
00098 
00099       /* a-axis of cubic crystal */
00100       PetscReal ax,ay,az;
00101       ax=gdata->propdat[NP*gdata->vertprop[j]+6];
00102       ay=gdata->propdat[NP*gdata->vertprop[j]+7];
00103       az=gdata->propdat[NP*gdata->vertprop[j]+8];
00104 
00105       /* b-axis of cubic crystal */
00106       PetscReal bx,by,bz;
00107       bx=gdata->propdat[NP*gdata->vertprop[j]+10];
00108       by=gdata->propdat[NP*gdata->vertprop[j]+11];
00109       bz=gdata->propdat[NP*gdata->vertprop[j]+12];
00110 
00111       /* c-axis of cubic crystal */
00112       PetscReal cx,cy,cz;
00113       cx=gdata->propdat[NP*gdata->vertprop[j]+13];
00114       cy=gdata->propdat[NP*gdata->vertprop[j]+14];
00115       cz=gdata->propdat[NP*gdata->vertprop[j]+15];
00116 
00117       /* get components of M */
00118       PetscReal mx,my,mz;
00119       mx=ta_M[ND*j+0];
00120       my=ta_M[ND*j+1];
00121       mz=ta_M[ND*j+2];
00122 
00123       /* direction cosines and their squares */
00124       PetscReal a1,a2,a3;
00125       a1=mx*ax+my*ay+mz*az;
00126       a2=mx*bx+my*by+mz*bz;
00127       a3=mx*cx+my*cy+mz*cz;
00128 
00129       PetscReal a12,a22,a32;
00130       a12=a1*a1;
00131       a22=a2*a2;
00132       a32=a3*a3;
00133 
00134       /* a_i ( K1 (a_j^2 a_k^2) + K2 a_j^2 a_k^2) and permutations */
00135       PetscReal b1,b2,b3;
00136       b1=-2.0*a1*(k1*(a22+a32)+k2*a22*a32);
00137       b2=-2.0*a2*(k1*(a12+a32)+k2*a12*a32);
00138       b3=-2.0*a3*(k1*(a12+a22)+k2*a12*a22);
00139 
00140       /* assemble cubic H components */
00141       ta_Hcubic[ND*j+0] = ax*b1+bx*b2+cx*b3;
00142       ta_Hcubic[ND*j+1] = ay*b1+by*b2+cy*b3;
00143       ta_Hcubic[ND*j+2] = az*b1+bz*b2+cz*b3;
00144 
00145 #ifdef EXCH
00146       Ecubic +=
00147         (k1*(a12*a22+a12*a32+a22*a32)+k2*a12*a22*a32)*ta_VMs3[ND*j+0];
00148     }
00149     else if (gdata->propdat[NP*gdata->vertprop[j]+3] != 0.0) {
00150       PetscReal k2,jsi;
00151       /* might as well do uniaxial k2 term... */
00152 
00153       jsi=gdata->propdat[NP*gdata->vertprop[j]+4];
00154       jsi /= (jsi*jsi+1.e-20);
00155 
00156       /* k2/Js */
00157       k2=gdata->propdat[NP*gdata->vertprop[j]+3]*jsi;
00158 
00159       /* easy axis of uniaxial crystal */
00160       PetscReal ax,ay,az;
00161       ax=gdata->propdat[NP*gdata->vertprop[j]+6];
00162       ay=gdata->propdat[NP*gdata->vertprop[j]+7];
00163       az=gdata->propdat[NP*gdata->vertprop[j]+8];
00164 
00165       PetscReal a1,a2;
00166       a1=ax*ta_M[ND*j+0]+ay*ta_M[ND*j+1]+az*ta_M[ND*j+2];  /* m.k */
00167       /* field magnitude: 4 k2 (m.k) (1-(m.k)^2) / Js */
00168       a2=4.0*k2*a1*(1.0-a1*a1);
00169 
00170       /* energy: k2(1-(m.k)^2)^2 * volume */
00171       Ecubic += k2*(1.0-a1*a1)*(1.0-a1*a1)*ta_VMs3[ND*j+0];
00172 
00173       /* assemble cubic H components */
00174       ta_Hcubic[ND*j+0] = ax*a2;
00175       ta_Hcubic[ND*j+1] = ay*a2;
00176       ta_Hcubic[ND*j+2] = az*a2;
00177 #endif
00178     }
00179     else {
00180       ta_Hcubic[ND*j+0] = 0.0;
00181       ta_Hcubic[ND*j+1] = 0.0;
00182       ta_Hcubic[ND*j+2] = 0.0;
00183     }
00184   }
00185 
00186   ierr = VecRestoreArray(gdata->M,&ta_M);CHKERRQ(ierr);
00187   ierr = VecRestoreArray(VHcubic,&ta_Hcubic);CHKERRQ(ierr);
00188 #ifdef EXCH
00189   ierr = VecRestoreArray(gdata->VMs3,&ta_VMs3);CHKERRQ(ierr);
00190 #endif
00191 
00192   ierr = VecAXPY(VHtotsum,1.0,VHcubic);CHKERRQ(ierr);
00193 
00194   MagparFunctionProfReturn(0);
00195 }
00196 
00197 
00198 int Hcubic_Energy(GridData *gdata,Vec VMom,PetscReal *energy)
00199 {
00200   PetscReal e;
00201 
00202   MagparFunctionInfoBegin;
00203 
00204   if (!fieldon) {
00205     *energy=0.0;
00206     MagparFunctionInfoReturn(0);
00207   }
00208 
00209 #ifdef EXCH
00210   ierr = PetscGlobalSum(&Ecubic,&e,PETSC_COMM_WORLD);CHKERRQ(ierr);
00211 #else
00212   /* only correct if K2=0, otherwise K2 term scaled by -1/6, not -1/4 */
00213   ierr = VecDot(VMom,VHcubic,&e);CHKERRQ(ierr);
00214   e /= -4.0;
00215 #endif
00216 
00217   *energy=e;
00218 
00219   MagparFunctionProfReturn(0);
00220 }
00221 

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