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: magset.c 2964 2010-02-05 19:31:03Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "init.h"
00028 #include "io/magpario.h"
00029 #include "util/util.h"
00030
00031 int MagSet(int pid,int init_mag,PetscReal init_magparm,GridData *gdata)
00032 {
00033 MagparFunctionLogBegin;
00034
00035 int rank,size;
00036 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00037 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00038
00039
00040 PetscPrintf(PETSC_COMM_WORLD,"pid: %i\n",pid+1);
00041 PetscPrintf(PETSC_COMM_WORLD,"init_mag: %i\n",init_mag);
00042 PetscPrintf(PETSC_COMM_WORLD,"init_magparm: %g\n",init_magparm);
00043
00044 if (init_mag==0) {
00045 char fmesh[256];
00046 PetscPrintf(PETSC_COMM_WORLD,"inp file to be read: %s.%04d.inp\n",gdata->simname,gdata->inp);
00047 ierr = PetscSNPrintf(fmesh,255,"%s.%04d.inp",gdata->simname,gdata->inp);CHKERRQ(ierr);
00048 gdata->inp++;
00049
00050 ierr = ReadINP(gdata,fmesh,gdata->M,1);CHKERRQ(ierr);
00051 MagparFunctionLogReturn(0);
00052 }
00053 else if (PetscAbsInt(init_mag)==14) {
00054 #ifdef PYTHON
00055 if (!rank) {
00056
00057
00058 if(!gdata->pModule) {
00059 SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Cannot execute Python function. Module %s has not been imported!\n",MAGPARSTDPYTHONMOD);
00060 }
00061
00062 #define INITMAG "initmag"
00063 PyObject *pFunc;
00064 pFunc=PyObject_GetAttrString(gdata->pModule,INITMAG);
00065 if(!pFunc) {
00066 SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Cannot find Python function %s!\n",INITMAG);
00067 }
00068 if (!PyCallable_Check(pFunc)) {
00069 SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Python object %s is not callable (not a function)!\n",INITMAG);
00070 }
00071
00072
00073 PyObject *pTuple;
00074 pTuple = PyTuple_New(ND);
00075
00076 PyObject *pVal=PETSC_NULL;
00077
00078
00079
00080
00081
00082
00083 PetscReal *ta_M;
00084 ierr = VecGetArray(gdata->M,&ta_M);CHKERRQ(ierr);
00085
00086
00087 for (int j=0; j<gdata->n_vert; j++) {
00088 PetscReal x,y,z;
00089 x=gdata->vertxyz[ND*j+0];
00090 y=gdata->vertxyz[ND*j+1];
00091 z=gdata->vertxyz[ND*j+2];
00092
00093 PyTuple_SetItem(pTuple,0,PyFloat_FromDouble(x));
00094 PyTuple_SetItem(pTuple,1,PyFloat_FromDouble(y));
00095 PyTuple_SetItem(pTuple,2,PyFloat_FromDouble(z));
00096
00097
00098 pVal = PyObject_CallObject(pFunc,pTuple);
00099
00100
00101 for (int i=0; i<ND; i++) {
00102 ta_M[ND*j+i]=PyFloat_AsDouble(PyTuple_GetItem(pVal,i));
00103 }
00104 }
00105
00106 ierr = VecRestoreArray(gdata->M,&ta_M);CHKERRQ(ierr);
00107
00108
00109 Py_DECREF(pFunc);
00110 Py_DECREF(pTuple);
00111 Py_DECREF(pVal);
00112 }
00113 #else
00114 SETERRQ(PETSC_ERR_ARG_CORRUPT,"magpar not linked with Python interpreter. Cannot execute Python scripts!\n");
00115 #endif
00116 }
00117
00118 if (rank) {
00119 VecAssemblyBegin(gdata->M);
00120 VecAssemblyEnd(gdata->M);
00121
00122
00123 RenormVec(gdata->M,1.0,0.0,ND);
00124
00125 MagparFunctionLogReturn(0);
00126 }
00127 else {
00128 PetscRandom rctx;
00129 ierr = PetscRandomCreate(PETSC_COMM_SELF,&rctx);CHKERRQ(ierr);
00130 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
00131
00132 ierr = VertProp(gdata);CHKERRQ(ierr);
00133
00134 for (int i=0; i<gdata->n_vert; i++) {
00135 if (pid>=0 && gdata->vertprop[i]!=pid) continue;
00136 switch(PetscAbsInt(init_mag)) {
00137 case 0:
00138
00139 break;
00140 case 1:
00141 ierr = VecSetValue(gdata->M,ND*i+0,1.0*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00142 ierr = VecSetValue(gdata->M,ND*i+1,0.0,INSERT_VALUES);CHKERRQ(ierr);
00143 ierr = VecSetValue(gdata->M,ND*i+2,0.0,INSERT_VALUES);CHKERRQ(ierr);
00144 break;
00145 case 2:
00146 ierr = VecSetValue(gdata->M,ND*i+0,0.0,INSERT_VALUES);CHKERRQ(ierr);
00147 ierr = VecSetValue(gdata->M,ND*i+1,1.0*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00148 ierr = VecSetValue(gdata->M,ND*i+2,0.0,INSERT_VALUES);CHKERRQ(ierr);
00149 break;
00150 case 3:
00151 ierr = VecSetValue(gdata->M,ND*i+0,0.0,INSERT_VALUES);CHKERRQ(ierr);
00152 ierr = VecSetValue(gdata->M,ND*i+1,0.0,INSERT_VALUES);CHKERRQ(ierr);
00153 ierr = VecSetValue(gdata->M,ND*i+2,1.0*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00154 break;
00155 case 4:
00156 ierr = VecSetValue(gdata->M,ND*i+0,0.57735027*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00157 ierr = VecSetValue(gdata->M,ND*i+1,0.57735027*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00158 ierr = VecSetValue(gdata->M,ND*i+2,0.57735027*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00159 break;
00160 case 5:
00161 ierr = VecSetValue(gdata->M,ND*i+0,(gdata->vertxyz[ND*i+0]-init_magparm)*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00162 ierr = VecSetValue(gdata->M,ND*i+1,(gdata->vertxyz[ND*i+1]-init_magparm)*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00163 ierr = VecSetValue(gdata->M,ND*i+2,1.0*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00164 break;
00165 case 6:
00166 ierr = VecSetValue(gdata->M,ND*i+0,sin(init_magparm)*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00167 ierr = VecSetValue(gdata->M,ND*i+1,0.0,INSERT_VALUES);CHKERRQ(ierr);
00168 ierr = VecSetValue(gdata->M,ND*i+2,cos(init_magparm)*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00169 break;
00170 case 7:
00171
00172
00173
00174
00175 {
00176 PetscReal x, y, r, R, A;
00177 PetscReal modmx,modmy,modmz;
00178
00179 x = gdata->vertxyz[ND*i+0];
00180 y = gdata->vertxyz[ND*i+1];
00181 r = sqrt(x*x + y*y);
00182 R = init_magparm;
00183
00184 A = 2.0*R*r/(R*R+r*r);
00185 if (r <= fabs(R)) {
00186 if (r > 1e-10) {
00187 modmx = A*y/r;
00188 modmy = -A*x/r;
00189 modmz = -sqrt(1.0-A*A);
00190 } else {
00191 modmx=0.0;
00192 modmy=0.0;
00193 modmz=-1.0;
00194 }
00195 } else {
00196 modmx=y/r;
00197 modmy=-x/r;
00198 modmz=0.0;
00199 }
00200 ierr = VecSetValue(gdata->M,ND*i+0,modmx*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00201 ierr = VecSetValue(gdata->M,ND*i+1,modmy*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00202 ierr = VecSetValue(gdata->M,ND*i+2,modmz*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00203 }
00204 break;
00205 case 8:
00206 PetscReal rval[ND],rnorm;
00207
00208 PetscRandomGetValue(rctx,&rval[0]);
00209 PetscRandomGetValue(rctx,&rval[1]);
00210 PetscRandomGetValue(rctx,&rval[2]);
00211 rnorm=my_dnrm2(ND,rval,1);
00212
00213 for (int k=0; k<ND; k++) {
00214 ierr = VecSetValue(
00215 gdata->M,
00216 ND*i+k,
00217 (rval[k]/rnorm*2.0-1.0)*PetscSign(init_mag),
00218 INSERT_VALUES
00219 );CHKERRQ(ierr);
00220 }
00221 break;
00222 case 9:
00223 if (gdata->vertxyz[ND*i+0] < init_magparm*0.95) {
00224 ierr = VecSetValue(gdata->M,ND*i+0, 0.0,INSERT_VALUES);CHKERRQ(ierr);
00225 ierr = VecSetValue(gdata->M,ND*i+1, 0.0,INSERT_VALUES);CHKERRQ(ierr);
00226 ierr = VecSetValue(gdata->M,ND*i+2,+1.0*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00227 }
00228 else if (gdata->vertxyz[ND*i+0] > init_magparm*1.05) {
00229 ierr = VecSetValue(gdata->M,ND*i+0, 0.0,INSERT_VALUES);CHKERRQ(ierr);
00230 ierr = VecSetValue(gdata->M,ND*i+1, 0.0,INSERT_VALUES);CHKERRQ(ierr);
00231 ierr = VecSetValue(gdata->M,ND*i+2,-1.0*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00232 }
00233 else {
00234 ierr = VecSetValue(gdata->M,ND*i+0, 0.1*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00235 ierr = VecSetValue(gdata->M,ND*i+1,-0.9*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00236 ierr = VecSetValue(gdata->M,ND*i+2, 0.1*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00237 }
00238 break;
00239 case 10:
00240 for (int k=0; k<ND; k++) {
00241 ierr = VecSetValue(
00242 gdata->M,
00243 ND*i+k,
00244 gdata->propdat[NP*gdata->vertprop[i]+6+k]*PetscSign(init_mag),
00245 INSERT_VALUES
00246 );CHKERRQ(ierr);
00247 }
00248 break;
00249 case 11:
00250 ierr = VecSetValue(gdata->M,ND*i+0,cos(init_magparm*PETSC_PI/180.0)*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00251 ierr = VecSetValue(gdata->M,ND*i+1,sin(init_magparm*PETSC_PI/180.0)*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00252 ierr = VecSetValue(gdata->M,ND*i+2,0.0,INSERT_VALUES);CHKERRQ(ierr);
00253 break;
00254 case 12:
00255 if (gdata->vertxyz[ND*i+0] < init_magparm*0.95) {
00256 ierr=VecSetValue(gdata->M,ND*i+0, 1.0*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00257 ierr=VecSetValue(gdata->M,ND*i+1, 0.0,INSERT_VALUES);CHKERRQ(ierr);
00258 ierr=VecSetValue(gdata->M,ND*i+2, 0.0,INSERT_VALUES);CHKERRQ(ierr);
00259 }
00260 else if (gdata->vertxyz[ND*i+0] > init_magparm*1.05) {
00261 ierr=VecSetValue(gdata->M,ND*i+0,-1.0*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00262 ierr=VecSetValue(gdata->M,ND*i+1, 0.0,INSERT_VALUES);CHKERRQ(ierr);
00263 ierr=VecSetValue(gdata->M,ND*i+2, 0.0,INSERT_VALUES);CHKERRQ(ierr);
00264 }
00265 else {
00266 ierr=VecSetValue(gdata->M,ND*i+0, 0.0,INSERT_VALUES);CHKERRQ(ierr);
00267 ierr=VecSetValue(gdata->M,ND*i+1, 1.0*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00268 ierr=VecSetValue(gdata->M,ND*i+2, 0.0,INSERT_VALUES);CHKERRQ(ierr);
00269 }
00270 break;
00271 case 13:
00272
00273
00274
00275
00276 {
00277 PetscReal x, y, r, R, A;
00278 PetscReal modmx,modmy,modmz;
00279
00280 x = gdata->vertxyz[ND*i+0];
00281 y = gdata->vertxyz[ND*i+1];
00282 r = sqrt(x*x + y*y);
00283 R = init_magparm;
00284
00285 A = 2.0*R*r/(R*R+r*r);
00286 if (r <= fabs(R)) {
00287 if (r > 1e-10) {
00288 modmx = A*y/r;
00289 modmy = -A*x/r;
00290 modmz = -sqrt(1.0-A*A);
00291 } else {
00292 modmx = 0.0;
00293 modmy = 0.0;
00294 modmz = -1.0;
00295 }
00296 } else {
00297 modmx = (x<0.0 ? +1.0 : -1.0);
00298 modmy = 0.0;
00299 modmz = 0.0;
00300 }
00301 ierr = VecSetValue(gdata->M,ND*i+0,modmx*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00302 ierr = VecSetValue(gdata->M,ND*i+1,modmy*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00303 ierr = VecSetValue(gdata->M,ND*i+2,modmz*PetscSign(init_mag),INSERT_VALUES);CHKERRQ(ierr);
00304 }
00305 break;
00306 case 14:
00307
00308 break;
00309 default:
00310 SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"init_mag value %i is illegal!\n",PetscAbsInt(init_mag));
00311 break;
00312 }
00313 }
00314
00315 ierr = PetscRandomDestroy(rctx);CHKERRQ(ierr);
00316 ierr = PetscFree(gdata->vertprop);CHKERRQ(ierr);
00317 gdata->vertprop=PETSC_NULL;
00318 }
00319
00320 ierr = VecAssemblyBegin(gdata->M);CHKERRQ(ierr);
00321 ierr = VecAssemblyEnd(gdata->M);CHKERRQ(ierr);
00322
00323
00324 RenormVec(gdata->M,1.0,0.0,ND);
00325
00326 MagparFunctionLogReturn(0);
00327 }