htot.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: htot.c 2966 2010-02-08 23:31:15Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "field.h"

#ifdef ADDONS
#include "addons/addons.h"
#endif

#ifdef SPINTORQ
#include "spintorq/spintorq.h"
#endif

static int doinit=1;
static Vec VHtotsum=PETSC_NULL;
static Vec VMom; 
int Htot_Init(GridData *gdata)
{
  MagparFunctionLogBegin;

  PetscPrintf(PETSC_COMM_WORLD,"Htot on\n");

  ierr = VecDuplicate(gdata->M,&gdata->VHtot);CHKERRQ(ierr);
  ierr = VecDuplicate(gdata->M,&VMom);CHKERRQ(ierr);
  ierr = VecDuplicate(gdata->M,&VHtotsum);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}


int Htot(GridData *gdata)
{
  MagparFunctionInfoBegin;

  /* initialize if necessary */
  if (doinit) {
    ierr = Htot_Init(gdata);CHKERRQ(ierr);
    doinit=0;
  }

  ierr = VecZeroEntries(VHtotsum);CHKERRQ(ierr);

  ierr = Hcubic(gdata,VHtotsum);CHKERRQ(ierr);
  ierr = Hdemag(gdata,VHtotsum);CHKERRQ(ierr);
  ierr = Helastic(gdata,VHtotsum);CHKERRQ(ierr);
  ierr = Hexchani(gdata,VHtotsum);CHKERRQ(ierr);
  ierr = Hexternal(gdata,VHtotsum);CHKERRQ(ierr);
#ifdef ADDONS
  ierr = HIlayerExch(gdata,VHtotsum);CHKERRQ(ierr);
  ierr = Hnonlin(gdata,VHtotsum);CHKERRQ(ierr);
  ierr = Hcancel(gdata,VHtotsum);CHKERRQ(ierr);
#endif
#ifdef SPINTORQ
  /* initialize vectors to NULL */
  gdata->VHspintorque=PETSC_NULL;
  gdata->VHspinforce=PETSC_NULL;

  ierr = Hspintorq(gdata);CHKERRQ(ierr);
  ierr = Hspinforce(gdata);CHKERRQ(ierr);

  if (gdata->VHspintorque!=PETSC_NULL && gdata->VHspinforce!=PETSC_NULL) {
    SETERRQ(PETSC_ERR_ARG_INCOMP,"Spin torque and spin force are switched on at the same time!\n");
  }
#endif

  ierr = VecCopy(VHtotsum,gdata->VHtot);CHKERRQ(ierr);

  MagparFunctionProfReturn(0);
}


int Htot_Gradient(GridData *gdata)
{
  Vec VHtemp;

  MagparFunctionInfoBegin;

  /* abusing VMom */
  VHtemp=VMom;

  /* initialize if necessary */
  if (doinit) {
    ierr = Htot_Init(gdata);CHKERRQ(ierr);
    doinit=0;
  }

#if 1
  /* gives correct result for Etot in example mumag3: 1.240375e+05 ; LLG: 1.240388e+05 */
  /* appears to converge faster for example fept, headfield_nonlin w/o demag (but higher Etot) ! */
  /*
  m=M*Ms*V
  grad(E)=dE/dm=dE/d(M*Ms*V)
  H=-dE/dm
  grad(E)=-H

  However, Ms and V are fixed and only M varies in our formulation:
  grad(E)=dE/dM=Ms*V*dE/d(M*Ms*V)=Ms*V*dE/dm=Ms*V*(-H)
  */

  ierr = Htot(gdata);CHKERRQ(ierr);

  /* This scaling with gdata->VMs3 is necessary for big writer models to converge properly.
     Causes wiggles in M in headfield_nonlin while iterating, but converges fine.
  */
  ierr = VecPointwiseMult(gdata->VHtot,gdata->VHtot,gdata->VMs3);CHKERRQ(ierr);
#else
  /* for comparison between hand-coded and TAO's finite difference gradient: magpar_examples5/cube_test_grad */
  /* gives incorrect result for Etot in example mumag3: 1.256195e+05 (too high) */

  ierr = VecZeroEntries(VHtotsum);CHKERRQ(ierr);

  /* calculate field/gradient contributions and apply fudge (?) factor */
  VecZeroEntries(VHtemp);Hcubic(gdata,VHtemp);     VecAXPY(VHtotsum,1.0,VHtemp); /* ok */
  VecZeroEntries(VHtemp);Hdemag(gdata,VHtemp);     VecAXPY(VHtotsum,0.5,VHtemp); /* ok */
  VecZeroEntries(VHtemp);Helastic(gdata,VHtemp);   VecAXPY(VHtotsum,2.0,VHtemp); /* ok */
  VecZeroEntries(VHtemp);Hexchani(gdata,VHtemp);   VecAXPY(VHtotsum,1.0,VHtemp); /* ok */
  VecZeroEntries(VHtemp);Hexternal(gdata,VHtemp);  VecAXPY(VHtotsum,1.0,VHtemp); /* ok */

#ifdef ADDONS
  VecZeroEntries(VHtemp);HIlayerExch(gdata,VHtemp);VecAXPY(VHtotsum,1.0,VHtemp); /* ok */
  VecZeroEntries(VHtemp);Hnonlin(gdata,VHtemp);    VecAXPY(VHtotsum,1.0,VHtemp); /* ok */
  VecZeroEntries(VHtemp);Hcancel(gdata,VHtemp);    VecAXPY(VHtotsum,1.0,VHtemp); /* TODO */
#endif

  ierr = VecCopy(VHtotsum,gdata->VHtot);CHKERRQ(ierr);
  ierr = VecPointwiseMult(gdata->VHtot,gdata->VHtot,gdata->VMs3);CHKERRQ(ierr);
#endif

  ierr = VecScale(gdata->VHtot,-1.0);CHKERRQ(ierr);

  PetscReal *ta_Htot;
  VecGetArray(gdata->VHtot,&ta_Htot);

  /* TODO: set up Vec and do VecPointwiseMult instead of the loop */
  for (int i=0;i<gdata->ln_vert;i++) {
    /* skip non-magnetic vertices and vertices with alpha >= RIGID_M_ALPHA */
    if (gdata->propdat[NP*gdata->vertprop[i]+4]<=0.0 ||
        gdata->propdat[NP*gdata->vertprop[i]+9]>=RIGID_M_ALPHA) {
      ta_Htot[ND*i+0]=0.0;
      ta_Htot[ND*i+1]=0.0;
      ta_Htot[ND*i+2]=0.0;
      continue;
    }
  }

  VecRestoreArray(gdata->VHtot,&ta_Htot);

  MagparFunctionProfReturn(0);
}


