hcubic.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2005-2009 Greg Parker, Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: hcubic.c 2681 2009-07-31 04:30:53Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "field.h"

static int       doinit=1;
static int       fieldon=0;
static Vec       VHcubic;
#ifdef EXCH
static PetscReal Ecubic;
#endif

int Hcubic_Init(GridData *gdata)
{
  MagparFunctionLogBegin;

  for (int i=0;i<gdata->n_prop;i++) {
    if (gdata->propdat[NP*i+16]==1 ||
        PetscAbsReal(gdata->propdat[NP*i+ 3]) > D_EPS) {
      /* FIXME: latter case not really cubic, but hcubic.c calculates the K_2 term
         anisotropy field calculation then off!? Hani calculated twice?
      */
      fieldon=1;
      break;
    }
  }

  if (fieldon) {
    PetscPrintf(PETSC_COMM_WORLD,"Hcubic on\n");
    ierr = VecDuplicate(gdata->M,&VHcubic);CHKERRQ(ierr);
  }
  else {
    PetscPrintf(PETSC_COMM_WORLD,"Hcubic off\n");
  }

  MagparFunctionLogReturn(0);
}


int Hcubic(GridData *gdata,Vec VHtotsum)
{
  MagparFunctionInfoBegin;

  /* initialize if necessary */
  if (doinit) {
    ierr = Hcubic_Init(gdata);CHKERRQ(ierr);
    doinit=0;
  }
  if (!fieldon) {
    MagparFunctionInfoReturn(0);
  }

  PetscReal *ta_Hcubic,*ta_M;
  ierr = VecGetArray(VHcubic,&ta_Hcubic);CHKERRQ(ierr);
  ierr = VecGetArray(gdata->M,&ta_M);CHKERRQ(ierr);

#ifdef EXCH
  PetscReal *ta_VMs3;
  ierr = VecGetArray(gdata->VMs3,&ta_VMs3);CHKERRQ(ierr);
  Ecubic=0.0;
#endif

  /* TODO: for nodes at interfaces, the assignment is "at will" */
  for (int j=0; j<gdata->ln_vert; j++) {
    if (gdata->propdat[NP*gdata->vertprop[j]+16]) {

      /* only collect cubic Hk if 'grain' has cubic Hk */
      PetscReal k1,k2,jsi;
      jsi=gdata->propdat[NP*gdata->vertprop[j]+4];
      /* FIXME: why 1e-20? */
      jsi /= (jsi*jsi+1.e-20);
      k1=gdata->propdat[NP*gdata->vertprop[j]+2]*jsi;
      k2=gdata->propdat[NP*gdata->vertprop[j]+3]*jsi;

      /* a-axis of cubic crystal */
      PetscReal ax,ay,az;
      ax=gdata->propdat[NP*gdata->vertprop[j]+6];
      ay=gdata->propdat[NP*gdata->vertprop[j]+7];
      az=gdata->propdat[NP*gdata->vertprop[j]+8];

      /* b-axis of cubic crystal */
      PetscReal bx,by,bz;
      bx=gdata->propdat[NP*gdata->vertprop[j]+10];
      by=gdata->propdat[NP*gdata->vertprop[j]+11];
      bz=gdata->propdat[NP*gdata->vertprop[j]+12];

      /* c-axis of cubic crystal */
      PetscReal cx,cy,cz;
      cx=gdata->propdat[NP*gdata->vertprop[j]+13];
      cy=gdata->propdat[NP*gdata->vertprop[j]+14];
      cz=gdata->propdat[NP*gdata->vertprop[j]+15];

      /* get components of M */
      PetscReal mx,my,mz;
      mx=ta_M[ND*j+0];
      my=ta_M[ND*j+1];
      mz=ta_M[ND*j+2];

      /* direction cosines and their squares */
      PetscReal a1,a2,a3;
      a1=mx*ax+my*ay+mz*az;
      a2=mx*bx+my*by+mz*bz;
      a3=mx*cx+my*cy+mz*cz;

      PetscReal a12,a22,a32;
      a12=a1*a1;
      a22=a2*a2;
      a32=a3*a3;

      /* a_i ( K1 (a_j^2 a_k^2) + K2 a_j^2 a_k^2) and permutations */
      PetscReal b1,b2,b3;
      b1=-2.0*a1*(k1*(a22+a32)+k2*a22*a32);
      b2=-2.0*a2*(k1*(a12+a32)+k2*a12*a32);
      b3=-2.0*a3*(k1*(a12+a22)+k2*a12*a22);

      /* assemble cubic H components */
      ta_Hcubic[ND*j+0] = ax*b1+bx*b2+cx*b3;
      ta_Hcubic[ND*j+1] = ay*b1+by*b2+cy*b3;
      ta_Hcubic[ND*j+2] = az*b1+bz*b2+cz*b3;

#ifdef EXCH
      Ecubic +=
        (k1*(a12*a22+a12*a32+a22*a32)+k2*a12*a22*a32)*ta_VMs3[ND*j+0];
    }
    else if (gdata->propdat[NP*gdata->vertprop[j]+3] != 0.0) {
      PetscReal k2,jsi;
      /* might as well do uniaxial k2 term... */

      jsi=gdata->propdat[NP*gdata->vertprop[j]+4];
      jsi /= (jsi*jsi+1.e-20);

      /* k2/Js */
      k2=gdata->propdat[NP*gdata->vertprop[j]+3]*jsi;

      /* easy axis of uniaxial crystal */
      PetscReal ax,ay,az;
      ax=gdata->propdat[NP*gdata->vertprop[j]+6];
      ay=gdata->propdat[NP*gdata->vertprop[j]+7];
      az=gdata->propdat[NP*gdata->vertprop[j]+8];

      PetscReal a1,a2;
      a1=ax*ta_M[ND*j+0]+ay*ta_M[ND*j+1]+az*ta_M[ND*j+2];  /* m.k */
      /* field magnitude: 4 k2 (m.k) (1-(m.k)^2) / Js */
      a2=4.0*k2*a1*(1.0-a1*a1);

      /* energy: k2(1-(m.k)^2)^2 * volume */
      Ecubic += k2*(1.0-a1*a1)*(1.0-a1*a1)*ta_VMs3[ND*j+0];

      /* assemble cubic H components */
      ta_Hcubic[ND*j+0] = ax*a2;
      ta_Hcubic[ND*j+1] = ay*a2;
      ta_Hcubic[ND*j+2] = az*a2;
#endif
    }
    else {
      ta_Hcubic[ND*j+0] = 0.0;
      ta_Hcubic[ND*j+1] = 0.0;
      ta_Hcubic[ND*j+2] = 0.0;
    }
  }

  ierr = VecRestoreArray(gdata->M,&ta_M);CHKERRQ(ierr);
  ierr = VecRestoreArray(VHcubic,&ta_Hcubic);CHKERRQ(ierr);
#ifdef EXCH
  ierr = VecRestoreArray(gdata->VMs3,&ta_VMs3);CHKERRQ(ierr);
#endif

  ierr = VecAXPY(VHtotsum,1.0,VHcubic);CHKERRQ(ierr);

  MagparFunctionProfReturn(0);
}


int Hcubic_Energy(GridData *gdata,Vec VMom,PetscReal *energy)
{
  PetscReal e;

  MagparFunctionInfoBegin;

  if (!fieldon) {
    *energy=0.0;
    MagparFunctionInfoReturn(0);
  }

#ifdef EXCH
  ierr = PetscGlobalSum(&Ecubic,&e,PETSC_COMM_WORLD);CHKERRQ(ierr);
#else
  /* only correct if K2=0, otherwise K2 term scaled by -1/6, not -1/4 */
  ierr = VecDot(VMom,VHcubic,&e);CHKERRQ(ierr);
  e /= -4.0;
#endif

  *energy=e;

  MagparFunctionProfReturn(0);
}