hexch_ani.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2010 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: hexch_ani.c 3023 2010-03-27 23:42:04Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "field.h"
#include "init/init.h"
#include "util/util.h"

static int       doinit=1;
static int       fieldon=0;

       Mat       Ad2E_dMidMk; 
static PetscReal Eanioffset;  
#ifdef EXCH
static Mat       AHexchonly;  
#endif

#define PREALLOC_DG 15
#define PREALLOC_OD 5

int Hexchani_Init(GridData *gdata)
{
  MagparFunctionLogBegin;

  int rank,size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  fieldon=0;

  Ad2E_dMidMk=PETSC_NULL;
  /* create gdata->VHexchani in any case
     otherwise we get SEGV in WriteAVS
     TODO: not most elegant way to solve this issue
  */
  ierr = VecDuplicate(gdata->M,&gdata->VHexchani);CHKERRQ(ierr);

  /* activate only if Ms>0 && (K1!=0 or A!=0) */
  for (int i=0;i<gdata->n_prop;i++) {
    if (gdata->propdat[NP*i+4]>0.0 && (gdata->propdat[NP*i+2]!=0.0 || gdata->propdat[NP*i+5]!=0.0)) {
      fieldon=1;
      break;
    }
  }
  if (!fieldon) {
    PetscPrintf(PETSC_COMM_WORLD,"Hexchani off\n");
    ierr = PetscOptionsSetValue("-precon","0");CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD,"Disabling preconditioning\n");

    gdata->VHexchani=PETSC_NULL;
#ifdef EXCH
    gdata->VHexchonly=PETSC_NULL;
#endif

    MagparFunctionLogReturn(0);
  }
  PetscPrintf(PETSC_COMM_WORLD,"Hexchani on\n");

  PetscPrintf(PETSC_COMM_WORLD,"Calculating sparsity pattern for matrix preallocation\n");

  /* get ownership range */
  int low,high;
  ierr = VecGetOwnershipRange(gdata->elevol,&low,&high);CHKERRQ(ierr);

  int *elevertall;

#ifdef HAVE_ELEVERTALL
  elevertall=gdata->elevertall;
