readkrn.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: readkrn.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "magpario.h"

#include <string.h>
#include <ctype.h>

int read_one_line(FILE *fd, char *line, int bufsize, char **comment_out)
{
  char *comment, *temp;
  int len, blank;

  do {
    if (fgets(line, bufsize, fd) == NULL) {
      line[0]=0; /* valid but empty */
      return(-1);
    }

    /* check that it's not too long */
    len = strlen(line);
    if (len==bufsize-1 && !(line[bufsize-1]=='\n' || feof(fd))) {
      /* line can't fit in 256-byte buffer */
      return(-1);
    }

    /* chop end-of-line \r\n or \n */
    if (len >= 2 && line[len-2]=='\r' && line[len-1]=='\n') {
      line[len-2] = 0;
    } else if (len >= 1 && line[len-1]=='\n') {
      line[len-1] = 0;
    }

    /* split comments */
    comment=strchr(line, '#');
    if (comment != NULL) {
      *comment=0;
      comment++;
    }

    /* skip blank lines (all space) */
    for (temp=line; *temp && isspace(*temp); temp++) {}
    if (*temp) { blank=0; } else { blank = 1; }
  } while (blank);

  if (comment_out) *comment_out = comment;
  return(1);
}

int ReadKrn(GridData *gdata)
{
  PetscTruth flg;

  MagparFunctionLogBegin;

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

  /* TODO: maginit needs it as along as M is a shared vector during serial init. */
  ierr = MPI_Bcast(&gdata->n_prop,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  PetscSynchronizedPrintf(PETSC_COMM_WORLD,
    "<%i>n_prop: %i\n",
    rank,gdata->n_prop
  );
  PetscSynchronizedFlush(PETSC_COMM_WORLD);

  ierr = PetscOptionsGetReal(PETSC_NULL,"-size",&gdata->lenscale,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    gdata->lenscale=1e-9;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -size not found, using default value: %g\n",
      gdata->lenscale
    );
  }

  PetscReal dum;
  ierr = PetscOptionsGetReal(PETSC_NULL,"-alpha",&dum,&flg);CHKERRQ(ierr);
  if (flg==PETSC_TRUE) {
    PetscPrintf(MPI_COMM_WORLD,
      "Warning: option -alpha will be ignored, "
      "using parameter from *.krn file instead!\n"
    );
  }

  /* SI value = internal (dimensionless) value * scaling factor */

  /* always scale fields and Js with 1.0 T, because the Js scaling parameter
     in allopt.txt only confused users
  */
  gdata->hscale=1.0;
  gdata->tscale=MU0/(gdata->hscale*GAMMA);
  gdata->escale=gdata->hscale*gdata->hscale/MU0;

  /* ***********************************************************************
     the parameters, which have been set above need NOT be broadcast
     to all other processors any more!
  */

  if (rank) {
    MagparFunctionLogReturn(0);
  }

  PetscPrintf(PETSC_COMM_WORLD,"\nscaling data:\n");
  PetscPrintf(PETSC_COMM_WORLD,"size: %g m\n",gdata->lenscale);

  PetscPrintf(MPI_COMM_WORLD,
    "magnetization, field: hscale: %g T = %g kA/m (hardcoded)\n",
    gdata->hscale,
    gdata->hscale/(MU0*1000.0)
  );
  PetscPrintf(PETSC_COMM_WORLD,"energy: escale: %g J/m^3 (hardcoded)\n",gdata->escale);
  PetscPrintf(PETSC_COMM_WORLD,"time:   tscale: %g ns (hardcoded)\n",gdata->tscale*1e9);

  /* TODO: read filename directly from option */
  char str[256];
  ierr = PetscSNPrintf(str,255,"%s%s",gdata->simname,".krn");CHKERRQ(ierr);

  FILE *fd;
  ierr = PetscFOpen(PETSC_COMM_WORLD,str,"r",&fd);CHKERRQ(ierr);
  if (!rank && !fd) {
    SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",str);
  }
  PetscPrintf(PETSC_COMM_WORLD,"Opening file %s\n",str);

/* TODO: valgrind reports lost memory:
==11675== 35,519 (14,679 direct, 20,840 indirect) bytes in 19 blocks are definitely lost in loss record 16 of 20
==11675==    at 0x4023D6E: malloc (vg_replace_malloc.c:207)
==11675==    by 0x84B9592: PetscMallocAlign(unsigned, int, char const*, char const*, char const*, void**) (in /home/scholz/work/magpar_trees/magpar/src/magpar.exe)
==11675==    by 0x80997F7: ReadKrn(GridData*) (readkrn.c:149)
==11675==    by 0x80786A7: SerInit(GridData*) (serinit.c:56)
==11675==    by 0x8063012: main (main.c:97)
*/

  ierr = PetscMalloc(gdata->n_prop*NP*sizeof(PetscReal),&gdata->propdat);CHKERRQ(ierr);
  ierr = PetscMemzero(gdata->propdat,gdata->n_prop*NP*sizeof(PetscReal));CHKERRQ(ierr);

  PetscReal t_larmorf=0.0;

  PetscPrintf(PETSC_COMM_WORLD,"\nscaled (dimensionless) material parameters:\n");
  for (int i=0; i<gdata->n_prop; i++) {
    char      line[256], *comment, *msg;
    PetscReal theta, phi, k1, k2, js, aexch, alpha;
    int       j, numchars;

    /* DRL: read a WHOLE LINE to avoid silent wraparound errors */
    read_one_line(fd, line, sizeof(line), &comment);

    j=sscanf(line,
#ifdef PETSC_USE_SINGLE
      "%f %f %f %f %f %f %f %n",
#else
      "%lf %lf %lf %lf %lf %lf %lf %n",
#endif
      &theta,&phi,&k1,&k2,&js,&aexch,&alpha,&numchars
    );
    if (j!=7) {
      SETERRQ2(PETSC_ERR_FILE_UNEXPECTED,"Data in %s corrupt: %i!\n",str,j);
    }

    /* read remainder of this line (possibly columns for cubic axes) */
    msg = &line[numchars];

    assert(NP==18+1);

#ifndef ADDONS
    /* set saturation magnetization js=0 if js<0
       unless ADDONS are compiled in
    */
    if (js<0.0) {
      js=0.0;
      PetscPrintf(MPI_COMM_WORLD,"Found Js<0 for material %i, resetting to Js=%g\n",i+1,js);
    }
#endif

    gdata->propdat[NP*i+ 0]=theta;
    gdata->propdat[NP*i+ 1]=phi;
    gdata->propdat[NP*i+ 2]=k1/gdata->escale;
    gdata->propdat[NP*i+ 3]=k2/gdata->escale;
    gdata->propdat[NP*i+ 4]=js/gdata->hscale;
    gdata->propdat[NP*i+ 5]=aexch/(gdata->escale*gdata->lenscale*gdata->lenscale);
    gdata->propdat[NP*i+ 9]=alpha;
    if (js == 0.0) {
      alpha=999;
      PetscPrintf(MPI_COMM_WORLD,"Found Js==0 for material %i, locking magnetization with alpha=%g\n",i+1,alpha);
    }

    PetscPrintf(MPI_COMM_WORLD,"----- property: %i -----\n",i);
    if (comment)
      PetscPrintf(MPI_COMM_WORLD, "comment:%s\n", comment);

    /* see if psi is there - then we have cubic anisotropy! */
    PetscReal psi;
    psi=PETSC_MIN;
#ifdef PETSC_USE_SINGLE
    j=sscanf(msg,"%f %n",&psi,&numchars);
#else
    j=sscanf(msg,"%lf %n",&psi,&numchars);
#endif
    if (j==1) {
      msg = &msg[numchars];
      /* specify the cubic anistropy axis via Euler angles,
         using the "x-convention"
         (let's pretend to be physicists- i.e.
         Goldstein's "Classical Mechanics",
         2nd ed. p143-148, Addison-Wesley (1980))
         e.g. http://mathworld.wolfram.com/EulerAngles.html
         first rotate about z-axis by phi,
         then rotate about new x-axis by theta,
         then rotate about new z-axis by psi
      */
      gdata->propdat[NP*i+16]=1;  /* cubic! */
      gdata->propdat[NP*i+ 6]= cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi);  /* 'x-axis' for cubic, x component */
      gdata->propdat[NP*i+ 7]= cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi);   /* y component */
      gdata->propdat[NP*i+ 8]= sin(psi)*sin(theta); /* z component */
      gdata->propdat[NP*i+10]=-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi);  /* 'y-axis' for cubic */
      gdata->propdat[NP*i+11]=-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi);
      gdata->propdat[NP*i+12]= cos(psi)*sin(theta);
      gdata->propdat[NP*i+13]= sin(phi)*sin(theta);  /* 'z-axis' for cubic */
      gdata->propdat[NP*i+14]=-cos(phi)*sin(theta);
      gdata->propdat[NP*i+15]= cos(theta);

      PetscPrintf(MPI_COMM_WORLD,
        "Euler angles: \n"
        " theta:  %g rad    phi: %g rad    psi :%g rad\n"
        " theta:  %g deg    phi: %g deg    psi :%g deg\n"
        "Cubic anisotropy axes: \n"
        " a-axis: (%g, %g, %g)\n"
        " b-axis: (%g, %g, %g)\n"
        " c-axis: (%g, %g, %g)\n"
        "K1:     %g  (cubic)\n"
        "K2:     %g  (cubic)\n",
        gdata->propdat[NP*i+ 0],
        gdata->propdat[NP*i+ 1],
        psi,
        gdata->propdat[NP*i+ 0]*180.0/PETSC_PI,
        gdata->propdat[NP*i+ 1]*180.0/PETSC_PI,
        psi*180.0/PETSC_PI,
        gdata->propdat[NP*i+ 6],
        gdata->propdat[NP*i+ 7],
        gdata->propdat[NP*i+ 8],
        gdata->propdat[NP*i+10],
        gdata->propdat[NP*i+11],
        gdata->propdat[NP*i+12],
        gdata->propdat[NP*i+13],
        gdata->propdat[NP*i+14],
        gdata->propdat[NP*i+15],
        gdata->propdat[NP*i+ 2],
        gdata->propdat[NP*i+ 3]
      );
    }
    else {
      char str2[256];
      j=sscanf(msg,"%s %n",&str2[0],&numchars);
      msg = &msg[numchars];

      clearerr(fd);    /* clear EOF and error indicators for fd */
      gdata->propdat[NP*i+16]=0;       /* uniaxial */
      gdata->propdat[NP*i+ 6]=sin(theta)*cos(phi);
      gdata->propdat[NP*i+ 7]=sin(theta)*sin(phi);
      gdata->propdat[NP*i+ 8]=cos(theta);

      PetscPrintf(MPI_COMM_WORLD,
        "theta:  %g rad = %g deg    phi: %g rad = %g deg\n"
        "c-axis: (%g, %g, %g)\n"
        "K1:     %g  (uniaxial)\n"
        "K2:     %g  (uniaxial)\n",
        gdata->propdat[NP*i+ 0],
        gdata->propdat[NP*i+ 0]*180.0/PETSC_PI,
        gdata->propdat[NP*i+ 1],
        gdata->propdat[NP*i+ 1]*180.0/PETSC_PI,
        gdata->propdat[NP*i+ 6],
        gdata->propdat[NP*i+ 7],
        gdata->propdat[NP*i+ 8],
        gdata->propdat[NP*i+ 2],
        gdata->propdat[NP*i+ 3]
      );
    }

    PetscReal mu,a;
