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

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

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

static int       doinit=1;
static int       fieldon=0;
static Vec       VHthis;
static Vec       VMpHext=PETSC_NULL;

static PetscReal hextamp;     
static PetscReal hextini;     
static PetscReal hexttheta;   
static PetscReal hextphi;     
static PetscReal hextxyz[ND]; 
static PetscReal hextsweep;   
static PetscReal hextstep;    
static PetscReal hextfinal;
static PetscReal tstart;

static PetscReal *tfdata=PETSC_NULL;
static PetscReal *hstep_fdata=PETSC_NULL;
static PetscReal hstep;


int MpHext(Vec M,PetscReal *vres)
{
  MagparFunctionInfoBegin;

  assert(VMpHext!=PETSC_NULL);

  /* calculate M//Hext */
  ierr = VecDot(M,VMpHext,vres);CHKERRQ(ierr);

  MagparFunctionInfoReturn(0);
}


int Hext_ho_hstep(GridData *gdata)
{
  MagparFunctionInfoBegin;

  PetscReal hextnewamp;

  if (!fieldon) {
    MagparFunctionInfoReturn(0);
  }

  if (hstep_fdata!=PETSC_NULL) {
    hstep+=1;
    /* current field amplitude = initial field * scaling factor */
    hextnewamp=hextini*get_scalf_hstep(hstep_fdata,hstep);
    hextstep=hextnewamp-hextamp;
  }
  else if ((hextstep < 0.0 && hextamp <= hextfinal) ||
           (hextstep > 0.0 && hextamp >= hextfinal)) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,
      "hextamp %g kA/m exceeds hextfinal %g kA/m\n",
      hextamp*gdata->hscale/(MU0*1000.0),hextfinal*gdata->hscale/(MU0*1000.0)
    );
    gdata->mode=99;
  }

  if (hextstep==0.0) {
    MagparFunctionInfoReturn(0);
  }

  hextamp += hextstep;
  gdata->equil=0;

  MagparFunctionInfoReturn(0);
}

int Hext_ho_hsweep(GridData *gdata)
{
  MagparFunctionInfoBegin;

  if (!fieldon || hextsweep==0) {
    MagparFunctionInfoReturn(0);
  }
  else if (hextsweep < 0.0 && hextamp <= hextfinal) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,
      "hextamp %g kA/m < hextfinal %g kA/m\n",
      hextamp*gdata->hscale/(MU0*1000.0),hextfinal*gdata->hscale/(MU0*1000.0)
    );
    gdata->mode=99;
  }
  else if (hextsweep > 0.0 && hextamp >= hextfinal) {
    PetscPrintf(PETSC_COMM_WORLD,
      "hextamp %g kA/m > hextfinal %g kA/m\n",
      hextamp*gdata->hscale/(MU0*1000.0),hextfinal*gdata->hscale/(MU0*1000.0)
    );
    gdata->mode=99;
  }
  else {
    hextamp = hextini + hextsweep*(gdata->time - tstart);
  }

  MagparFunctionInfoReturn(0);
}

int Hext_ho_ht_Init()
{
  MagparFunctionLogBegin;

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

  char       str[256];
  PetscTruth flg;
  ierr = PetscOptionsGetString(PETSC_NULL,"-hext_ho_htfile",str,256,&flg);CHKERRQ(ierr);

  if (flg!=PETSC_TRUE) {
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -hext_ho_htfile not found, continuing with scaling=1\n",
      str
    );
    MagparFunctionLogReturn(0);
  }

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

  ierr = readht(fd,&tfdata,2);CHKERRQ(ierr);

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

  MagparFunctionLogReturn(0);
}

int Hext_ho_hstep_Init()
{
  MagparFunctionLogBegin;

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

  char str[256];
  PetscTruth flg;
  ierr = PetscOptionsGetString(PETSC_NULL,"-hext_ho_hstepfile",str,256,&flg);CHKERRQ(ierr);

  if (flg!=PETSC_TRUE) {
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -hext_ho_hstepfile not found, continuing with scaling=1\n",
      str
    );
    MagparFunctionLogReturn(0);
  }

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

  ierr = read_hstep_file(fd,&hstep_fdata,2);CHKERRQ(ierr);
  hstep = 0;

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

  MagparFunctionLogReturn(0);
}