#else
  if (size>1) {
    ierr = PetscMalloc(NV*gdata->n_ele*sizeof(int),&elevertall);CHKERRQ(ierr);
    ierr = PetscMemcpy(elevertall+NV*low,gdata->elevert,NV*gdata->ln_ele*sizeof(int));CHKERRQ(ierr);

    int *p;
    p=elevertall;

    int recvcount[size];
    recvcount[rank]=NV*gdata->ln_ele;
    for (int i=0;i<size;i++) {
      ierr = MPI_Bcast(recvcount+i,1,MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
      /* we could also use MPI_Send/Recv to get everything just to the first processor */
      ierr = MPI_Bcast(p,recvcount[i],MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
      p+=recvcount[i];
    }
    assert(p==elevertall+NV*gdata->n_ele);
  }
  else {
    elevertall=gdata->elevert;
  }
#endif

  int *ia,*ja;
#ifdef HAVE_M2IJ
  ia=gdata->mesh2nodal_ia;
  ja=gdata->mesh2nodal_ja;
#else
  ierr = Mesh2Nodal(gdata->n_ele,gdata->n_vert,elevertall,&ia,&ja);CHKERRQ(ierr);
#endif

  /* get ownership range */
  ierr = VecGetOwnershipRange(gdata->M,&low,&high);CHKERRQ(ierr);
  low/=ND;
  high/=ND;

  int *d_nnz;
  int *o_nnz;
  ierr = PetscMalloc(ND*(high-low)*sizeof(int),&d_nnz);CHKERRQ(ierr);
  ierr = PetscMalloc(ND*(high-low)*sizeof(int),&o_nnz);CHKERRQ(ierr);
  /* calculate number of nonzeroes for each line */
  for (int i=low;i<high;i++) {
    for (int j=0;j<ND;j++) {
      d_nnz[ND*(i-low)+j]=ND;
      o_nnz[ND*(i-low)+j]=0;
      for (int l=ia[i];l<ia[i+1];l++) {
        if (ja[l]>=low && ja[l]<high) {
          d_nnz[ND*(i-low)+j]+=ND;
        }
        else {
          o_nnz[ND*(i-low)+j]+=ND;
        }
      }
    }
  }

  PetscPrintf(PETSC_COMM_WORLD,"Creating matrices\n");
  /* create matrix to calculate Ad2E_dMidMk */
/* normal MPIAIJ format: */
  ierr = MatCreateMPIAIJ(
    PETSC_COMM_WORLD,
    ND*gdata->ln_vert,ND*gdata->ln_vert,
    ND*gdata->n_vert,ND*gdata->n_vert,
    0,d_nnz,0,o_nnz,
    &Ad2E_dMidMk
  );CHKERRQ(ierr);
  ierr = MatSetFromOptions(Ad2E_dMidMk);CHKERRQ(ierr);
/* block MPIAIJ format:
  ierr = MatCreateMPIBAIJ(
    PETSC_COMM_WORLD,
    ND,
    ND*gdata->ln_vert,ND*gdata->ln_vert,
    ND*gdata->n_vert,ND*gdata->n_vert,
    PREALLOC_DG,PETSC_NULL,PREALLOC_OD,PETSC_NULL,
    &Ad2E_dMidMk
  );CHKERRQ(ierr);
  ierr = MatSetFromOptions(Ad2E_dMidMk);CHKERRQ(ierr);
*/

#ifdef EXCH
  /* calculate number of nonzeroes for each line */
  for (int i=low;i<high;i++) {
    int k;
    k=ia[i+1]-ia[i]+1;
    for (int j=0;j<ND;j++) {
      if (i>=low || i<high) {
        d_nnz[ND*(i-low)+j]=k;
      }
      else {
        o_nnz[ND*(i-low)+j]=k;
      }
    }
  }

  /* matrix for exchange field */
  ierr = MatCreateMPIAIJ(
    PETSC_COMM_WORLD,
    ND*gdata->ln_vert,ND*gdata->ln_vert,
    ND*gdata->n_vert,ND*gdata->n_vert,
    0,d_nnz,0,o_nnz,
    &AHexchonly
  );CHKERRQ(ierr);
  ierr = MatSetFromOptions(AHexchonly);CHKERRQ(ierr);

  ierr = VecDuplicate(gdata->M,&gdata->VHexchonly);CHKERRQ(ierr);
#endif

#ifndef HAVE_M2IJ
  ierr = PetscFree(ia);CHKERRQ(ierr);
  ierr = PetscFree(ja);CHKERRQ(ierr);
#endif
  ierr = PetscFree(d_nnz);CHKERRQ(ierr);
  ierr = PetscFree(o_nnz);CHKERRQ(ierr);
#ifndef HAVE_ELEVERTALL
  if (size>1) {
    ierr = PetscFree(elevertall);CHKERRQ(ierr);
  }
#endif

  PetscPrintf(PETSC_COMM_WORLD,"Calculating matrix elements\n");

  PetscReal *ta_elevol;
  ierr = VecGetArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);

  for (int i=0; i<gdata->ln_ele; i++) {
    /* skip elements with Ms<=0 */
    if (gdata->propdat[NP*gdata->eleprop[i]+4]<=0.0) {
      continue;
    }

    /* get element volume */
    PetscReal t_elevol;
    t_elevol=ta_elevol[i];

    /* get vertex ids */
    int *t_v;
    t_v=gdata->elevert+i*NV;

    /* Get Vertex coordinates */
    PetscReal *x1,*x2,*x3,*x4;
    x1=gdata->vertxyz+(t_v[0])*ND;
    x2=gdata->vertxyz+(t_v[1])*ND;
    x3=gdata->vertxyz+(t_v[2])*ND;
    x4=gdata->vertxyz+(t_v[3])*ND;

    PetscReal D_etaj[NV][ND];
    tetgrad(x1,x2,x3,x4,D_etaj);

    /* set matrix elements */

    /* prefactor for exchange field: pre_exch=2*A */
    PetscReal   pre_exch;
    pre_exch=2.0*gdata->propdat[NP*gdata->eleprop[i]+5];
/* must not divide by Js to account for different materials:
   |M|=1 -> Js has to be considered in local fields !!!
      /gdata->propdat[NP*gdata->eleprop[i]+4];
*/
    for (int l=0;l<NV;l++) {
      for (int j=0;j<NV;j++) {
        PetscReal matele;

        /* set matrix element for jacobian */
        matele=-t_elevol*my_ddot(ND,D_etaj[j],1,D_etaj[l],1)*pre_exch;
        if (PetscAbsReal(matele)>D_EPS) {
          for (int k=0;k<ND;k++) {
            ierr = MatSetValue(
              Ad2E_dMidMk,
              ND*t_v[j]+k,
              ND*t_v[l]+k,
              matele,
              ADD_VALUES
            );CHKERRQ(ierr);
#ifdef EXCH
            ierr = MatSetValue(
              AHexchonly,
              ND*t_v[j]+k,
              ND*t_v[l]+k,
              matele,
              ADD_VALUES
            );CHKERRQ(ierr);
#endif
          }
        }
      }
    }
  }

  ierr = MatAssemblyBegin(Ad2E_dMidMk,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Ad2E_dMidMk,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);

#ifdef EXCH
  ierr = MatAssemblyBegin(AHexchonly,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(AHexchonly,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
#endif

  ierr = PetscGetTime(&t_t2);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,
    "exchange matrix calculation took %g s = %g min\n",
    t_t2-t_t1,
    (t_t2-t_t1)/60.0
  );


  Eanioffset = 0.0;

  PetscReal   pre_ani;
  for (int i=0;i<gdata->ln_ele;i++) {
    /* skip elements with Ms<=0 */
    if (gdata->propdat[NP*gdata->eleprop[i]+4]<=0.0) {
      continue;
    }
    /* skip elements with K1==0 */
    if (PetscAbsReal(gdata->propdat[NP*gdata->eleprop[i]+2])<D_EPS) {
      continue;
    }
    /* skip matrix element if material has cubic anisotropy */
    if (gdata->propdat[NP*gdata->eleprop[i]+16]) {
      continue;
    }

    /* prefactor for anisotropy field: pre_ani=2*K1*V */
    pre_ani=2.0*ta_elevol[i]*gdata->propdat[NP*gdata->eleprop[i]+2];
    /* division by Js to account for different materials done below (vertMomrec3) */

    Eanioffset += pre_ani/2.0;

    /* TODO: more elegant with MatSetValuesLocal !? */
    for (int k=0;k<ND;k++) {
      for (int l=0;l<ND;l++) {
        PetscReal matele;

        matele=
          pre_ani*
          gdata->propdat[NP*gdata->eleprop[i]+6+k]*
          gdata->propdat[NP*gdata->eleprop[i]+6+l];

        if (PetscAbsReal(matele)<D_EPS) continue;

        for (int m=0;m<NV;m++) {
          /* calculate Int(phi[i]*phi[j])

             i==j ? Int(phi[i]*phi[j])=V/10 : Int(phi[i]*phi[j])=V/20
             source: Advanced Finite Element Methods (ASEN 6367) - Spring 2009
                     Aerospace Engineering Sciences - University of Colorado at Boulder
                     Chapter 16: THE LINEAR TETRAHEDRON
                     page 10, Eq. (16.35)
                     http://www.colorado.edu/engineering/cas/courses.d/AFEM.d/AFEM.Ch16.d/AFEM.Ch16.pdf

                     http://smirnov.mae.wvu.edu/hedra/poisson.pdf

                     Wenjie Chen, D. R. Fredkin, and T. R. Koehler,
                     "A New Finite Element Method in Micromagnetics,"
                     IEEE Trans. Magn. 29 (1993) 2124.
                     http://ieeexplore.ieee.org/iel1/20/5774/00221033.pdf

                     see also hdemag.c
          */
          for (int j=0;j<NV;j++) {
            PetscReal matele2;

            if (j==m) {
              matele2=matele/10.0;
            }
            else {
              matele2=matele/20.0;
            }

            ierr = MatSetValue(
              Ad2E_dMidMk,
              gdata->elevert[NV*i+j]*ND+k,
              gdata->elevert[NV*i+m]*ND+l,
              matele2,
              ADD_VALUES
            );CHKERRQ(ierr);
          }
        }
      }
    }
  }

  ierr = MatAssemblyBegin(Ad2E_dMidMk,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Ad2E_dMidMk,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);

  PetscLogDouble t_t3;
  ierr = PetscGetTime(&t_t3);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,
    "anisotropy matrix calculation took %g s = %g min\n",
    t_t3-t_t2,
    (t_t3-t_t2)/60.0
  );

  ierr = VecRestoreArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);

  PetscReal tmp;
  ierr = PetscGlobalSum(&Eanioffset,&tmp,PETSC_COMM_WORLD);CHKERRQ(ierr);
  Eanioffset=tmp;

  PetscPrintf(PETSC_COMM_WORLD,"Eanioffset: %g = %g J = %g J/m^3\n",
    Eanioffset,
    Eanioffset*gdata->escale,
    Eanioffset*gdata->escale/gdata->totvol
  );

  Vec vertMomrec3;
  ierr = VecDuplicate(gdata->M,&vertMomrec3);CHKERRQ(ierr);
  ierr = VecCopy(gdata->VMs3,vertMomrec3);CHKERRQ(ierr);
  ierr = VecReciprocal(vertMomrec3);CHKERRQ(ierr);

  PetscInfo(0,"Ad2E_dMidMk matrix assembly complete\n");
  ierr = MatAssemblyBegin(Ad2E_dMidMk,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Ad2E_dMidMk,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  ierr = MatDiagonalScale(Ad2E_dMidMk,vertMomrec3,PETSC_NULL);CHKERRQ(ierr);
  ierr = PrintMatInfoAll(Ad2E_dMidMk);CHKERRQ(ierr);

#ifdef EXCH
  PetscInfo(0,"AHexchonly matrix assembly complete\n");
  ierr = MatAssemblyBegin(AHexchonly,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(AHexchonly,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatDiagonalScale(AHexchonly,vertMomrec3,PETSC_NULL);CHKERRQ(ierr);
  ierr = PrintMatInfoAll(AHexchonly);CHKERRQ(ierr);
#endif

  ierr = VecDestroy(vertMomrec3);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}


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

  if (doinit) {
    ierr = Hexchani_Init(gdata);CHKERRQ(ierr);
    doinit=0;
  }
  if (!fieldon) {
    MagparFunctionInfoReturn(0);
  }

  /* calculate field */
  ierr = MatMult(Ad2E_dMidMk,gdata->M,gdata->VHexchani);CHKERRQ(ierr);

  ierr = VecAXPY(VHtotsum,1.0,gdata->VHexchani);CHKERRQ(ierr);

#ifdef EXCH
  /* calculate exchange field separately */
  ierr = MatMult(AHexchonly,gdata->M,gdata->VHexchonly);CHKERRQ(ierr);
#endif

  MagparFunctionProfReturn(0);
}


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

  MagparFunctionInfoBegin;

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

  /* calculate exchange+anisotropy energy */
  ierr = VecDot(VMom,gdata->VHexchani,&e);CHKERRQ(ierr);
  e /= -2.0;
  e += Eanioffset;

  *energy=e;

  MagparFunctionProfReturn(0);
}


#ifdef EXCH
int Hexchonly_Energy(GridData *gdata,Vec VMom,PetscReal *energy)
{
  PetscReal e;

  MagparFunctionInfoBegin;

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

  /* calculate exchange energy */
  ierr = VecDot(VMom,gdata->VHexchonly,&e);CHKERRQ(ierr);
  e /= -2.0;

  *energy=e;

  MagparFunctionProfReturn(0);
}
#endif