#ifdef PETSC_USE_SINGLE
    j=sscanf(msg,"%f %f",&mu,&a);
#else
    j=sscanf(msg,"%lf %lf",&mu,&a);
#endif
    if (j==2) {
      gdata->propdat[NP*i+17]=mu; /* permeability */
      gdata->propdat[NP*i+18]=a;  /* a parameter for knee */
      PetscPrintf(MPI_COMM_WORLD,
        "mu:     %g\n"
        "a:      %g\n",
        gdata->propdat[NP*i+17],
        gdata->propdat[NP*i+18]
      );
      if ((gdata->propdat[NP*i+ 4]>0.0  && gdata->propdat[NP*i+17]==1.0) ||
          (gdata->propdat[NP*i+ 4]==0.0 && gdata->propdat[NP*i+17]!=1.0) ||
          (gdata->propdat[NP*i+17]<0.0)
      ) {
        SETERRQ2(PETSC_ERR_ARG_INCOMP,
          "Warning: Incorrect material with Ms=%g and mu=%g\n",
          gdata->propdat[NP*i+ 4],
          gdata->propdat[NP*i+17]
        );
      }
    }

    PetscPrintf(MPI_COMM_WORLD,
      "Js:     %g\n"
      "A:      %g\n"
      "alpha:  %g\n"
      "H_ani   = 2*K1/Js: %g A/m = %g T\n"
      "l_exch1 = sqrt(2*mu0*A/Js^2):      %g nm\n",
      gdata->propdat[NP*i+ 4],
      gdata->propdat[NP*i+ 5],
      gdata->propdat[NP*i+ 9],
      2.0*k1/js,
      2.0*MU0*k1/js,
      sqrt(2.0*aexch*MU0/(js*js))*1e9
    );
    if (gdata->propdat[NP*i+ 9]<0) {
      PetscPrintf(MPI_COMM_WORLD,
        "Warning: alpha=%g < 0! Using modified (large) alpha at t<0!\n",
        gdata->propdat[NP*i+ 9]
      );
    }
    if (gdata->propdat[NP*i+ 4]==0.0) {
      PetscPrintf(MPI_COMM_WORLD,"Found material with Js==0!\n");
    }
    /* l_exch1 and l_exch2, both are called "exchange length"!
       cf. Hubert, Schaefer, Magnetic Domains, p. 153
    */

    if (PetscAbsReal(k1)>D_EPS) {
      PetscPrintf(MPI_COMM_WORLD,
        "l_exch2 = sqrt(abs(A/K1)):         %g nm\n",
        sqrt(PetscAbsReal(aexch/k1))*1e9
      );

      PetscPrintf(MPI_COMM_WORLD,
        "l_Bloch = Pi*sqrt(abs(A/K1)):      %g nm\n",
        PETSC_PI*sqrt(PetscAbsReal(aexch/k1))*1e9
      );
    }

    PetscPrintf(MPI_COMM_WORLD,
      "E_Bloch = 4.0*sqrt(abs(A*K1)):     %g J/m^2\n",
      4.0*sqrt(PetscAbsReal(aexch*k1))
    );

    /* larmor frequency in anisotropy field */
    PetscReal t_larmorf2;
    if (k1 > D_EPS) {
      t_larmorf2=GAMMA*2.0*k1/js;
      PetscPrintf(PETSC_COMM_WORLD,"f_Larmor_ani: %g GHz -> T=1/f=%g ns\n",t_larmorf2/(2.0*PETSC_PI*1e9),1.0/(t_larmorf2/(2.0*PETSC_PI*1e9)));
      t_larmorf=PetscMax(t_larmorf,t_larmorf2);
    }

    /* larmor frequency in demag field
       (assumes demag factor N=1, e.g. thin film)
    */
    t_larmorf2=GAMMA*js/MU0;
    PetscPrintf(PETSC_COMM_WORLD,"f_Larmor_dem: %g GHz -> T=1/f=%g ns\n",t_larmorf2/(2.0*PETSC_PI*1e9),1.0/(t_larmorf2/(2.0*PETSC_PI*1e9)));
    t_larmorf=PetscMax(t_larmorf,t_larmorf2);
  }

  ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}