int Hext_ho_Init(GridData *gdata)
{
  PetscTruth   flg;

  MagparFunctionLogBegin;

  ierr = PetscOptionsGetReal(PETSC_NULL,"-hextini",&hextini,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    hextini=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -hextini not found. Using default value %g\n",
      hextini
    );
  }
  ierr = PetscOptionsGetReal(PETSC_NULL,"-htheta",&hexttheta,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    hexttheta=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -htheta not found. Using default value %g\n",
      hexttheta
    );
  }
  ierr = PetscOptionsGetReal(PETSC_NULL,"-hphi",&hextphi,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    hextphi=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -hphi not found. Using default value %g\n",
      hextphi
    );
  }
  ierr = PetscOptionsGetReal(PETSC_NULL,"-hstep",&hextstep,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    hextstep=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -hstep not found. Using default value %g\n",
      hextstep
    );
  }
  ierr = PetscOptionsGetReal(PETSC_NULL,"-hsweep",&hextsweep,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    hextsweep=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -hsweep not found. Using default value %g\n",
      hextsweep
    );
  }

  if (PetscAbsReal(hextstep)>D_EPS && PetscAbsReal(hextsweep)>D_EPS) {
    SETERRQ(PETSC_ERR_ARG_INCOMP,"hstep and hsweep != 0.0\n");
  }

  if (gdata->mode==1 && PetscAbsReal(hextsweep)>D_EPS) {
    SETERRQ(PETSC_ERR_ARG_INCOMP,"mode == 1 and hsweep != 0.0\n");
  }

  hextxyz[0]=sin(hexttheta)*cos(hextphi);
  hextxyz[1]=sin(hexttheta)*sin(hextphi);
  hextxyz[2]=cos(hexttheta);

  /* create VMpHext */
  ierr = VecDuplicate(gdata->M,&VMpHext);CHKERRQ(ierr);

  PetscReal *ta_elevol;
  ierr = VecGetArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);
  for (int i=0; i<gdata->ln_ele; i++) {
    for (int j=0; j<NV; j++) {
      for (int k=0; k<ND; k++) {
        PetscReal matele;
        matele=gdata->propdat[NP*gdata->eleprop[i]+4]*ta_elevol[i]/NV*hextxyz[k];
        ierr = VecSetValue(
          VMpHext,
          ND*gdata->elevert[NV*i+j]+k,
          matele,
          ADD_VALUES
        );CHKERRQ(ierr);
      }
    }
  }
  ierr = VecRestoreArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);

  ierr = VecAssemblyBegin(VMpHext);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(VMpHext);CHKERRQ(ierr);

  PetscReal Javg, matele;
  ierr = VecSum(VMpHext,&matele);CHKERRQ(ierr);

  ierr = VecSum(gdata->VMs3,&Javg);CHKERRQ(ierr);
  Javg /= ND*gdata->totvol;
  PetscPrintf(PETSC_COMM_WORLD,"average magnetization: %g T\n",Javg*gdata->hscale);

  PetscPrintf(PETSC_COMM_WORLD,
    "Javg: %g  vol: %g  matele: %g  matele/vol: %g\n",
    Javg,
    gdata->totvol,
    matele,
    matele/gdata->totvol
  );

  Javg = 1.0/(Javg*gdata->totvol);
  ierr = VecScale(VMpHext,Javg);CHKERRQ(ierr);

  if (hextini==0.0 && hextstep==0.0 && hextsweep==0.0) {
    fieldon=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Options -hextini, -hextstep, and -hextsweep == 0 (field off).\n"
    );
    PetscPrintf(PETSC_COMM_WORLD,"Hext_ho off\n");
    MagparFunctionLogReturn(0);
  }

  fieldon=1;
  PetscPrintf(PETSC_COMM_WORLD,"Hext_ho on\n");

  ierr = PetscOptionsGetReal(PETSC_NULL,"-hfinal",&hextfinal,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    hextfinal=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -hfinal not found. Using default value %i\n",
      hextfinal
    );
  }

  ierr = VecDuplicate(gdata->M,&VHthis);CHKERRQ(ierr);


  /* convert from kA/m to SI units A/m */
  hextini   *= 1000.0;
  hextstep  *= 1000.0;
  hextfinal *= 1000.0;
  hextsweep *= 1000.0;

  PetscPrintf(MPI_COMM_WORLD,
    "Hext:\n"
    "theta: %g rad = %g deg    phi: %g rad = %g deg\n"
    "e_H: (%g, %g, %g)\n"
    "|Hini|: %g A/m = %g T\n"
    "Hini: (%g, %g, %g) T\n"
    "Hstep: %g A/m = %g T\n"
    "Hsweep: %g A/(m*ns) = %g T/ns\n"
    "hextfinal: %g A/m = %g T\n",
    hexttheta,hexttheta*180.0/PETSC_PI,
    hextphi,  hextphi*180.0/PETSC_PI,
    hextxyz[0],hextxyz[1],hextxyz[2],
    hextini,  hextini*MU0,
    hextxyz[0]*hextini*MU0,hextxyz[1]*hextini*MU0,hextxyz[2]*hextini*MU0,
    hextstep, hextstep*MU0,
    hextsweep,hextsweep*MU0,
    hextfinal,hextfinal*MU0
  );

  /* rescale external field */
  hextini   /= gdata->hscale/MU0;
  hextamp=hextini;
  hextstep  /= gdata->hscale/MU0;
  hextfinal /= gdata->hscale/MU0;
  hextsweep /= gdata->hscale/(1e9*MU0*gdata->tscale);

  PetscPrintf(MPI_COMM_WORLD,
    "Hext_scaled:\n"
    "Hini: %g\n"
    "Hstep: %g\n"
    "Hsweep: %g\n"
    "Hfinal: %g\n",
    hextini,
    hextstep,
    hextsweep,
    hextfinal
  );

  /* larmor frequency in external field */
  if (hextini > D_EPS) {
    PetscReal  t_larmorf2;

    t_larmorf2=GAMMA/MU0*hextini;
    PetscPrintf(PETSC_COMM_WORLD,"f_Larmor_ext: %g GHz -> T=1/f=%g ns\n",t_larmorf2/(2.0*PETSC_PI*1e9),1.0/(t_larmorf2/(2.0*PETSC_PI*1e9)));
  }

  tstart=gdata->time;

  ierr = Hext_ho_ht_Init();CHKERRQ(ierr);
  ierr = Hext_ho_hstep_Init();CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}