int CalcEnergy(Vec M,GridData *gdata)
{
  MagparFunctionInfoBegin;

  /* calculate "magnetic moments" VMom on vertices */
  ierr = VecPointwiseMult(VMom,gdata->VMs3,gdata->M);CHKERRQ(ierr);

  PetscReal e,etot;
  etot=0;

  ierr = Hcubic_Energy(gdata,M,&e);     CHKERRQ(ierr);etot += e;gdata->Eexchani  = e/gdata->totvol; /* part of anisotropy energy */
  ierr = Hdemag_Energy(gdata,M,&e);     CHKERRQ(ierr);etot += e;gdata->Edem      = e/gdata->totvol;
  ierr = Helastic_Energy(gdata,M,&e);   CHKERRQ(ierr);etot += e;
  ierr = Hexchani_Energy(gdata,M,&e);   CHKERRQ(ierr);etot += e;gdata->Eexchani += e/gdata->totvol; /* add to anisotropy+exchange energy */
#ifdef EXCH
  ierr = Hexchonly_Energy(gdata,M,&e);  CHKERRQ(ierr);          gdata->Eexchonly = e/gdata->totvol; /* don't add to etot!!! */
#endif
  ierr = Hexternal_Energy(gdata,M,&e);  CHKERRQ(ierr);etot += e;gdata->Eext      = e/gdata->totvol;
#ifdef ADDONS
  ierr = Hcancel_Energy(gdata,M,&e);    CHKERRQ(ierr);etot += e;gdata->Eext     += e/gdata->totvol;
  ierr = HIlayerExch_Energy(gdata,M,&e);CHKERRQ(ierr);etot += e;gdata->Eexchani += e/gdata->totvol; /* add to anisotropy+exchange energy */
  ierr = Hnonlin_Energy(gdata,M,&e);    CHKERRQ(ierr);etot += e;gdata->Eexchani += e/gdata->totvol; /* add to anisotropy+exchange energy */
#ifdef EXCH
                                                                gdata->Eexchonly+= e/gdata->totvol; /* add to exchange only energy */
#endif
#endif
#ifdef SPINTORQ
/* TODO: implement
  ierr = Hspintorq_Energy(gdata,M,&e);  CHKERRQ(ierr);etot += e;
  ierr = Hspinforce_Energy(gdata,M,&e); CHKERRQ(ierr);etot += e;
*/
#endif

  gdata->Etot=etot/gdata->totvol;

  MagparFunctionProfReturn(0);
}


int Htot_Energy(GridData *gdata)
{
  MagparFunctionInfoBegin;

  /* calculate "magnetic moments" VMom on vertices */
  ierr = VecPointwiseMult(VMom,gdata->VMs3,gdata->M);CHKERRQ(ierr);

  ierr = CalcEnergy(VMom,gdata);CHKERRQ(ierr);

  MagparFunctionProfReturn(0);
}


int Htot_EminiEnergy(GridData *gdata, PetscReal *f)
{
  MagparFunctionInfoBegin;

#ifdef NO_VMS3
  /* calculate "pseudo" energy for energy minmizer */
  ierr = CalcEnergy(gdata->M,gdata);CHKERRQ(ierr);
  *f = gdata->Etot*gdata->totvol;
#endif

  /* calculate real energy with "magnetic moments" VMom on vertices */
  ierr = VecPointwiseMult(VMom,gdata->VMs3,gdata->M);CHKERRQ(ierr);
  ierr = CalcEnergy(VMom,gdata);CHKERRQ(ierr);

#ifndef NO_VMS3
#endif
  if (gdata->mode==53) {
    *f = gdata->Etot*gdata->totvol;
  }

  PetscPrintf(PETSC_COMM_WORLD,"debug Etot=%g\n",*f);
/*
  gdata->inp=PetscAbs(gdata->inp);
  ierr = WriteLog(gdata);CHKERRQ(ierr);
  ierr = WriteSet(gdata);CHKERRQ(ierr);
*/

  MagparFunctionProfReturn(0);
}