00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 static char const Id[] = "$Id: hext_ho.c 2681 2009-07-31 04:30:53Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "field.h"
00028 #include "util/util.h"
00029
00030 #ifdef ADDONS
00031 #include "addons/addons.h"
00032 #endif
00033
00034 static int doinit=1;
00035 static int fieldon=0;
00036 static Vec VHthis;
00037 static Vec VMpHext=PETSC_NULL;
00038
00039 static PetscReal hextamp;
00040 static PetscReal hextini;
00041 static PetscReal hexttheta;
00042 static PetscReal hextphi;
00043 static PetscReal hextxyz[ND];
00044 static PetscReal hextsweep;
00045 static PetscReal hextstep;
00046 static PetscReal hextfinal;
00047 static PetscReal tstart;
00048
00049 static PetscReal *tfdata=PETSC_NULL;
00050 static PetscReal *hstep_fdata=PETSC_NULL;
00051 static PetscReal hstep;
00052
00053
00054 int MpHext(Vec M,PetscReal *vres)
00055 {
00056 MagparFunctionInfoBegin;
00057
00058 assert(VMpHext!=PETSC_NULL);
00059
00060
00061 ierr = VecDot(M,VMpHext,vres);CHKERRQ(ierr);
00062
00063 MagparFunctionInfoReturn(0);
00064 }
00065
00066
00067 int Hext_ho_hstep(GridData *gdata)
00068 {
00069 MagparFunctionInfoBegin;
00070
00071 PetscReal hextnewamp;
00072
00073 if (!fieldon) {
00074 MagparFunctionInfoReturn(0);
00075 }
00076
00077 if (hstep_fdata!=PETSC_NULL) {
00078 hstep+=1;
00079
00080 hextnewamp=hextini*get_scalf_hstep(hstep_fdata,hstep);
00081 hextstep=hextnewamp-hextamp;
00082 }
00083 else if ((hextstep < 0.0 && hextamp <= hextfinal) ||
00084 (hextstep > 0.0 && hextamp >= hextfinal)) {
00085 ierr = PetscPrintf(PETSC_COMM_WORLD,
00086 "hextamp %g kA/m exceeds hextfinal %g kA/m\n",
00087 hextamp*gdata->hscale/(MU0*1000.0),hextfinal*gdata->hscale/(MU0*1000.0)
00088 );
00089 gdata->mode=99;
00090 }
00091
00092 if (hextstep==0.0) {
00093 MagparFunctionInfoReturn(0);
00094 }
00095
00096 hextamp += hextstep;
00097 gdata->equil=0;
00098
00099 MagparFunctionInfoReturn(0);
00100 }
00101
00102 int Hext_ho_hsweep(GridData *gdata)
00103 {
00104 MagparFunctionInfoBegin;
00105
00106 if (!fieldon || hextsweep==0) {
00107 MagparFunctionInfoReturn(0);
00108 }
00109 else if (hextsweep < 0.0 && hextamp <= hextfinal) {
00110 ierr = PetscPrintf(PETSC_COMM_WORLD,
00111 "hextamp %g kA/m < hextfinal %g kA/m\n",
00112 hextamp*gdata->hscale/(MU0*1000.0),hextfinal*gdata->hscale/(MU0*1000.0)
00113 );
00114 gdata->mode=99;
00115 }
00116 else if (hextsweep > 0.0 && hextamp >= hextfinal) {
00117 PetscPrintf(PETSC_COMM_WORLD,
00118 "hextamp %g kA/m > hextfinal %g kA/m\n",
00119 hextamp*gdata->hscale/(MU0*1000.0),hextfinal*gdata->hscale/(MU0*1000.0)
00120 );
00121 gdata->mode=99;
00122 }
00123 else {
00124 hextamp = hextini + hextsweep*(gdata->time - tstart);
00125 }
00126
00127 MagparFunctionInfoReturn(0);
00128 }
00129
00130 int Hext_ho_ht_Init()
00131 {
00132 MagparFunctionLogBegin;
00133
00134 int rank,size;
00135 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00136 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00137
00138 char str[256];
00139 PetscTruth flg;
00140 ierr = PetscOptionsGetString(PETSC_NULL,"-hext_ho_htfile",str,256,&flg);CHKERRQ(ierr);
00141
00142 if (flg!=PETSC_TRUE) {
00143 PetscPrintf(PETSC_COMM_WORLD,
00144 "Option -hext_ho_htfile not found, continuing with scaling=1\n",
00145 str
00146 );
00147 MagparFunctionLogReturn(0);
00148 }
00149
00150 FILE *fd=PETSC_NULL;
00151 if (!rank) {
00152 ierr = PetscFOpen(PETSC_COMM_WORLD,str,"r",&fd);CHKERRQ(ierr);
00153 if (!fd){
00154 SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",str);
00155 }
00156 }
00157
00158 ierr = readht(fd,&tfdata,2);CHKERRQ(ierr);
00159
00160 ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
00161
00162 MagparFunctionLogReturn(0);
00163 }
00164
00165 int Hext_ho_hstep_Init()
00166 {
00167 MagparFunctionLogBegin;
00168
00169 int rank;
00170 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00171
00172 char str[256];
00173 PetscTruth flg;
00174 ierr = PetscOptionsGetString(PETSC_NULL,"-hext_ho_hstepfile",str,256,&flg);CHKERRQ(ierr);
00175
00176 if (flg!=PETSC_TRUE) {
00177 PetscPrintf(PETSC_COMM_WORLD,
00178 "Option -hext_ho_hstepfile not found, continuing with scaling=1\n",
00179 str
00180 );
00181 MagparFunctionLogReturn(0);
00182 }
00183
00184 FILE *fd=PETSC_NULL;
00185 if (!rank) {
00186 ierr = PetscFOpen(PETSC_COMM_WORLD,str,"r",&fd);CHKERRQ(ierr);
00187 if (!fd){
00188 SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",str);
00189 }
00190 }
00191
00192 ierr = read_hstep_file(fd,&hstep_fdata,2);CHKERRQ(ierr);
00193 hstep = 0;
00194
00195 ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
00196
00197 MagparFunctionLogReturn(0);
00198 }
00199
00200
00201 int Hext_ho_Init(GridData *gdata)
00202 {
00203 PetscTruth flg;
00204
00205 MagparFunctionLogBegin;
00206
00207 ierr = PetscOptionsGetReal(PETSC_NULL,"-hextini",&hextini,&flg);CHKERRQ(ierr);
00208 if (flg!=PETSC_TRUE) {
00209 hextini=0;
00210 PetscPrintf(PETSC_COMM_WORLD,
00211 "Option -hextini not found. Using default value %g\n",
00212 hextini
00213 );
00214 }
00215 ierr = PetscOptionsGetReal(PETSC_NULL,"-htheta",&hexttheta,&flg);CHKERRQ(ierr);
00216 if (flg!=PETSC_TRUE) {
00217 hexttheta=0;
00218 PetscPrintf(PETSC_COMM_WORLD,
00219 "Option -htheta not found. Using default value %g\n",
00220 hexttheta
00221 );
00222 }
00223 ierr = PetscOptionsGetReal(PETSC_NULL,"-hphi",&hextphi,&flg);CHKERRQ(ierr);
00224 if (flg!=PETSC_TRUE) {
00225 hextphi=0;
00226 PetscPrintf(PETSC_COMM_WORLD,
00227 "Option -hphi not found. Using default value %g\n",
00228 hextphi
00229 );
00230 }
00231 ierr = PetscOptionsGetReal(PETSC_NULL,"-hstep",&hextstep,&flg);CHKERRQ(ierr);
00232 if (flg!=PETSC_TRUE) {
00233 hextstep=0;
00234 PetscPrintf(PETSC_COMM_WORLD,
00235 "Option -hstep not found. Using default value %g\n",
00236 hextstep
00237 );
00238 }
00239 ierr = PetscOptionsGetReal(PETSC_NULL,"-hsweep",&hextsweep,&flg);CHKERRQ(ierr);
00240 if (flg!=PETSC_TRUE) {
00241 hextsweep=0;
00242 PetscPrintf(PETSC_COMM_WORLD,
00243 "Option -hsweep not found. Using default value %g\n",
00244 hextsweep
00245 );
00246 }
00247
00248 if (PetscAbsReal(hextstep)>D_EPS && PetscAbsReal(hextsweep)>D_EPS) {
00249 SETERRQ(PETSC_ERR_ARG_INCOMP,"hstep and hsweep != 0.0\n");
00250 }
00251
00252 if (gdata->mode==1 && PetscAbsReal(hextsweep)>D_EPS) {
00253 SETERRQ(PETSC_ERR_ARG_INCOMP,"mode == 1 and hsweep != 0.0\n");
00254 }
00255
00256 hextxyz[0]=sin(hexttheta)*cos(hextphi);
00257 hextxyz[1]=sin(hexttheta)*sin(hextphi);
00258 hextxyz[2]=cos(hexttheta);
00259
00260
00261 ierr = VecDuplicate(gdata->M,&VMpHext);CHKERRQ(ierr);
00262
00263 PetscReal *ta_elevol;
00264 ierr = VecGetArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);
00265 for (int i=0; i<gdata->ln_ele; i++) {
00266 for (int j=0; j<NV; j++) {
00267 for (int k=0; k<ND; k++) {
00268 PetscReal matele;
00269 matele=gdata->propdat[NP*gdata->eleprop[i]+4]*ta_elevol[i]/NV*hextxyz[k];
00270 ierr = VecSetValue(
00271 VMpHext,
00272 ND*gdata->elevert[NV*i+j]+k,
00273 matele,
00274 ADD_VALUES
00275 );CHKERRQ(ierr);
00276 }
00277 }
00278 }
00279 ierr = VecRestoreArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);
00280
00281 ierr = VecAssemblyBegin(VMpHext);CHKERRQ(ierr);
00282 ierr = VecAssemblyEnd(VMpHext);CHKERRQ(ierr);
00283
00284 PetscReal Javg, matele;
00285 ierr = VecSum(VMpHext,&matele);CHKERRQ(ierr);
00286
00287 ierr = VecSum(gdata->VMs3,&Javg);CHKERRQ(ierr);
00288 Javg /= ND*gdata->totvol;
00289 PetscPrintf(PETSC_COMM_WORLD,"average magnetization: %g T\n",Javg*gdata->hscale);
00290
00291 PetscPrintf(PETSC_COMM_WORLD,
00292 "Javg: %g vol: %g matele: %g matele/vol: %g\n",
00293 Javg,
00294 gdata->totvol,
00295 matele,
00296 matele/gdata->totvol
00297 );
00298
00299 Javg = 1.0/(Javg*gdata->totvol);
00300 ierr = VecScale(VMpHext,Javg);CHKERRQ(ierr);
00301
00302 if (hextini==0.0 && hextstep==0.0 && hextsweep==0.0) {
00303 fieldon=0;
00304 PetscPrintf(PETSC_COMM_WORLD,
00305 "Options -hextini, -hextstep, and -hextsweep == 0 (field off).\n"
00306 );
00307 PetscPrintf(PETSC_COMM_WORLD,"Hext_ho off\n");
00308 MagparFunctionLogReturn(0);
00309 }
00310
00311 fieldon=1;
00312 PetscPrintf(PETSC_COMM_WORLD,"Hext_ho on\n");
00313
00314 ierr = PetscOptionsGetReal(PETSC_NULL,"-hfinal",&hextfinal,&flg);CHKERRQ(ierr);
00315 if (flg!=PETSC_TRUE) {
00316 hextfinal=0;
00317 PetscPrintf(PETSC_COMM_WORLD,
00318 "Option -hfinal not found. Using default value %i\n",
00319 hextfinal
00320 );
00321 }
00322
00323 ierr = VecDuplicate(gdata->M,&VHthis);CHKERRQ(ierr);
00324
00325
00326
00327 hextini *= 1000.0;
00328 hextstep *= 1000.0;
00329 hextfinal *= 1000.0;
00330 hextsweep *= 1000.0;
00331
00332 PetscPrintf(MPI_COMM_WORLD,
00333 "Hext:\n"
00334 "theta: %g rad = %g deg phi: %g rad = %g deg\n"
00335 "e_H: (%g, %g, %g)\n"
00336 "|Hini|: %g A/m = %g T\n"
00337 "Hini: (%g, %g, %g) T\n"
00338 "Hstep: %g A/m = %g T\n"
00339 "Hsweep: %g A/(m*ns) = %g T/ns\n"
00340 "hextfinal: %g A/m = %g T\n",
00341 hexttheta,hexttheta*180.0/PETSC_PI,
00342 hextphi, hextphi*180.0/PETSC_PI,
00343 hextxyz[0],hextxyz[1],hextxyz[2],
00344 hextini, hextini*MU0,
00345 hextxyz[0]*hextini*MU0,hextxyz[1]*hextini*MU0,hextxyz[2]*hextini*MU0,
00346 hextstep, hextstep*MU0,
00347 hextsweep,hextsweep*MU0,
00348 hextfinal,hextfinal*MU0
00349 );
00350
00351
00352 hextini /= gdata->hscale/MU0;
00353 hextamp=hextini;
00354 hextstep /= gdata->hscale/MU0;
00355 hextfinal /= gdata->hscale/MU0;
00356 hextsweep /= gdata->hscale/(1e9*MU0*gdata->tscale);
00357
00358 PetscPrintf(MPI_COMM_WORLD,
00359 "Hext_scaled:\n"
00360 "Hini: %g\n"
00361 "Hstep: %g\n"
00362 "Hsweep: %g\n"
00363 "Hfinal: %g\n",
00364 hextini,
00365 hextstep,
00366 hextsweep,
00367 hextfinal
00368 );
00369
00370
00371 if (hextini > D_EPS) {
00372 PetscReal t_larmorf2;
00373
00374 t_larmorf2=GAMMA/MU0*hextini;
00375 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)));
00376 }
00377
00378 tstart=gdata->time;
00379
00380 ierr = Hext_ho_ht_Init();CHKERRQ(ierr);
00381 ierr = Hext_ho_hstep_Init();CHKERRQ(ierr);
00382
00383 MagparFunctionLogReturn(0);
00384 }
00385
00386
00387 int Hext_ho(GridData *gdata,Vec VHextsum,PetscReal *hext)
00388 {
00389 MagparFunctionInfoBegin;
00390
00391
00392 if (doinit) {
00393 ierr = Hext_ho_Init(gdata);CHKERRQ(ierr);
00394 doinit=0;
00395 }
00396 if (!fieldon) {
00397 MagparFunctionInfoReturn(0);
00398 }
00399
00400
00401 PetscReal nstime;
00402 nstime=gdata->time*gdata->tscale*1e9;
00403
00404 if (tfdata!=PETSC_NULL) {
00405
00406 hextamp=hextini*getscalf(tfdata,nstime);
00407 }
00408
00409 if (PetscAbsReal(hextsweep)>D_EPS) {
00410 ierr = Hext_ho_hsweep(gdata);CHKERRQ(ierr);
00411 }
00412
00413 static PetscReal hextampbak=PETSC_MAX;
00414 static PetscReal hextxyzbak[ND]={PETSC_MAX,PETSC_MAX,PETSC_MAX};
00415 if (hextxyz[0]!=hextxyzbak[0] ||
00416 hextxyz[1]!=hextxyzbak[1] ||
00417 hextxyz[2]!=hextxyzbak[2] ||
00418 hextamp!=hextampbak) {
00419
00420 ierr = VecSetVec(VHthis,hextxyz,ND);CHKERRQ(ierr);
00421 hextxyzbak[0]=hextxyz[0];
00422 hextxyzbak[1]=hextxyz[1];
00423 hextxyzbak[2]=hextxyz[2];
00424
00425 ierr = VecScale(VHthis,hextamp);CHKERRQ(ierr);
00426 hextampbak=hextamp;
00427 }
00428
00429 ierr = VecAXPY(VHextsum,1.0,VHthis);CHKERRQ(ierr);
00430
00431 *hext=hextamp;
00432
00433 MagparFunctionProfReturn(0);
00434 }
00435