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: readkrn.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "magpario.h"
00028
00029 #include <string.h>
00030 #include <ctype.h>
00031
00032 int read_one_line(FILE *fd, char *line, int bufsize, char **comment_out)
00033 {
00034 char *comment, *temp;
00035 int len, blank;
00036
00037 do {
00038 if (fgets(line, bufsize, fd) == NULL) {
00039 line[0]=0;
00040 return(-1);
00041 }
00042
00043
00044 len = strlen(line);
00045 if (len==bufsize-1 && !(line[bufsize-1]=='\n' || feof(fd))) {
00046
00047 return(-1);
00048 }
00049
00050
00051 if (len >= 2 && line[len-2]=='\r' && line[len-1]=='\n') {
00052 line[len-2] = 0;
00053 } else if (len >= 1 && line[len-1]=='\n') {
00054 line[len-1] = 0;
00055 }
00056
00057
00058 comment=strchr(line, '#');
00059 if (comment != NULL) {
00060 *comment=0;
00061 comment++;
00062 }
00063
00064
00065 for (temp=line; *temp && isspace(*temp); temp++) {}
00066 if (*temp) { blank=0; } else { blank = 1; }
00067 } while (blank);
00068
00069 if (comment_out) *comment_out = comment;
00070 return(1);
00071 }
00072
00073 int ReadKrn(GridData *gdata)
00074 {
00075 PetscTruth flg;
00076
00077 MagparFunctionLogBegin;
00078
00079 int rank,size;
00080 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00081 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00082
00083
00084 ierr = MPI_Bcast(&gdata->n_prop,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00085 PetscSynchronizedPrintf(PETSC_COMM_WORLD,
00086 "<%i>n_prop: %i\n",
00087 rank,gdata->n_prop
00088 );
00089 PetscSynchronizedFlush(PETSC_COMM_WORLD);
00090
00091 ierr = PetscOptionsGetReal(PETSC_NULL,"-size",&gdata->lenscale,&flg);CHKERRQ(ierr);
00092 if (flg!=PETSC_TRUE) {
00093 gdata->lenscale=1e-9;
00094 PetscPrintf(PETSC_COMM_WORLD,
00095 "Option -size not found, using default value: %g\n",
00096 gdata->lenscale
00097 );
00098 }
00099
00100 PetscReal dum;
00101 ierr = PetscOptionsGetReal(PETSC_NULL,"-alpha",&dum,&flg);CHKERRQ(ierr);
00102 if (flg==PETSC_TRUE) {
00103 PetscPrintf(MPI_COMM_WORLD,
00104 "Warning: option -alpha will be ignored, "
00105 "using parameter from *.krn file instead!\n"
00106 );
00107 }
00108
00109
00110
00111
00112
00113
00114 gdata->hscale=1.0;
00115 gdata->tscale=MU0/(gdata->hscale*GAMMA);
00116 gdata->escale=gdata->hscale*gdata->hscale/MU0;
00117
00118
00119
00120
00121
00122
00123 if (rank) {
00124 MagparFunctionLogReturn(0);
00125 }
00126
00127 PetscPrintf(PETSC_COMM_WORLD,"\nscaling data:\n");
00128 PetscPrintf(PETSC_COMM_WORLD,"size: %g m\n",gdata->lenscale);
00129
00130 PetscPrintf(MPI_COMM_WORLD,
00131 "magnetization, field: hscale: %g T = %g kA/m (hardcoded)\n",
00132 gdata->hscale,
00133 gdata->hscale/(MU0*1000.0)
00134 );
00135 PetscPrintf(PETSC_COMM_WORLD,"energy: escale: %g J/m^3 (hardcoded)\n",gdata->escale);
00136 PetscPrintf(PETSC_COMM_WORLD,"time: tscale: %g ns (hardcoded)\n",gdata->tscale*1e9);
00137
00138
00139 char str[256];
00140 ierr = PetscSNPrintf(str,255,"%s%s",gdata->simname,".krn");CHKERRQ(ierr);
00141
00142 FILE *fd;
00143 ierr = PetscFOpen(PETSC_COMM_WORLD,str,"r",&fd);CHKERRQ(ierr);
00144 if (!rank && !fd) {
00145 SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",str);
00146 }
00147 PetscPrintf(PETSC_COMM_WORLD,"Opening file %s\n",str);
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158 ierr = PetscMalloc(gdata->n_prop*NP*sizeof(PetscReal),&gdata->propdat);CHKERRQ(ierr);
00159 ierr = PetscMemzero(gdata->propdat,gdata->n_prop*NP*sizeof(PetscReal));CHKERRQ(ierr);
00160
00161 PetscReal t_larmorf=0.0;
00162
00163 PetscPrintf(PETSC_COMM_WORLD,"\nscaled (dimensionless) material parameters:\n");
00164 for (int i=0; i<gdata->n_prop; i++) {
00165 char line[256], *comment, *msg;
00166 PetscReal theta, phi, k1, k2, js, aexch, alpha;
00167 int j, numchars;
00168
00169
00170 read_one_line(fd, line, sizeof(line), &comment);
00171
00172 j=sscanf(line,
00173 #ifdef PETSC_USE_SINGLE
00174 "%f %f %f %f %f %f %f %n",
00175 #else
00176 "%lf %lf %lf %lf %lf %lf %lf %n",
00177 #endif
00178 &theta,&phi,&k1,&k2,&js,&aexch,&alpha,&numchars
00179 );
00180 if (j!=7) {
00181 SETERRQ2(PETSC_ERR_FILE_UNEXPECTED,"Data in %s corrupt: %i!\n",str,j);
00182 }
00183
00184
00185 msg = &line[numchars];
00186
00187 assert(NP==18+1);
00188
00189 #ifndef ADDONS
00190
00191
00192
00193 if (js<0.0) {
00194 js=0.0;
00195 PetscPrintf(MPI_COMM_WORLD,"Found Js<0 for material %i, resetting to Js=%g\n",i+1,js);
00196 }
00197 #endif
00198
00199 gdata->propdat[NP*i+ 0]=theta;
00200 gdata->propdat[NP*i+ 1]=phi;
00201 gdata->propdat[NP*i+ 2]=k1/gdata->escale;
00202 gdata->propdat[NP*i+ 3]=k2/gdata->escale;
00203 gdata->propdat[NP*i+ 4]=js/gdata->hscale;
00204 gdata->propdat[NP*i+ 5]=aexch/(gdata->escale*gdata->lenscale*gdata->lenscale);
00205 gdata->propdat[NP*i+ 9]=alpha;
00206 if (js == 0.0) {
00207 alpha=999;
00208 PetscPrintf(MPI_COMM_WORLD,"Found Js==0 for material %i, locking magnetization with alpha=%g\n",i+1,alpha);
00209 }
00210
00211 PetscPrintf(MPI_COMM_WORLD,"----- property: %i -----\n",i);
00212 if (comment)
00213 PetscPrintf(MPI_COMM_WORLD, "comment:%s\n", comment);
00214
00215
00216 PetscReal psi;
00217 psi=PETSC_MIN;
00218 #ifdef PETSC_USE_SINGLE
00219 j=sscanf(msg,"%f %n",&psi,&numchars);
00220 #else
00221 j=sscanf(msg,"%lf %n",&psi,&numchars);
00222 #endif
00223 if (j==1) {
00224 msg = &msg[numchars];
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235 gdata->propdat[NP*i+16]=1;
00236 gdata->propdat[NP*i+ 6]= cos(psi)*cos(phi)-cos(theta)*sin(phi)*sin(psi);
00237 gdata->propdat[NP*i+ 7]= cos(psi)*sin(phi)+cos(theta)*cos(phi)*sin(psi);
00238 gdata->propdat[NP*i+ 8]= sin(psi)*sin(theta);
00239 gdata->propdat[NP*i+10]=-sin(psi)*cos(phi)-cos(theta)*sin(phi)*cos(psi);
00240 gdata->propdat[NP*i+11]=-sin(psi)*sin(phi)+cos(theta)*cos(phi)*cos(psi);
00241 gdata->propdat[NP*i+12]= cos(psi)*sin(theta);
00242 gdata->propdat[NP*i+13]= sin(phi)*sin(theta);
00243 gdata->propdat[NP*i+14]=-cos(phi)*sin(theta);
00244 gdata->propdat[NP*i+15]= cos(theta);
00245
00246 PetscPrintf(MPI_COMM_WORLD,
00247 "Euler angles: \n"
00248 " theta: %g rad phi: %g rad psi :%g rad\n"
00249 " theta: %g deg phi: %g deg psi :%g deg\n"
00250 "Cubic anisotropy axes: \n"
00251 " a-axis: (%g, %g, %g)\n"
00252 " b-axis: (%g, %g, %g)\n"
00253 " c-axis: (%g, %g, %g)\n"
00254 "K1: %g (cubic)\n"
00255 "K2: %g (cubic)\n",
00256 gdata->propdat[NP*i+ 0],
00257 gdata->propdat[NP*i+ 1],
00258 psi,
00259 gdata->propdat[NP*i+ 0]*180.0/PETSC_PI,
00260 gdata->propdat[NP*i+ 1]*180.0/PETSC_PI,
00261 psi*180.0/PETSC_PI,
00262 gdata->propdat[NP*i+ 6],
00263 gdata->propdat[NP*i+ 7],
00264 gdata->propdat[NP*i+ 8],
00265 gdata->propdat[NP*i+10],
00266 gdata->propdat[NP*i+11],
00267 gdata->propdat[NP*i+12],
00268 gdata->propdat[NP*i+13],
00269 gdata->propdat[NP*i+14],
00270 gdata->propdat[NP*i+15],
00271 gdata->propdat[NP*i+ 2],
00272 gdata->propdat[NP*i+ 3]
00273 );
00274 }
00275 else {
00276 char str2[256];
00277 j=sscanf(msg,"%s %n",&str2[0],&numchars);
00278 msg = &msg[numchars];
00279
00280 clearerr(fd);
00281 gdata->propdat[NP*i+16]=0;
00282 gdata->propdat[NP*i+ 6]=sin(theta)*cos(phi);
00283 gdata->propdat[NP*i+ 7]=sin(theta)*sin(phi);
00284 gdata->propdat[NP*i+ 8]=cos(theta);
00285
00286 PetscPrintf(MPI_COMM_WORLD,
00287 "theta: %g rad = %g deg phi: %g rad = %g deg\n"
00288 "c-axis: (%g, %g, %g)\n"
00289 "K1: %g (uniaxial)\n"
00290 "K2: %g (uniaxial)\n",
00291 gdata->propdat[NP*i+ 0],
00292 gdata->propdat[NP*i+ 0]*180.0/PETSC_PI,
00293 gdata->propdat[NP*i+ 1],
00294 gdata->propdat[NP*i+ 1]*180.0/PETSC_PI,
00295 gdata->propdat[NP*i+ 6],
00296 gdata->propdat[NP*i+ 7],
00297 gdata->propdat[NP*i+ 8],
00298 gdata->propdat[NP*i+ 2],
00299 gdata->propdat[NP*i+ 3]
00300 );
00301 }
00302
00303 PetscReal mu,a;
00304 #ifdef PETSC_USE_SINGLE
00305 j=sscanf(msg,"%f %f",&mu,&a);
00306 #else
00307 j=sscanf(msg,"%lf %lf",&mu,&a);
00308 #endif
00309 if (j==2) {
00310 gdata->propdat[NP*i+17]=mu;
00311 gdata->propdat[NP*i+18]=a;
00312 PetscPrintf(MPI_COMM_WORLD,
00313 "mu: %g\n"
00314 "a: %g\n",
00315 gdata->propdat[NP*i+17],
00316 gdata->propdat[NP*i+18]
00317 );
00318 if ((gdata->propdat[NP*i+ 4]>0.0 && gdata->propdat[NP*i+17]==1.0) ||
00319 (gdata->propdat[NP*i+ 4]==0.0 && gdata->propdat[NP*i+17]!=1.0) ||
00320 (gdata->propdat[NP*i+17]<0.0)
00321 ) {
00322 SETERRQ2(PETSC_ERR_ARG_INCOMP,
00323 "Warning: Incorrect material with Ms=%g and mu=%g\n",
00324 gdata->propdat[NP*i+ 4],
00325 gdata->propdat[NP*i+17]
00326 );
00327 }
00328 }
00329
00330 PetscPrintf(MPI_COMM_WORLD,
00331 "Js: %g\n"
00332 "A: %g\n"
00333 "alpha: %g\n"
00334 "H_ani = 2*K1/Js: %g A/m = %g T\n"
00335 "l_exch1 = sqrt(2*mu0*A/Js^2): %g nm\n",
00336 gdata->propdat[NP*i+ 4],
00337 gdata->propdat[NP*i+ 5],
00338 gdata->propdat[NP*i+ 9],
00339 2.0*k1/js,
00340 2.0*MU0*k1/js,
00341 sqrt(2.0*aexch*MU0/(js*js))*1e9
00342 );
00343 if (gdata->propdat[NP*i+ 9]<0) {
00344 PetscPrintf(MPI_COMM_WORLD,
00345 "Warning: alpha=%g < 0! Using modified (large) alpha at t<0!\n",
00346 gdata->propdat[NP*i+ 9]
00347 );
00348 }
00349 if (gdata->propdat[NP*i+ 4]==0.0) {
00350 PetscPrintf(MPI_COMM_WORLD,"Found material with Js==0!\n");
00351 }
00352
00353
00354
00355
00356 if (PetscAbsReal(k1)>D_EPS) {
00357 PetscPrintf(MPI_COMM_WORLD,
00358 "l_exch2 = sqrt(abs(A/K1)): %g nm\n",
00359 sqrt(PetscAbsReal(aexch/k1))*1e9
00360 );
00361
00362 PetscPrintf(MPI_COMM_WORLD,
00363 "l_Bloch = Pi*sqrt(abs(A/K1)): %g nm\n",
00364 PETSC_PI*sqrt(PetscAbsReal(aexch/k1))*1e9
00365 );
00366 }
00367
00368 PetscPrintf(MPI_COMM_WORLD,
00369 "E_Bloch = 4.0*sqrt(abs(A*K1)): %g J/m^2\n",
00370 4.0*sqrt(PetscAbsReal(aexch*k1))
00371 );
00372
00373
00374 PetscReal t_larmorf2;
00375 if (k1 > D_EPS) {
00376 t_larmorf2=GAMMA*2.0*k1/js;
00377 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)));
00378 t_larmorf=PetscMax(t_larmorf,t_larmorf2);
00379 }
00380
00381
00382
00383
00384 t_larmorf2=GAMMA*js/MU0;
00385 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)));
00386 t_larmorf=PetscMax(t_larmorf,t_larmorf2);
00387 }
00388
00389 ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
00390
00391 MagparFunctionLogReturn(0);
00392 }
00393