helastic.c

Go to the documentation of this file.
00001 /*
00002     This file is part of magpar.
00003 
00004     Copyright (C) 2006-2009 Ahmet Kaya
00005     Copyright (C) 2005-2009 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: helastic.c 2844 2009-11-05 19:32:09Z 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 PetscReal *rval;
00033 
00034 #define NPE 6
00035 static PetscReal *propdat;
00036 static Vec       VHelastic;
00037 
00038 int ReadProp(int n_prop)
00039 {
00040   PetscTruth   flg;
00041 
00042   MagparFunctionLogBegin;
00043 
00044   int rank,size;
00045   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00046   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00047 
00048   char str[256];
00049   ierr = PetscOptionsGetString(PETSC_NULL,"-helastic_propfile",str,256,&flg);CHKERRQ(ierr);
00050   if (flg!=PETSC_TRUE) {
00051     SETERRQ1(PETSC_ERR_ARG_INCOMP,"Problem with option -helastic_propfile %s\n",str);
00052   }
00053 
00054   ierr = PetscMalloc(n_prop*NPE*sizeof(PetscReal),&propdat);CHKERRQ(ierr);
00055 
00056   if (!rank) {
00057     FILE *fd;
00058     ierr = PetscFOpen(PETSC_COMM_WORLD,str,"r",&fd);CHKERRQ(ierr);
00059     if (!fd) {
00060       SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",str);
00061     }
00062     PetscPrintf(PETSC_COMM_WORLD,"Opening file %s\n",str);
00063 
00064     for (int i=0; i<n_prop; i++) {
00065       for (int j=0; j<NPE; j++) {
00066         int k;
00067 #ifdef PETSC_USE_SINGLE
00068         k=fscanf(fd,"%f",propdat+NPE*i+j);
00069 #else
00070         k=fscanf(fd,"%lf",propdat+NPE*i+j);
00071 #endif
00072         if (k!=1) {
00073           SETERRQ2(PETSC_ERR_FILE_UNEXPECTED,"Data in %s corrupt: %i!\n",str,j);
00074         }
00075       }
00076       /* read remainder of this line */
00077       fgets(str,256,fd);
00078 
00079       PetscPrintf(MPI_COMM_WORLD,"----- property: %i -----\n",i);
00080       PetscPrintf(MPI_COMM_WORLD,
00081         "texture:    %g\n"
00082         "lambda100:  %g\n"
00083         "lambda111:  %g\n"
00084         "sigmaX:     %g\n"
00085         "sigmaY:     %g\n"
00086         "sigmaZ:     %g\n",
00087         propdat[NPE*i+0],
00088         propdat[NPE*i+1],
00089         propdat[NPE*i+2],
00090         propdat[NPE*i+3],
00091         propdat[NPE*i+4],
00092         propdat[NPE*i+5]
00093       );
00094       assert(NPE==6);
00095     }
00096 
00097     ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
00098   }
00099   ierr = MPI_Bcast(propdat,NPE*n_prop,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00100 
00101   MagparFunctionLogReturn(0);
00102 }
00103 
00104 
00105 int Helastic_Init(GridData *gdata)
00106 {
00107   PetscTruth flg;
00108 
00109   MagparFunctionLogBegin;
00110 
00111   ierr = PetscOptionsHasName(PETSC_NULL,"-helastic_propfile",&flg);CHKERRQ(ierr);
00112   if (flg!=PETSC_TRUE) {
00113     PetscPrintf(PETSC_COMM_WORLD,
00114       "Option -helastic_propfile not found, using default value: %i (field off)\n",
00115       fieldon
00116     );
00117   }
00118   else {
00119     fieldon=1;
00120   }
00121 
00122   if (!fieldon) {
00123     PetscPrintf(PETSC_COMM_WORLD,"Helastic off\n");
00124     MagparFunctionLogReturn(0);
00125   }
00126   PetscPrintf(PETSC_COMM_WORLD,"Helastic on\n");
00127 
00128   /* read material properties from file */
00129   ierr = ReadProp(gdata->n_prop);CHKERRQ(ierr);
00130 
00131   ierr = VecDuplicate(gdata->M,&VHelastic);CHKERRQ(ierr);
00132   ierr = VecZeroEntries(VHelastic);CHKERRQ(ierr);
00133 
00134 
00135   /* create arrays of random numbers for random texture */
00136 
00137   ierr = PetscMalloc(ND*gdata->ln_vert*sizeof(PetscReal),&rval);CHKERRQ(ierr);
00138 
00139   PetscRandom r;
00140   ierr = PetscRandomCreate(PETSC_COMM_WORLD,&r);CHKERRQ(ierr);
00141   ierr = PetscRandomSetFromOptions(r);CHKERRQ(ierr);
00142 
00143   for (int i=0;i<gdata->ln_vert;i++) {
00144     PetscReal  value1,value2,value3,value4;
00145     PetscReal  ran_num1,ran_num2,ran_num3;
00146     PetscReal  ww;
00147     PetscReal  xx1,xx2,xx3,xx4;
00148 
00149     do {
00150       PetscRandomGetValue(r,&value1);
00151       PetscRandomGetValue(r,&value2);
00152 
00153       xx1 = 2.0*value1-1.0;
00154       xx2 = 2.0*value2-1.0;
00155       ww = xx1*xx1+xx2*xx2;
00156     } while (ww >= 1.0);
00157 
00158     ww=sqrt((-2.0*log(ww))/ww);
00159     ran_num1=xx1*ww;
00160     ran_num2=xx2*ww;
00161 
00162     do {
00163       PetscRandomGetValue(r,&value3);
00164       PetscRandomGetValue(r,&value4);
00165 
00166       xx3 = 2.0*value3-1.0;
00167       xx4 = 2.0*value4-1.0;
00168       ww = xx3*xx3+xx4*xx4;
00169     } while (ww >= 1.0);
00170 
00171     ww=sqrt((-2.0*log(ww))/ww);
00172     ran_num3 = xx3*ww;
00173 
00174     rval[ND*i+0]=(ran_num1+4.0)/(8.0);
00175     rval[ND*i+1]=(ran_num2+4.0)/(8.0);
00176     rval[ND*i+2]=(ran_num3+4.0)/(8.0);
00177   }
00178 
00179   PetscRandomDestroy(r);
00180   /* end of create random number */
00181 
00182   MagparFunctionLogReturn(0);
00183 }
00184 
00185 
00186 int Helastic(GridData *gdata,Vec VHtotsum)
00187 {
00188   MagparFunctionInfoBegin;
00189 
00190   if (doinit) {
00191     ierr = Helastic_Init(gdata);CHKERRQ(ierr);
00192     doinit=0;
00193   }
00194   if (!fieldon) {
00195     MagparFunctionInfoReturn(0);
00196   }
00197 
00198   PetscReal *ta_Helastic,*ta_M;
00199   ierr = VecGetArray(VHelastic,&ta_Helastic);CHKERRQ(ierr);
00200   ierr = VecGetArray(gdata->M,&ta_M);CHKERRQ(ierr);
00201 
00202   /* TODO: for nodes at interfaces, the assignment is "at will"
00203      not any more: higher Ms wins
00204   */
00205 
00206   for (int j=0; j<gdata->ln_vert; j++) {
00207     int celltexture;
00208     PetscReal xA,xB,xC;
00209 
00210     if (propdat[NPE*gdata->vertprop[j]+0]<=0.0) continue;
00211 
00212     /* FIXME: safer to use round or rint?? */
00213     celltexture=(int)propdat[NPE*gdata->vertprop[j]+0];
00214     switch (celltexture) {
00215       case 1:
00216         xA=0.0;
00217         xB=0.0;
00218         xC=0.0;
00219       break;
00220       case 2:
00221         xA=0.0;
00222         xB=rval[ND*j+0]*2.0*PETSC_PI;
00223         xC=rval[ND*j+1]*2.0*PETSC_PI;
00224       break;
00225       case 3:
00226         xA=PETSC_PI/2.0;
00227         xB=3.0*PETSC_PI/4.0;
00228         xC=0.0;
00229       break;
00230       case 4:
00231         xA=PETSC_PI/2.0;
00232         xB=3.0*PETSC_PI/4.0;
00233         xC=rval[ND*j+0]*2.0*PETSC_PI;
00234       break;
00235       case 5:
00236         xA=PETSC_PI/2.0;
00237         xB=asin(sqrt(2.0/3.0));
00238         xC=PETSC_PI;
00239       break;
00240       case 6:
00241         xA=asin(sqrt(2.0/3.0));
00242         xB=3*PETSC_PI/4.0;
00243         xC=rval[ND*j+0]*2.0*PETSC_PI;
00244       break;
00245       case 7:
00246         xA=rval[ND*j+0]*PETSC_PI;
00247         xB=rval[ND*j+1]*2.0*PETSC_PI;
00248         xC=rval[ND*j+2]*2.0*PETSC_PI;
00249       break;
00250       default:
00251         SETERRQ1(PETSC_ERR_ARG_INCOMP,
00252           "Unknown value for texture: %i\n",
00253           celltexture
00254         );
00255       break;
00256     }
00257 
00258     /* Find the matrix elements from angles */
00259     PetscReal c1,c2,c3;
00260     c1=cos(xA);
00261     c2=cos(xB);
00262     c3=cos(xC);
00263 
00264     PetscReal s1,s2,s3;
00265     s1=sin(xA);
00266     s2=sin(xB);
00267     s3=sin(xC);
00268 
00269     /* matrix A^-1 */
00270     PetscReal a11,a12,a13,a21,a22,a23,a31,a32,a33;
00271 
00272     a11= c2*c3-c1*s2*s3;
00273     a21=-c2*s3-c1*s2*c3;
00274     a31=s1*s2;
00275     a12= s2*c3+c1*c2*s3;
00276     a22=-s2*s3+c1*c2*c3;
00277     a32=-s1*c2;
00278     a13=s1*s3;
00279     a23=s1*c3;
00280     a33=c1;
00281 
00282     /* below operations are for finding inverse Matrix */
00283 
00284     PetscReal detA;
00285     detA=a11*a22*a33+a13*a21*a32+a12*a23*a31-a11*a23*a32-a13*a22*a31-a12*a21*a33;
00286 
00287     /* matrix A */
00288     PetscReal b11,b12,b13,b21,b22,b23,b31,b32,b33;
00289     b11 = (a22*a33 - a23*a32)/detA;
00290     b12 = (a23*a31 - a21*a33)/detA;
00291     b13 = (a21*a32 - a22*a31)/detA;
00292     b21 = (a13*a32 - a12*a33)/detA;
00293     b22 = (a11*a33 - a13*a31)/detA;
00294     b23 = (a12*a31 - a11*a32)/detA;
00295     b31 = (a12*a23 - a13*a22)/detA;
00296     b32 = (a13*a21 - a11*a23)/detA;
00297     b33 = (a11*a22 - a12*a21)/detA;
00298 
00299     /* get components of M */
00300 
00301     PetscReal mx,my,mz;
00302     mx=ta_M[ND*j+0];
00303     my=ta_M[ND*j+1];
00304     mz=ta_M[ND*j+2];
00305 
00306     PetscReal mxp,myp,mzp;
00307     mxp=b11*mx+b21*my+b31*mz;
00308     myp=b12*mx+b22*my+b32*mz;
00309     mzp=b13*mx+b23*my+b33*mz;
00310 
00311     PetscReal lam100,lam111;
00312     lam100=propdat[NPE*gdata->vertprop[j]+1];
00313     lam111=propdat[NPE*gdata->vertprop[j]+2];
00314 
00315     PetscReal sigx,sigy,sigz;
00316     sigx=propdat[NPE*gdata->vertprop[j]+3];
00317     sigy=propdat[NPE*gdata->vertprop[j]+4];
00318     sigz=propdat[NPE*gdata->vertprop[j]+5];
00319 
00320     /* direction cosines of Z stress in xtal coor. sys. */
00321 /*
00322     PetscReal norm;
00323     norm=sqrt(sigx*sigx+sigy*sigy+sigz*sigz);
00324 
00325     PetscReal gamma1g,gamma2g,gamma3g;
00326     gamma1g=sigx/norm;
00327     gamma2g=sigy/norm;
00328     gamma3g=sigz/norm;
00329 
00330     gamma1=b11*gamma1g+b21*gamma2g+b31*gamma3g;
00331     gamma2=b12*gamma1g+b22*gamma2g+b32*gamma3g;
00332     gamma3=b13*gamma1g+b23*gamma2g+b33*gamma3g;
00333 */
00334 
00335     PetscReal msat;
00336     msat=gdata->propdat[NP*gdata->vertprop[j]+4];
00337 
00338     PetscReal gamma1,gamma2,gamma3;
00339     PetscReal elasHkXp,elasHkYp,elasHkZp;
00340 
00341     gamma1=b31;
00342     gamma2=b32;
00343     gamma3=b33;
00344     elasHkXp =  3.0/2.0*lam100*sigz/(msat)*(2.0*mxp*(gamma1*gamma1))+3.0*lam111*sigz/(msat)*(myp*gamma2+mzp*gamma3)*gamma1;
00345     elasHkYp =  3.0/2.0*lam100*sigz/(msat)*(2.0*myp*(gamma2*gamma2))+3.0*lam111*sigz/(msat)*(mzp*gamma3+mxp*gamma1)*gamma2;
00346     elasHkZp =  3.0/2.0*lam100*sigz/(msat)*(2.0*mzp*(gamma3*gamma3))+3.0*lam111*sigz/(msat)*(mxp*gamma1+myp*gamma2)*gamma3;
00347 
00348     gamma1=b21;
00349     gamma2=b22;
00350     gamma3=b23;
00351     elasHkXp += 3.0/2.0*lam100*sigy/(msat)*(2.0*mxp*(gamma1*gamma1))+3.0*lam111*sigy/(msat)*(myp*gamma2+mzp*gamma3)*gamma1;
00352     elasHkYp += 3.0/2.0*lam100*sigy/(msat)*(2.0*myp*(gamma2*gamma2))+3.0*lam111*sigy/(msat)*(mzp*gamma3+mxp*gamma1)*gamma2;
00353     elasHkZp += 3.0/2.0*lam100*sigy/(msat)*(2.0*mzp*(gamma3*gamma3))+3.0*lam111*sigy/(msat)*(mxp*gamma1+myp*gamma2)*gamma3;
00354 
00355     gamma1=b11;
00356     gamma2=b12;
00357     gamma3=b13;
00358     elasHkXp += 3.0/2.0*lam100*sigx/(msat)*(2.0*mxp*(gamma1*gamma1))+3.0*lam111*sigx/(msat)*(myp*gamma2+mzp*gamma3)*gamma1;
00359     elasHkYp += 3.0/2.0*lam100*sigx/(msat)*(2.0*myp*(gamma2*gamma2))+3.0*lam111*sigx/(msat)*(mzp*gamma3+mxp*gamma1)*gamma2;
00360     elasHkZp += 3.0/2.0*lam100*sigx/(msat)*(2.0*mzp*(gamma3*gamma3))+3.0*lam111*sigx/(msat)*(mxp*gamma1+myp*gamma2)*gamma3;
00361 
00362     /* convert into global coordinate system */
00363     PetscReal elasHkX,elasHkY,elasHkZ;
00364     elasHkX=a11*elasHkXp+a12*elasHkYp+a13*elasHkZp;
00365     elasHkY=a21*elasHkXp+a22*elasHkYp+a23*elasHkZp;
00366     elasHkZ=a31*elasHkXp+a32*elasHkYp+a33*elasHkZp;
00367 
00368     /* assemble cubic H components */
00369     ta_Helastic[ND*j+0] = elasHkX/gdata->hscale*MU0;
00370     ta_Helastic[ND*j+1] = elasHkY/gdata->hscale*MU0;
00371     ta_Helastic[ND*j+2] = elasHkZ/gdata->hscale*MU0;
00372   }
00373 
00374   ierr = VecRestoreArray(gdata->M,&ta_M);CHKERRQ(ierr);
00375   ierr = VecRestoreArray(VHelastic,&ta_Helastic);CHKERRQ(ierr);
00376 
00377   ierr = VecAXPY(VHtotsum,1.0,VHelastic);CHKERRQ(ierr);
00378 
00379 /* DEBUG: verify results
00380   PetscPrintf(PETSC_COMM_WORLD,
00381     "deb01: time : %g and elasticfields are.: %lf,%lf ,%lf \n",
00382     gdata->time*gdata->tscale*1e9,elasHkX,elasHkY,elasHkZ
00383   );
00384   ierr = VecView(VHelastic,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
00385 */
00386 
00387   MagparFunctionProfReturn(0);
00388 }
00389 
00390 
00391 int Helastic_Energy(GridData *gdata,Vec VMom,PetscReal *energy)
00392 {
00393   PetscReal e;
00394 
00395   MagparFunctionInfoBegin;
00396 
00397   if (!fieldon) {
00398     *energy=0.0;
00399     MagparFunctionInfoReturn(0);
00400   }
00401 
00402   ierr = VecDot(VMom,VHelastic,&e);CHKERRQ(ierr);
00403   e *= -1.0;
00404 
00405   *energy=e;
00406 
00407   MagparFunctionProfReturn(0);
00408 }
00409 

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