00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 static char const Id[] = "$Id: helastic.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
00026 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00027
00028 #include "field.h"
00029
00030 static int doinit=1;
00031 static int fieldon=0;
00032 static PetscReal *rval;
00033
00034 #define NPE 6
00035 static PetscReal *propdat;
00036 static Vec VHelastic;
00037
00038 int ReadProp(int n_prop)
00039 {
00040 PetscTruth flg;
00041
00042 MagparFunctionLogBegin;
00043
00044 int rank,size;
00045 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00046 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00047
00048 char str[256];
00049 ierr = PetscOptionsGetString(PETSC_NULL,"-helastic_propfile",str,256,&flg);CHKERRQ(ierr);
00050 if (flg!=PETSC_TRUE) {
00051 SETERRQ1(PETSC_ERR_ARG_INCOMP,"Problem with option -helastic_propfile %s\n",str);
00052 }
00053
00054 ierr = PetscMalloc(n_prop*NPE*sizeof(PetscReal),&propdat);CHKERRQ(ierr);
00055
00056 if (!rank) {
00057 FILE *fd;
00058 ierr = PetscFOpen(PETSC_COMM_WORLD,str,"r",&fd);CHKERRQ(ierr);
00059 if (!fd) {
00060 SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",str);
00061 }
00062 PetscPrintf(PETSC_COMM_WORLD,"Opening file %s\n",str);
00063
00064 for (int i=0; i<n_prop; i++) {
00065 for (int j=0; j<NPE; j++) {
00066 int k;
00067 #ifdef PETSC_USE_SINGLE
00068 k=fscanf(fd,"%f",propdat+NPE*i+j);
00069 #else
00070 k=fscanf(fd,"%lf",propdat+NPE*i+j);
00071 #endif
00072 if (k!=1) {
00073 SETERRQ2(PETSC_ERR_FILE_UNEXPECTED,"Data in %s corrupt: %i!\n",str,j);
00074 }
00075 }
00076
00077 fgets(str,256,fd);
00078
00079 PetscPrintf(MPI_COMM_WORLD,"----- property: %i -----\n",i);
00080 PetscPrintf(MPI_COMM_WORLD,
00081 "texture: %g\n"
00082 "lambda100: %g\n"
00083 "lambda111: %g\n"
00084 "sigmaX: %g\n"
00085 "sigmaY: %g\n"
00086 "sigmaZ: %g\n",
00087 propdat[NPE*i+0],
00088 propdat[NPE*i+1],
00089 propdat[NPE*i+2],
00090 propdat[NPE*i+3],
00091 propdat[NPE*i+4],
00092 propdat[NPE*i+5]
00093 );
00094 assert(NPE==6);
00095 }
00096
00097 ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
00098 }
00099 ierr = MPI_Bcast(propdat,NPE*n_prop,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00100
00101 MagparFunctionLogReturn(0);
00102 }
00103
00104
00105 int Helastic_Init(GridData *gdata)
00106 {
00107 PetscTruth flg;
00108
00109 MagparFunctionLogBegin;
00110
00111 ierr = PetscOptionsHasName(PETSC_NULL,"-helastic_propfile",&flg);CHKERRQ(ierr);
00112 if (flg!=PETSC_TRUE) {
00113 PetscPrintf(PETSC_COMM_WORLD,
00114 "Option -helastic_propfile not found, using default value: %i (field off)\n",
00115 fieldon
00116 );
00117 }
00118 else {
00119 fieldon=1;
00120 }
00121
00122 if (!fieldon) {
00123 PetscPrintf(PETSC_COMM_WORLD,"Helastic off\n");
00124 MagparFunctionLogReturn(0);
00125 }
00126 PetscPrintf(PETSC_COMM_WORLD,"Helastic on\n");
00127
00128
00129 ierr = ReadProp(gdata->n_prop);CHKERRQ(ierr);
00130
00131 ierr = VecDuplicate(gdata->M,&VHelastic);CHKERRQ(ierr);
00132 ierr = VecZeroEntries(VHelastic);CHKERRQ(ierr);
00133
00134
00135
00136
00137 ierr = PetscMalloc(ND*gdata->ln_vert*sizeof(PetscReal),&rval);CHKERRQ(ierr);
00138
00139 PetscRandom r;
00140 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&r);CHKERRQ(ierr);
00141 ierr = PetscRandomSetFromOptions(r);CHKERRQ(ierr);
00142
00143 for (int i=0;i<gdata->ln_vert;i++) {
00144 PetscReal value1,value2,value3,value4;
00145 PetscReal ran_num1,ran_num2,ran_num3;
00146 PetscReal ww;
00147 PetscReal xx1,xx2,xx3,xx4;
00148
00149 do {
00150 PetscRandomGetValue(r,&value1);
00151 PetscRandomGetValue(r,&value2);
00152
00153 xx1 = 2.0*value1-1.0;
00154 xx2 = 2.0*value2-1.0;
00155 ww = xx1*xx1+xx2*xx2;
00156 } while (ww >= 1.0);
00157
00158 ww=sqrt((-2.0*log(ww))/ww);
00159 ran_num1=xx1*ww;
00160 ran_num2=xx2*ww;
00161
00162 do {
00163 PetscRandomGetValue(r,&value3);
00164 PetscRandomGetValue(r,&value4);
00165
00166 xx3 = 2.0*value3-1.0;
00167 xx4 = 2.0*value4-1.0;
00168 ww = xx3*xx3+xx4*xx4;
00169 } while (ww >= 1.0);
00170
00171 ww=sqrt((-2.0*log(ww))/ww);
00172 ran_num3 = xx3*ww;
00173
00174 rval[ND*i+0]=(ran_num1+4.0)/(8.0);
00175 rval[ND*i+1]=(ran_num2+4.0)/(8.0);
00176 rval[ND*i+2]=(ran_num3+4.0)/(8.0);
00177 }
00178
00179 PetscRandomDestroy(r);
00180
00181
00182 MagparFunctionLogReturn(0);
00183 }
00184
00185
00186 int Helastic(GridData *gdata,Vec VHtotsum)
00187 {
00188 MagparFunctionInfoBegin;
00189
00190 if (doinit) {
00191 ierr = Helastic_Init(gdata);CHKERRQ(ierr);
00192 doinit=0;
00193 }
00194 if (!fieldon) {
00195 MagparFunctionInfoReturn(0);
00196 }
00197
00198 PetscReal *ta_Helastic,*ta_M;
00199 ierr = VecGetArray(VHelastic,&ta_Helastic);CHKERRQ(ierr);
00200 ierr = VecGetArray(gdata->M,&ta_M);CHKERRQ(ierr);
00201
00202
00203
00204
00205
00206 for (int j=0; j<gdata->ln_vert; j++) {
00207 int celltexture;
00208 PetscReal xA,xB,xC;
00209
00210 if (propdat[NPE*gdata->vertprop[j]+0]<=0.0) continue;
00211
00212
00213 celltexture=(int)propdat[NPE*gdata->vertprop[j]+0];
00214 switch (celltexture) {
00215 case 1:
00216 xA=0.0;
00217 xB=0.0;
00218 xC=0.0;
00219 break;
00220 case 2:
00221 xA=0.0;
00222 xB=rval[ND*j+0]*2.0*PETSC_PI;
00223 xC=rval[ND*j+1]*2.0*PETSC_PI;
00224 break;
00225 case 3:
00226 xA=PETSC_PI/2.0;
00227 xB=3.0*PETSC_PI/4.0;
00228 xC=0.0;
00229 break;
00230 case 4:
00231 xA=PETSC_PI/2.0;
00232 xB=3.0*PETSC_PI/4.0;
00233 xC=rval[ND*j+0]*2.0*PETSC_PI;
00234 break;
00235 case 5:
00236 xA=PETSC_PI/2.0;
00237 xB=asin(sqrt(2.0/3.0));
00238 xC=PETSC_PI;
00239 break;
00240 case 6:
00241 xA=asin(sqrt(2.0/3.0));
00242 xB=3*PETSC_PI/4.0;
00243 xC=rval[ND*j+0]*2.0*PETSC_PI;
00244 break;
00245 case 7:
00246 xA=rval[ND*j+0]*PETSC_PI;
00247 xB=rval[ND*j+1]*2.0*PETSC_PI;
00248 xC=rval[ND*j+2]*2.0*PETSC_PI;
00249 break;
00250 default:
00251 SETERRQ1(PETSC_ERR_ARG_INCOMP,
00252 "Unknown value for texture: %i\n",
00253 celltexture
00254 );
00255 break;
00256 }
00257
00258
00259 PetscReal c1,c2,c3;
00260 c1=cos(xA);
00261 c2=cos(xB);
00262 c3=cos(xC);
00263
00264 PetscReal s1,s2,s3;
00265 s1=sin(xA);
00266 s2=sin(xB);
00267 s3=sin(xC);
00268
00269
00270 PetscReal a11,a12,a13,a21,a22,a23,a31,a32,a33;
00271
00272 a11= c2*c3-c1*s2*s3;
00273 a21=-c2*s3-c1*s2*c3;
00274 a31=s1*s2;
00275 a12= s2*c3+c1*c2*s3;
00276 a22=-s2*s3+c1*c2*c3;
00277 a32=-s1*c2;
00278 a13=s1*s3;
00279 a23=s1*c3;
00280 a33=c1;
00281
00282
00283
00284 PetscReal detA;
00285 detA=a11*a22*a33+a13*a21*a32+a12*a23*a31-a11*a23*a32-a13*a22*a31-a12*a21*a33;
00286
00287
00288 PetscReal b11,b12,b13,b21,b22,b23,b31,b32,b33;
00289 b11 = (a22*a33 - a23*a32)/detA;
00290 b12 = (a23*a31 - a21*a33)/detA;
00291 b13 = (a21*a32 - a22*a31)/detA;
00292 b21 = (a13*a32 - a12*a33)/detA;
00293 b22 = (a11*a33 - a13*a31)/detA;
00294 b23 = (a12*a31 - a11*a32)/detA;
00295 b31 = (a12*a23 - a13*a22)/detA;
00296 b32 = (a13*a21 - a11*a23)/detA;
00297 b33 = (a11*a22 - a12*a21)/detA;
00298
00299
00300
00301 PetscReal mx,my,mz;
00302 mx=ta_M[ND*j+0];
00303 my=ta_M[ND*j+1];
00304 mz=ta_M[ND*j+2];
00305
00306 PetscReal mxp,myp,mzp;
00307 mxp=b11*mx+b21*my+b31*mz;
00308 myp=b12*mx+b22*my+b32*mz;
00309 mzp=b13*mx+b23*my+b33*mz;
00310
00311 PetscReal lam100,lam111;
00312 lam100=propdat[NPE*gdata->vertprop[j]+1];
00313 lam111=propdat[NPE*gdata->vertprop[j]+2];
00314
00315 PetscReal sigx,sigy,sigz;
00316 sigx=propdat[NPE*gdata->vertprop[j]+3];
00317 sigy=propdat[NPE*gdata->vertprop[j]+4];
00318 sigz=propdat[NPE*gdata->vertprop[j]+5];
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335 PetscReal msat;
00336 msat=gdata->propdat[NP*gdata->vertprop[j]+4];
00337
00338 PetscReal gamma1,gamma2,gamma3;
00339 PetscReal elasHkXp,elasHkYp,elasHkZp;
00340
00341 gamma1=b31;
00342 gamma2=b32;
00343 gamma3=b33;
00344 elasHkXp = 3.0/2.0*lam100*sigz/(msat)*(2.0*mxp*(gamma1*gamma1))+3.0*lam111*sigz/(msat)*(myp*gamma2+mzp*gamma3)*gamma1;
00345 elasHkYp = 3.0/2.0*lam100*sigz/(msat)*(2.0*myp*(gamma2*gamma2))+3.0*lam111*sigz/(msat)*(mzp*gamma3+mxp*gamma1)*gamma2;
00346 elasHkZp = 3.0/2.0*lam100*sigz/(msat)*(2.0*mzp*(gamma3*gamma3))+3.0*lam111*sigz/(msat)*(mxp*gamma1+myp*gamma2)*gamma3;
00347
00348 gamma1=b21;
00349 gamma2=b22;
00350 gamma3=b23;
00351 elasHkXp += 3.0/2.0*lam100*sigy/(msat)*(2.0*mxp*(gamma1*gamma1))+3.0*lam111*sigy/(msat)*(myp*gamma2+mzp*gamma3)*gamma1;
00352 elasHkYp += 3.0/2.0*lam100*sigy/(msat)*(2.0*myp*(gamma2*gamma2))+3.0*lam111*sigy/(msat)*(mzp*gamma3+mxp*gamma1)*gamma2;
00353 elasHkZp += 3.0/2.0*lam100*sigy/(msat)*(2.0*mzp*(gamma3*gamma3))+3.0*lam111*sigy/(msat)*(mxp*gamma1+myp*gamma2)*gamma3;
00354
00355 gamma1=b11;
00356 gamma2=b12;
00357 gamma3=b13;
00358 elasHkXp += 3.0/2.0*lam100*sigx/(msat)*(2.0*mxp*(gamma1*gamma1))+3.0*lam111*sigx/(msat)*(myp*gamma2+mzp*gamma3)*gamma1;
00359 elasHkYp += 3.0/2.0*lam100*sigx/(msat)*(2.0*myp*(gamma2*gamma2))+3.0*lam111*sigx/(msat)*(mzp*gamma3+mxp*gamma1)*gamma2;
00360 elasHkZp += 3.0/2.0*lam100*sigx/(msat)*(2.0*mzp*(gamma3*gamma3))+3.0*lam111*sigx/(msat)*(mxp*gamma1+myp*gamma2)*gamma3;
00361
00362
00363 PetscReal elasHkX,elasHkY,elasHkZ;
00364 elasHkX=a11*elasHkXp+a12*elasHkYp+a13*elasHkZp;
00365 elasHkY=a21*elasHkXp+a22*elasHkYp+a23*elasHkZp;
00366 elasHkZ=a31*elasHkXp+a32*elasHkYp+a33*elasHkZp;
00367
00368
00369 ta_Helastic[ND*j+0] = elasHkX/gdata->hscale*MU0;
00370 ta_Helastic[ND*j+1] = elasHkY/gdata->hscale*MU0;
00371 ta_Helastic[ND*j+2] = elasHkZ/gdata->hscale*MU0;
00372 }
00373
00374 ierr = VecRestoreArray(gdata->M,&ta_M);CHKERRQ(ierr);
00375 ierr = VecRestoreArray(VHelastic,&ta_Helastic);CHKERRQ(ierr);
00376
00377 ierr = VecAXPY(VHtotsum,1.0,VHelastic);CHKERRQ(ierr);
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387 MagparFunctionProfReturn(0);
00388 }
00389
00390
00391 int Helastic_Energy(GridData *gdata,Vec VMom,PetscReal *energy)
00392 {
00393 PetscReal e;
00394
00395 MagparFunctionInfoBegin;
00396
00397 if (!fieldon) {
00398 *energy=0.0;
00399 MagparFunctionInfoReturn(0);
00400 }
00401
00402 ierr = VecDot(VMom,VHelastic,&e);CHKERRQ(ierr);
00403 e *= -1.0;
00404
00405 *energy=e;
00406
00407 MagparFunctionProfReturn(0);
00408 }
00409