int Hext_ho(GridData *gdata,Vec VHextsum,PetscReal *hext)
{
  MagparFunctionInfoBegin;

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

  /* get time in nanoseconds */
  PetscReal nstime;
  nstime=gdata->time*gdata->tscale*1e9;

  if (tfdata!=PETSC_NULL) {
    /* current field amplitude =  initial field * scaling factor */
    hextamp=hextini*getscalf(tfdata,nstime);
  }

  if (PetscAbsReal(hextsweep)>D_EPS) {
    ierr = Hext_ho_hsweep(gdata);CHKERRQ(ierr);
  }

  static PetscReal hextampbak=PETSC_MAX;
  static PetscReal hextxyzbak[ND]={PETSC_MAX,PETSC_MAX,PETSC_MAX};
  if (hextxyz[0]!=hextxyzbak[0] ||
      hextxyz[1]!=hextxyzbak[1] ||
      hextxyz[2]!=hextxyzbak[2] ||
      hextamp!=hextampbak) {

    ierr = VecSetVec(VHthis,hextxyz,ND);CHKERRQ(ierr);
    hextxyzbak[0]=hextxyz[0];
    hextxyzbak[1]=hextxyz[1];
    hextxyzbak[2]=hextxyz[2];

    ierr = VecScale(VHthis,hextamp);CHKERRQ(ierr);
    hextampbak=hextamp;
  }

  ierr = VecAXPY(VHextsum,1.0,VHthis);CHKERRQ(ierr);

  *hext=hextamp;

  MagparFunctionProfReturn(0);
}