magset.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2009 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: magset.c 2753 2009-08-31 14:49:58Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "init.h"
#include "io/magpario.h"
#include "util/util.h"

int MagSet(int pid,int init_mag,PetscReal init_magparm,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);


  PetscPrintf(PETSC_COMM_WORLD,"pid:          %i\n",pid+1);
  PetscPrintf(PETSC_COMM_WORLD,"init_mag:     %i\n",init_mag);
  PetscPrintf(PETSC_COMM_WORLD,"init_magparm: %g\n",init_magparm);

  if (init_mag==0) {
    char fmesh[256];
    PetscPrintf(PETSC_COMM_WORLD,"inp file to be read: %s.%04d.inp\n",gdata->simname,gdata->inp);
    ierr = PetscSNPrintf(fmesh,255,"%s.%04d.inp",gdata->simname,gdata->inp);CHKERRQ(ierr);
    gdata->inp++;

    ierr = ReadINP(gdata,fmesh,gdata->M,1);CHKERRQ(ierr);
    MagparFunctionLogReturn(0);
  }
  else if (PetscAbsInt(init_mag)==14) {
#ifdef PYTHON
    if (!rank) {
      /* Python stuff */

      if(!gdata->pModule) {
        SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Cannot execute Python function. Module %s has not been imported!\n",MAGPARSTDPYTHONMOD);
      }

#define INITMAG "initmag"
      PyObject *pFunc;
      pFunc=PyObject_GetAttrString(gdata->pModule,INITMAG);  /* new reference */
      if(!pFunc) {
        SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Cannot find Python function %s!\n",INITMAG);
      }
      if (!PyCallable_Check(pFunc)) {
        SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Python object %s is not callable (not a function)!\n",INITMAG);
      }

      /* Arguments for the module */
      PyObject *pTuple;
      pTuple = PyTuple_New(ND); /* new reference */

      PyObject *pVal=PETSC_NULL;

      /* this example sets the external field in a certain area to a
         non-zero value and everywhere else to zero
      */

      /* make the magnetization vector accessible */
      PetscReal *ta_M;
      ierr = VecGetArray(gdata->M,&ta_M);CHKERRQ(ierr);

      /* loop over all nodes */
      for (int j=0; j<gdata->n_vert; j++) {
        PetscReal x,y,z;
        x=gdata->vertxyz[ND*j+0];
        y=gdata->vertxyz[ND*j+1];
        z=gdata->vertxyz[ND*j+2];

        PyTuple_SetItem(pTuple,0,PyFloat_FromDouble(x)); /* The tuple steals the reference */
        PyTuple_SetItem(pTuple,1,PyFloat_FromDouble(y));
        PyTuple_SetItem(pTuple,2,PyFloat_FromDouble(z));

        /* Call the Python function */
        pVal = PyObject_CallObject(pFunc,pTuple);  /* new reference */

        /* set magnetization */
        for (int i=0; i<ND; i++) {
          ta_M[ND*j+i]=PyFloat_AsDouble(PyTuple_GetItem(pVal,i));
        }
      }

      ierr = VecRestoreArray(gdata->M,&ta_M);CHKERRQ(ierr);

      /* decrement references */
      Py_DECREF(pFunc);
      Py_DECREF(pTuple);
      Py_DECREF(pVal);
    }
#else
    SETERRQ(PETSC_ERR_ARG_CORRUPT,"magpar not linked with Python interpreter. Cannot execute Python scripts!\n");
#endif
  }

  if (rank) {
    VecAssemblyBegin(gdata->M);
    VecAssemblyEnd(gdata->M);

    /* rescale magnetization (required by init_mag==5) */
    RenormVec(gdata->M,1.0,0.0,ND);

    if (init_mag<0) {
      ierr = VecScale(gdata->M,-1.0);CHKERRQ(ierr);
    }
    MagparFunctionLogReturn(0);
  }
  else {
    PetscRandom rctx;
    ierr = PetscRandomCreate(PETSC_COMM_SELF,&rctx);CHKERRQ(ierr);
    ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
    /* assign property ids to vertices */
    ierr = VertProp(gdata);CHKERRQ(ierr);

    for (int i=0; i<gdata->n_vert; i++) {
      if (pid>=0 && gdata->vertprop[i]!=pid) continue;
      switch(PetscAbsInt(init_mag)) {
        case 0:
          /* taken care of before switch */
        break;
        case 1: /* Mx=1 */
          ierr = VecSetValue(gdata->M,ND*i+0,1.0,INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+1,0.0,INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+2,0.0,INSERT_VALUES);CHKERRQ(ierr);
        break;
        case 2: /* My=1 */
          ierr = VecSetValue(gdata->M,ND*i+0,0.0,INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+1,1.0,INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+2,0.0,INSERT_VALUES);CHKERRQ(ierr);
        break;
        case 3: /* Mz=1 */
          ierr = VecSetValue(gdata->M,ND*i+0,0.0,INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+1,0.0,INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+2,1.0,INSERT_VALUES);CHKERRQ(ierr);
        break;
        case 4: /* Mx=My=Mz=0.57735027 */
          ierr = VecSetValue(gdata->M,ND*i+0,0.57735027,INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+1,0.57735027,INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+2,0.57735027,INSERT_VALUES);CHKERRQ(ierr);
        break;
        case 5: /* artificial flower state */
          ierr = VecSetValue(gdata->M,ND*i+0,gdata->vertxyz[ND*i+0]-init_magparm,INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+1,gdata->vertxyz[ND*i+1]-init_magparm,INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+2,1.0,INSERT_VALUES);CHKERRQ(ierr);
        break;
        case 6: /* set magnetization in x-z plane to theta=init_magpar */
          ierr = VecSetValue(gdata->M,ND*i+0,sin(init_magparm),INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+1,0.0,INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+2,cos(init_magparm),INSERT_VALUES);CHKERRQ(ierr);
        break;
        case 7: /* vortex centered at (0,0) with radius R=init_magparm:
                 * if R>0, then vortex is clockwise in xy plane.
                 * if R<0, then vortex is counter-clockwise in xy plane.
                 * center of vortex points in -z direction.
                 */
          {
            PetscReal x, y, r, R, A;
            PetscReal modmx,modmy,modmz;

            x = gdata->vertxyz[ND*i+0];
            y = gdata->vertxyz[ND*i+1];
            r = sqrt(x*x + y*y);
            R = init_magparm;

            A = 2.0*R*r/(R*R+r*r);
            if (r <= fabs(R)) {
              if (r > 1e-10) {
                modmx = A*y/r;
                modmy = -A*x/r;
                modmz = -sqrt(1.0-A*A);
              } else {
                modmx=0.0;
                modmy=0.0;
                modmz=-1.0;
              }
            } else {
              modmx=y/r;
              modmy=-x/r;
              modmz=0.0;
            }
            ierr = VecSetValue(gdata->M,ND*i+0,modmx,INSERT_VALUES);CHKERRQ(ierr);
            ierr = VecSetValue(gdata->M,ND*i+1,modmy,INSERT_VALUES);CHKERRQ(ierr);
            ierr = VecSetValue(gdata->M,ND*i+2,modmz,INSERT_VALUES);CHKERRQ(ierr);
          }
        break;
        case 8: /* random magnetization - inaccurate: not uniform on unit sphere! */
          PetscReal rval[ND],rnorm;

          PetscRandomGetValue(rctx,&rval[0]);
          PetscRandomGetValue(rctx,&rval[1]);
          PetscRandomGetValue(rctx,&rval[2]);
          rnorm=my_dnrm2(ND,rval,1);

          for (int k=0; k<ND; k++) {
            ierr = VecSetValue(
              gdata->M,
              ND*i+k,
              rval[k]/rnorm*2.0-1.0,
              INSERT_VALUES
            );CHKERRQ(ierr);
          }
        break;
        case 9: /* Bloch wall at x=+9.5 */
          if (gdata->vertxyz[ND*i+0] < init_magparm*0.95) {
            ierr = VecSetValue(gdata->M,ND*i+0, 0.0,INSERT_VALUES);CHKERRQ(ierr);
            ierr = VecSetValue(gdata->M,ND*i+1, 0.0,INSERT_VALUES);CHKERRQ(ierr);
            ierr = VecSetValue(gdata->M,ND*i+2,+1.0,INSERT_VALUES);CHKERRQ(ierr);
          }
          else if (gdata->vertxyz[ND*i+0] > init_magparm*1.05) {
            ierr = VecSetValue(gdata->M,ND*i+0, 0.0,INSERT_VALUES);CHKERRQ(ierr);
            ierr = VecSetValue(gdata->M,ND*i+1, 0.0,INSERT_VALUES);CHKERRQ(ierr);
            ierr = VecSetValue(gdata->M,ND*i+2,-1.0,INSERT_VALUES);CHKERRQ(ierr);
          }
          else {
            ierr = VecSetValue(gdata->M,ND*i+0, 0.1,INSERT_VALUES);CHKERRQ(ierr);
            ierr = VecSetValue(gdata->M,ND*i+1,-0.9,INSERT_VALUES);CHKERRQ(ierr);
            ierr = VecSetValue(gdata->M,ND*i+2, 0.1,INSERT_VALUES);CHKERRQ(ierr);
          }
        break;
        case 10: /* magnetization parallel to anisotropy axes */
          for (int k=0; k<ND; k++) {
            ierr = VecSetValue(
              gdata->M,
              ND*i+k,
              gdata->propdat[NP*gdata->vertprop[i]+6+k],
              INSERT_VALUES
            );CHKERRQ(ierr);
          }
        break;
        case 11: /* set magnetization in x-y plane to alpha=init_magparm */
          ierr = VecSetValue(gdata->M,ND*i+0,cos(init_magparm*M_PI/180.0),INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+1,sin(init_magparm*M_PI/180.0),INSERT_VALUES);CHKERRQ(ierr);
          ierr = VecSetValue(gdata->M,ND*i+2,0.0,INSERT_VALUES);CHKERRQ(ierr);
        break;
        case 12: /* Head-to-head transverse wall at x=init_magparm */
          if (gdata->vertxyz[ND*i+0] < init_magparm*0.95) {
            ierr=VecSetValue(gdata->M,ND*i+0, 1.0,INSERT_VALUES);CHKERRQ(ierr);
            ierr=VecSetValue(gdata->M,ND*i+1, 0.0,INSERT_VALUES);CHKERRQ(ierr);
            ierr=VecSetValue(gdata->M,ND*i+2, 0.0,INSERT_VALUES);CHKERRQ(ierr);
          }
          else if (gdata->vertxyz[ND*i+0] > init_magparm*1.05) {
            ierr=VecSetValue(gdata->M,ND*i+0,-1.0,INSERT_VALUES);CHKERRQ(ierr);
            ierr=VecSetValue(gdata->M,ND*i+1, 0.0,INSERT_VALUES);CHKERRQ(ierr);
            ierr=VecSetValue(gdata->M,ND*i+2, 0.0,INSERT_VALUES);CHKERRQ(ierr);
          }
          else {
            ierr=VecSetValue(gdata->M,ND*i+0, 0.0,INSERT_VALUES);CHKERRQ(ierr);
            ierr=VecSetValue(gdata->M,ND*i+1, 1.0,INSERT_VALUES);CHKERRQ(ierr);
            ierr=VecSetValue(gdata->M,ND*i+2, 0.0,INSERT_VALUES);CHKERRQ(ierr);
          }
        break;
        case 13: /* Head-to-head vortex wall at (0,0) with radius R=init_magparm:
                  * if R>0, then vortex is clockwise in xy plane.
                  * if R<0, then vortex is counter-clockwise in xy plane.
                  * center of vortex points in -z direction.
                  */
          {
            PetscReal x, y, r, R, A;
            PetscReal modmx,modmy,modmz;

            x = gdata->vertxyz[ND*i+0];
            y = gdata->vertxyz[ND*i+1];
            r = sqrt(x*x + y*y);
            R = init_magparm;

            A = 2.0*R*r/(R*R+r*r);
            if (r <= fabs(R)) {
              if (r > 1e-10) {
                modmx = A*y/r;
                modmy = -A*x/r;
                modmz = -sqrt(1.0-A*A);
              } else {
                modmx = 0.0;
                modmy = 0.0;
                modmz = -1.0;
              }
            } else {
              modmx = (x<0.0 ? +1.0 : -1.0);
              modmy = 0.0;
              modmz = 0.0;
            }
            ierr = VecSetValue(gdata->M,ND*i+0,modmx,INSERT_VALUES);CHKERRQ(ierr);
            ierr = VecSetValue(gdata->M,ND*i+1,modmy,INSERT_VALUES);CHKERRQ(ierr);
            ierr = VecSetValue(gdata->M,ND*i+2,modmz,INSERT_VALUES);CHKERRQ(ierr);
          }
        break;
        case 14:
          /* taken care of before switch */
        break;
        default:
          SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"init_mag value %i is illegal!\n",PetscAbsInt(init_mag));
        break;
      }
    }

    ierr = PetscRandomDestroy(rctx);CHKERRQ(ierr);
    ierr = PetscFree(gdata->vertprop);CHKERRQ(ierr);
    gdata->vertprop=PETSC_NULL;
  }

  ierr = VecAssemblyBegin(gdata->M);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(gdata->M);CHKERRQ(ierr);

  /* rescale magnetization (required by init_mag==5) */
  RenormVec(gdata->M,1.0,0.0,ND);

  if (init_mag<0) {
    ierr = VecScale(gdata->M,-1.0);CHKERRQ(ierr);
  }

  MagparFunctionLogReturn(0);
}