MagparWiki: attachment:hext_ho.c of BugTracker/ExitBeforeEndWithHstepfile

   1 /*
   2     This file is part of magpar.
   3 
   4     Copyright (C) 2002-2009 Werner Scholz
   5 
   6     www:   http://magnet.atp.tuwien.ac.at/scholz/magpar/
   7     email: magpar(at)magnet.atp.tuwien.ac.at
   8 
   9     magpar is free software; you can redistribute it and/or modify
  10     it under the terms of the GNU General Public License as published by
  11     the Free Software Foundation; either version 2 of the License, or
  12     (at your option) any later version.
  13 
  14     magpar is distributed in the hope that it will be useful,
  15     but WITHOUT ANY WARRANTY; without even the implied warranty of
  16     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17     GNU General Public License for more details.
  18 
  19     You should have received a copy of the GNU General Public License
  20     along with magpar; if not, write to the Free Software
  21     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  22 */
  23 
  24 static char const Id[] = "$Id: hext_ho.c 2416 2009-01-08 05:47:52Z scholz $\n\n";
  25 static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";
  26 
  27 #include "griddata.h"
  28 #include "field.h"
  29 #include "util/util.h"
  30 
  31 #ifdef ADDONS
  32 #include "addons/addons.h"
  33 #endif
  34 
  35 static int       doinit=1;
  36 static int       fieldon=0;
  37 static Vec       VHthis;
  38 static Vec       VMpHext=PETSC_NULL;
  39 
  40 static PetscReal hextamp;     /**< strength of the external field */
  41 static PetscReal hextini;     /**< initial strength of the external field */
  42 static PetscReal hexttheta;   /**< azimutal angle */
  43 static PetscReal hextphi;     /**< polar angle */
  44 static PetscReal hextxyz[ND]; /**< direction of the external field in cartesian coordinates */
  45 static PetscReal hextsweep;   /**< sweeping rate of the external field */
  46 static PetscReal hextstep;    /**< step size of the external field */
  47 static PetscReal hextfinal;
  48 static PetscReal tstart;
  49 
  50 static PetscReal *tfdata=PETSC_NULL;
  51 static PetscReal *hstep_fdata=PETSC_NULL;
  52 static PetscReal hstep;
  53 
  54 
  55 int MpHext(Vec M,PetscReal *vres)
  56 {
  57   MagparFunctionInfoBegin;
  58 
  59   assert(VMpHext!=PETSC_NULL);
  60 
  61   /* calculate M//Hext */
  62   ierr = VecDot(M,VMpHext,vres);CHKERRQ(ierr);
  63 
  64   MagparFunctionInfoReturn(0);
  65 }
  66 
  67 
  68 int Hext_ho_hstep(GridData *gdata)
  69 {
  70   MagparFunctionInfoBegin;
  71 
  72   PetscReal hextnewamp;
  73 
  74   if (!fieldon) {
  75     MagparFunctionInfoReturn(0);
  76   }
  77 
  78   if (hstep_fdata!=PETSC_NULL) {
  79     hstep+=1;
  80     /* current field amplitude = initial field * scaling factor */
  81     hextnewamp=hextini*get_scalf_hstep(hstep_fdata,hstep);
  82     hextstep=hextnewamp-hextamp;
  83   }
  84   else if ((hextstep < 0.0 && hextamp <= hextfinal) ||
  85            (hextstep > 0.0 && hextamp >= hextfinal)) {
  86     ierr = PetscPrintf(PETSC_COMM_WORLD,
  87       "hextamp %g kA/m exceeds hextfinal %g kA/m\n",
  88       hextamp*gdata->hscale/(MU0*1000.0),hextfinal*gdata->hscale/(MU0*1000.0)
  89     );
  90     gdata->mode=99;
  91   }
  92 
  93   if (hextstep==0.0) {
  94     MagparFunctionInfoReturn(0);
  95   }
  96 
  97   hextamp += hextstep;
  98   gdata->equil=0;
  99 
 100   MagparFunctionInfoReturn(0);
 101 }
 102 
 103 int Hext_ho_hsweep(GridData *gdata)
 104 {
 105   MagparFunctionInfoBegin;
 106 
 107   if (!fieldon || hextsweep==0) {
 108     MagparFunctionInfoReturn(0);
 109   }
 110   else if (hextsweep < 0.0 && hextamp <= hextfinal) {
 111     ierr = PetscPrintf(PETSC_COMM_WORLD,
 112       "hextamp %g kA/m < hextfinal %g kA/m\n",
 113       hextamp*gdata->hscale/(MU0*1000.0),hextfinal*gdata->hscale/(MU0*1000.0)
 114     );
 115     gdata->mode=99;
 116   }
 117   else if (hextsweep > 0.0 && hextamp >= hextfinal) {
 118     PetscPrintf(PETSC_COMM_WORLD,
 119       "hextamp %g kA/m > hextfinal %g kA/m\n",
 120       hextamp*gdata->hscale/(MU0*1000.0),hextfinal*gdata->hscale/(MU0*1000.0)
 121     );
 122     gdata->mode=99;
 123   }
 124   else {
 125     hextamp = hextini + hextsweep*(gdata->time - tstart);
 126   }
 127 
 128   MagparFunctionInfoReturn(0);
 129 }
 130 
 131 int Hext_ho_ht_Init()
 132 {
 133   MagparFunctionLogBegin;
 134 
 135   int rank,size;
 136   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
 137   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
 138 
 139   char       str[256];
 140   PetscTruth flg;
 141   ierr = PetscOptionsGetString(PETSC_NULL,"-hext_ho_htfile",str,256,&flg);CHKERRQ(ierr);
 142 
 143   if (flg!=PETSC_TRUE) {
 144     PetscPrintf(PETSC_COMM_WORLD,
 145       "Option -hext_ho_htfile not found, continuing with scaling=1\n",
 146       str
 147     );
 148     MagparFunctionLogReturn(0);
 149   }
 150 
 151   FILE *fd=PETSC_NULL;
 152   if (!rank) {
 153     ierr = PetscFOpen(PETSC_COMM_WORLD,str,"r",&fd);CHKERRQ(ierr);
 154     if (!fd){
 155       SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",str);
 156     }
 157   }
 158 
 159   ierr = readht(fd,&tfdata,2);CHKERRQ(ierr);
 160 
 161   ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
 162 
 163   MagparFunctionLogReturn(0);
 164 }
 165 
 166 int Hext_ho_hstep_Init()
 167 {
 168   MagparFunctionLogBegin;
 169 
 170   int rank;
 171   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
 172 
 173   char str[256];
 174   PetscTruth flg;
 175   ierr = PetscOptionsGetString(PETSC_NULL,"-hext_ho_hstepfile",str,256,&flg);CHKERRQ(ierr);
 176 
 177   if (flg!=PETSC_TRUE) {
 178     PetscPrintf(PETSC_COMM_WORLD,
 179       "Option -hext_ho_hstepfile not found, continuing with scaling=1\n",
 180       str
 181     );
 182     MagparFunctionLogReturn(0);
 183   }
 184 
 185   FILE *fd=PETSC_NULL;
 186   if (!rank) {
 187     ierr = PetscFOpen(PETSC_COMM_WORLD,str,"r",&fd);CHKERRQ(ierr);
 188     if (!fd){
 189       SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",str);
 190     }
 191   }
 192 
 193   ierr = read_hstep_file(fd,&hstep_fdata,2);CHKERRQ(ierr);
 194   hstep = 0;
 195 
 196   ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
 197 
 198   MagparFunctionLogReturn(0);
 199 }
 200 
 201 
 202 int Hext_ho_Init(GridData *gdata)
 203 {
 204   PetscTruth   flg;
 205 
 206   MagparFunctionLogBegin;
 207 
 208   ierr = PetscOptionsGetReal(PETSC_NULL,"-hextini",&hextini,&flg);CHKERRQ(ierr);
 209   if (flg!=PETSC_TRUE) {
 210     hextini=0;
 211     PetscPrintf(PETSC_COMM_WORLD,
 212       "Option -hextini not found. Using default value %g\n",
 213       hextini
 214     );
 215   }
 216   ierr = PetscOptionsGetReal(PETSC_NULL,"-htheta",&hexttheta,&flg);CHKERRQ(ierr);
 217   if (flg!=PETSC_TRUE) {
 218     hexttheta=0;
 219     PetscPrintf(PETSC_COMM_WORLD,
 220       "Option -htheta not found. Using default value %g\n",
 221       hexttheta
 222     );
 223   }
 224   ierr = PetscOptionsGetReal(PETSC_NULL,"-hphi",&hextphi,&flg);CHKERRQ(ierr);
 225   if (flg!=PETSC_TRUE) {
 226     hextphi=0;
 227     PetscPrintf(PETSC_COMM_WORLD,
 228       "Option -hphi not found. Using default value %g\n",
 229       hextphi
 230     );
 231   }
 232   ierr = PetscOptionsGetReal(PETSC_NULL,"-hstep",&hextstep,&flg);CHKERRQ(ierr);
 233   if (flg!=PETSC_TRUE) {
 234     hextstep=0;
 235     PetscPrintf(PETSC_COMM_WORLD,
 236       "Option -hstep not found. Using default value %g\n",
 237       hextstep
 238     );
 239   }
 240   ierr = PetscOptionsGetReal(PETSC_NULL,"-hsweep",&hextsweep,&flg);CHKERRQ(ierr);
 241   if (flg!=PETSC_TRUE) {
 242     hextsweep=0;
 243     PetscPrintf(PETSC_COMM_WORLD,
 244       "Option -hsweep not found. Using default value %g\n",
 245       hextsweep
 246     );
 247   }
 248 
 249   if (PetscAbsReal(hextstep)>D_EPS && PetscAbsReal(hextsweep)>D_EPS) {
 250     SETERRQ(PETSC_ERR_ARG_INCOMP,"hstep and hsweep != 0.0\n");
 251   }
 252 
 253   if (gdata->mode==1 && PetscAbsReal(hextsweep)>D_EPS) {
 254     SETERRQ(PETSC_ERR_ARG_INCOMP,"mode == 1 and hsweep != 0.0\n");
 255   }
 256 
 257   hextxyz[0]=sin(hexttheta)*cos(hextphi);
 258   hextxyz[1]=sin(hexttheta)*sin(hextphi);
 259   hextxyz[2]=cos(hexttheta);
 260 
 261   /* create VMpHext */
 262   ierr = VecDuplicate(gdata->M,&VMpHext);CHKERRQ(ierr);
 263 
 264   PetscReal *ta_elevol;
 265   ierr = VecGetArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);
 266   for (int i=0; i<gdata->ln_ele; i++) {
 267     for (int j=0; j<NV; j++) {
 268       for (int k=0; k<ND; k++) {
 269         PetscReal matele;
 270         matele=gdata->propdat[NP*gdata->eleprop[i]+4]*ta_elevol[i]/NV*hextxyz[k];
 271         ierr = VecSetValue(
 272           VMpHext,
 273           ND*gdata->elevert[NV*i+j]+k,
 274           matele,
 275           ADD_VALUES
 276         );CHKERRQ(ierr);
 277       }
 278     }
 279   }
 280   ierr = VecRestoreArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);
 281 
 282   ierr = VecAssemblyBegin(VMpHext);CHKERRQ(ierr);
 283   ierr = VecAssemblyEnd(VMpHext);CHKERRQ(ierr);
 284 
 285   PetscReal Javg, matele;
 286   ierr = VecSum(VMpHext,&matele);CHKERRQ(ierr);
 287 
 288   ierr = VecSum(gdata->VMs3,&Javg);CHKERRQ(ierr);
 289   Javg /= ND*gdata->totvol;
 290   PetscPrintf(PETSC_COMM_WORLD,"average magnetization: %g T\n",Javg*gdata->hscale);
 291 
 292   PetscPrintf(PETSC_COMM_WORLD,
 293     "Javg: %g  vol: %g  matele: %g  matele/vol: %g\n",
 294     Javg,
 295     gdata->totvol,
 296     matele,
 297     matele/gdata->totvol
 298   );
 299 
 300   Javg = 1.0/(Javg*gdata->totvol);
 301   ierr = VecScale(VMpHext,Javg);CHKERRQ(ierr);
 302 
 303   if (hextini==0.0 && hextstep==0.0 && hextsweep==0.0) {
 304     fieldon=0;
 305     PetscPrintf(PETSC_COMM_WORLD,
 306       "Options -hextini, -hextstep, and -hextsweep == 0 (field off).\n"
 307     );
 308     PetscPrintf(PETSC_COMM_WORLD,"Hext_ho off\n");
 309     MagparFunctionLogReturn(0);
 310   }
 311 
 312   fieldon=1;
 313   PetscPrintf(PETSC_COMM_WORLD,"Hext_ho on\n");
 314 
 315   ierr = PetscOptionsGetReal(PETSC_NULL,"-hfinal",&hextfinal,&flg);CHKERRQ(ierr);
 316   if (flg!=PETSC_TRUE) {
 317     hextfinal=0;
 318     PetscPrintf(PETSC_COMM_WORLD,
 319       "Option -hfinal not found. Using default value %i\n",
 320       hextfinal
 321     );
 322   }
 323 
 324   ierr = VecDuplicate(gdata->M,&VHthis);CHKERRQ(ierr);
 325 
 326 
 327   /* convert from kA/m to SI units A/m */
 328   hextini   *= 1000.0;
 329   hextstep  *= 1000.0;
 330   hextfinal *= 1000.0;
 331   hextsweep *= 1000.0;
 332 
 333   PetscPrintf(MPI_COMM_WORLD,
 334     "Hext:\n"
 335     "theta: %g rad = %g deg    phi: %g rad = %g deg\n"
 336     "e_H: (%g, %g, %g)\n"
 337     "|Hini|: %g A/m = %g T\n"
 338     "Hini: (%g, %g, %g) T\n"
 339     "Hstep: %g A/m = %g T\n"
 340     "Hsweep: %g A/(m*ns) = %g T/ns\n"
 341     "hextfinal: %g A/m = %g T\n",
 342     hexttheta,hexttheta*180.0/PETSC_PI,
 343     hextphi,  hextphi*180.0/PETSC_PI,
 344     hextxyz[0],hextxyz[1],hextxyz[2],
 345     hextini,  hextini*MU0,
 346     hextxyz[0]*hextini*MU0,hextxyz[1]*hextini*MU0,hextxyz[2]*hextini*MU0,
 347     hextstep, hextstep*MU0,
 348     hextsweep,hextsweep*MU0,
 349     hextfinal,hextfinal*MU0
 350   );
 351 
 352   /* rescale external field */
 353   hextini   /= gdata->hscale/MU0;
 354   hextamp=hextini;
 355   hextstep  /= gdata->hscale/MU0;
 356   hextfinal /= gdata->hscale/MU0;
 357   hextsweep /= gdata->hscale/(1e9*MU0*gdata->tscale);
 358 
 359   PetscPrintf(MPI_COMM_WORLD,
 360     "Hext_scaled:\n"
 361     "Hini: %g\n"
 362     "Hstep: %g\n"
 363     "Hsweep: %g\n"
 364     "Hfinal: %g\n",
 365     hextini,
 366     hextstep,
 367     hextsweep,
 368     hextfinal
 369   );
 370 
 371   /* larmor frequency in external field */
 372   if (hextini > D_EPS) {
 373     PetscReal  t_larmorf2;
 374 
 375     t_larmorf2=GAMMA/MU0*hextini;
 376     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)));
 377   }
 378 
 379   tstart=gdata->time;
 380 
 381   ierr = Hext_ho_ht_Init();CHKERRQ(ierr);
 382   ierr = Hext_ho_hstep_Init();CHKERRQ(ierr);
 383 
 384   MagparFunctionLogReturn(0);
 385 }
 386 
 387 
 388 int Hext_ho(GridData *gdata,Vec VHextsum,PetscReal *hext)
 389 {
 390   MagparFunctionInfoBegin;
 391 
 392   /* initialize if necessary */
 393   if (doinit) {
 394     ierr = Hext_ho_Init(gdata);CHKERRQ(ierr);
 395     doinit=0;
 396   }
 397   if (!fieldon) {
 398     MagparFunctionInfoReturn(0);
 399   }
 400 
 401   /* get time in nanoseconds */
 402   PetscReal nstime;
 403   nstime=gdata->time*gdata->tscale*1e9;
 404 
 405   if (tfdata!=PETSC_NULL) {
 406     /* current field amplitude =  initial field * scaling factor */
 407     hextamp=hextini*getscalf(tfdata,nstime);
 408   }
 409 
 410   if (PetscAbsReal(hextsweep)>D_EPS) {
 411     ierr = Hext_ho_hsweep(gdata);CHKERRQ(ierr);
 412   }
 413 
 414   static PetscReal hextampbak=PETSC_MAX;
 415   static PetscReal hextxyzbak[ND]={PETSC_MAX,PETSC_MAX,PETSC_MAX};
 416   if (hextxyz[0]!=hextxyzbak[0] ||
 417       hextxyz[1]!=hextxyzbak[1] ||
 418       hextxyz[2]!=hextxyzbak[2] ||
 419       hextamp!=hextampbak) {
 420 
 421     ierr = VecSetVec(VHthis,hextxyz,ND);CHKERRQ(ierr);
 422     hextxyzbak[0]=hextxyz[0];
 423     hextxyzbak[1]=hextxyz[1];
 424     hextxyzbak[2]=hextxyz[2];
 425 
 426     ierr = VecScale(VHthis,hextamp);CHKERRQ(ierr);
 427     hextampbak=hextamp;
 428   }
 429 
 430   ierr = VecAXPY(VHextsum,1.0,VHthis);CHKERRQ(ierr);
 431 
 432   *hext=hextamp;
 433 
 434   MagparFunctionProfReturn(0);
 435 }