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: writelog_pid.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 #include "util/util.h"
00029 #include "field/field.h"
00030
00031 #ifdef ADDONS
00032 #include "addons/addons.h"
00033 #endif
00034
00035 #ifdef SUNDIALS_VERSION
00036 #include "llg/llg.h"
00037 #endif
00038
00039 static int doinit=1;
00040 static int fieldon=0;
00041
00042 static FILE **logfile;
00043 static Vec Mpid,Mbak;
00044 static Vec Hpid;
00045 static Vec Xcomp;
00046 static Vec *maskpidvol;
00047 static Vec *maskpidone;
00048 static PetscReal *volpid;
00049
00050 int WriteLogPidInit(GridData *gdata)
00051 {
00052 MagparFunctionLogBegin;
00053
00054 int rank,size;
00055 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00056 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00057
00058 PetscTruth flg;
00059 ierr = PetscOptionsGetInt(PETSC_NULL,"-logpid",(PetscInt*)&fieldon,&flg);CHKERRQ(ierr);
00060 if (flg!=PETSC_TRUE) {
00061 fieldon=0;
00062 PetscPrintf(PETSC_COMM_WORLD,
00063 "Option -logpid not found, using default value: %i\n",
00064 fieldon
00065 );
00066 }
00067
00068 if (gdata->n_prop==1) {
00069 PetscPrintf(PETSC_COMM_WORLD,"Only %i pid found\n",gdata->n_prop);
00070 fieldon=0;
00071 }
00072
00073 if (!fieldon) {
00074 PetscPrintf(PETSC_COMM_WORLD,"Not storing *.log_pid files\n");
00075 MagparFunctionLogReturn(0);
00076 }
00077
00078 PetscPrintf(PETSC_COMM_WORLD,"Creating *.log_pid files\n");
00079
00080 ierr = VecDuplicate(gdata->M,&Mpid);CHKERRQ(ierr);
00081 ierr = VecDuplicate(gdata->M,&Mbak);CHKERRQ(ierr);
00082
00083 ierr = VecDuplicate(gdata->M,&Hpid);CHKERRQ(ierr);
00084
00085 ierr = VecCreate(PETSC_COMM_WORLD,&Xcomp);CHKERRQ(ierr);
00086 ierr = VecSetSizes(Xcomp,gdata->ln_vert,gdata->n_vert);CHKERRQ(ierr);
00087 ierr = VecSetFromOptions(Xcomp);CHKERRQ(ierr);
00088
00089 ierr = PetscMalloc(gdata->n_prop*sizeof(FILE*),&logfile);CHKERRQ(ierr);
00090 ierr = PetscMalloc(gdata->n_prop*sizeof(Vec),&maskpidvol);CHKERRQ(ierr);
00091 for (int i=0;i<gdata->n_prop;i++) {
00092 ierr = VecDuplicate(gdata->M,&maskpidvol[i]);CHKERRQ(ierr);
00093 }
00094 PetscReal *ta_evol,*ta_vvol;
00095 VecGetArray(gdata->elevol,&ta_evol);
00096 for (int i=0;i<gdata->ln_ele;i++) {
00097 for (int l=0;l<NV;l++) {
00098 for (int k=0;k<ND;k++) {
00099 ierr = VecSetValue(
00100 maskpidvol[gdata->eleprop[i]],
00101 ND*gdata->elevert[NV*i+l]+k,
00102 ta_evol[i]/NV,
00103 ADD_VALUES
00104 );CHKERRQ(ierr);
00105 }
00106 }
00107 }
00108 VecRestoreArray(gdata->elevol,&ta_evol);
00109
00110 ierr = PetscMalloc(gdata->n_prop*sizeof(PetscReal),&volpid);CHKERRQ(ierr);
00111 for (int i=0;i<gdata->n_prop;i++) {
00112 ierr = VecAssemblyBegin(maskpidvol[i]);CHKERRQ(ierr);
00113 ierr = VecAssemblyEnd(maskpidvol[i]);CHKERRQ(ierr);
00114
00115 VecStrideNorm(maskpidvol[i],0,NORM_1,&volpid[i]);
00116 PetscPrintf(PETSC_COMM_WORLD,"pid: %i volume: %g\n",i,volpid[i]);
00117 }
00118
00119 ierr = PetscMalloc(gdata->n_prop*sizeof(Vec),&maskpidone);CHKERRQ(ierr);
00120 for (int i=0;i<gdata->n_prop;i++) {
00121 PetscReal *ta_mpv;
00122
00123 ierr = VecDuplicate(gdata->M,&maskpidone[i]);CHKERRQ(ierr);
00124 ierr = VecCopy(maskpidvol[i],maskpidone[i]);CHKERRQ(ierr);
00125 VecGetArray(maskpidone[i],&ta_mpv);
00126 VecGetArray(gdata->vertvol,&ta_vvol);
00127 for (int j=0;j<gdata->ln_vert;j++) {
00128 for (int k=0;k<ND;k++) {
00129 ta_mpv[ND*j+k]/=ta_vvol[j];
00130 }
00131 }
00132 VecRestoreArray(maskpidone[i],&ta_mpv);
00133 VecRestoreArray(gdata->vertvol,&ta_vvol);
00134
00135
00136
00137
00138
00139
00140 }
00141
00142 if (rank) {
00143 MagparFunctionLogReturn(0);
00144 }
00145
00146 for (int i=0;i<gdata->n_prop;i++) {
00147
00148 if (volpid[i]<D_EPS) {
00149 PetscPrintf(PETSC_COMM_WORLD,"volume of pid: %i is zero, skipping output\n",i);
00150 continue;
00151 }
00152
00153 char fmesh[256];
00154 ierr = PetscSNPrintf(fmesh,255,"%s%s%03i",gdata->simname,".log_",i+1);CHKERRQ(ierr);
00155 ierr = PetscFOpen(PETSC_COMM_WORLD,fmesh,"a",&logfile[i]);CHKERRQ(ierr);
00156 if (!logfile[i]) {
00157 SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",fmesh);
00158 }
00159 PetscPrintf(PETSC_COMM_WORLD,"Opening file %s\n",fmesh);
00160
00161 PetscGetDate(fmesh,255);
00162 ierr = PetscFPrintf(PETSC_COMM_WORLD,logfile[i],"#.date: %s\n",fmesh);
00163 ierr = PetscFPrintf(PETSC_COMM_WORLD,logfile[i],"#.volume_with_pid:%i\n",i+1);
00164 ierr = PetscFPrintf(PETSC_COMM_WORLD,logfile[i],"#.volume:%g\n",volpid[i]);
00165
00166
00167 ierr = PetscFPrintf(PETSC_COMM_WORLD,logfile[i],
00168 "#.%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n",
00169 "............1:",
00170 " 2:",
00171 " 3:",
00172 " 4:",
00173 " 5:",
00174 " 6:",
00175 " 7:",
00176 " 8:",
00177 " 9:",
00178 "10:",
00179 "11:",
00180 "12:",
00181 "13:"
00182 );CHKERRQ(ierr);
00183
00184 ierr = PetscFPrintf(PETSC_COMM_WORLD,logfile[i],
00185 "#.%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n",
00186 "..........time",
00187 "Etot",
00188 "Mavg_x",
00189 "Mavg_y",
00190 "Mavg_z",
00191 "|Mavg|",
00192 "Havg_x",
00193 "Havg_y",
00194 "Havg_z",
00195 "Edem",
00196 #ifdef EXCH
00197 "Eanionly",
00198 #else
00199 "Eexchani",
00200 #endif
00201 "Eext",
00202 #ifdef EXCH
00203 "Eexchonly"
00204 #else
00205 "-"
00206 #endif
00207 );CHKERRQ(ierr);
00208
00209 ierr = PetscFPrintf(PETSC_COMM_WORLD,logfile[i],
00210 "#.%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n",
00211 "..........(ns)",
00212 "(J/m^3)",
00213 "(1)",
00214 "(1)",
00215 "(1)",
00216 "(1)",
00217 "(T)",
00218 "(T)",
00219 "(T)",
00220 "(J/m^3)",
00221 "(J/m^3)",
00222 "(J/m^3)",
00223 #ifdef EXCH
00224 "(J/m^3)"
00225 #else
00226 "-"
00227 #endif
00228 );CHKERRQ(ierr);
00229 }
00230
00231 MagparFunctionLogReturn(0);
00232 }
00233
00234
00235 int WriteLogPid(GridData *gdata)
00236 {
00237 MagparFunctionInfoBegin;
00238
00239 int rank,size;
00240 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00241 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00242
00243 if (doinit) {
00244 ierr = WriteLogPidInit(gdata);CHKERRQ(ierr);
00245 doinit=0;
00246 }
00247 if (!fieldon) {
00248 MagparFunctionInfoReturn(0);
00249 }
00250
00251 VecCopy(gdata->M,Mbak);
00252
00253 for (int i=0;i<gdata->n_prop;i++) {
00254
00255 if (volpid[i]<D_EPS) continue;
00256
00257 ierr = VecPointwiseMult(Mpid,Mbak,maskpidvol[i]);CHKERRQ(ierr);
00258
00259
00260 PetscReal Mavg[ND];
00261 for (int j=0;j<ND;j++) {
00262
00263 ierr = VecStrideGather(Mpid,j,Xcomp,INSERT_VALUES);CHKERRQ(ierr);
00264 ierr = VecSum(Xcomp,&Mavg[j]);CHKERRQ(ierr);
00265 }
00266 my_dscal(ND,1.0/volpid[i],Mavg,1);
00267
00268 ierr = VecPointwiseMult(Hpid,gdata->VHtot,maskpidvol[i]);CHKERRQ(ierr);
00269
00270
00271 PetscReal Havg[ND];
00272 for (int j=0;j<ND;j++) {
00273 ierr = VecStrideGather(Hpid,j,Xcomp,INSERT_VALUES);CHKERRQ(ierr);
00274 ierr = VecSum(Xcomp,&Havg[j]);CHKERRQ(ierr);
00275 }
00276 my_dscal(ND,1.0/volpid[i],Havg,1);
00277
00278 ierr = VecPointwiseMult(Mpid,Mbak,maskpidone[i]);CHKERRQ(ierr);
00279 VecCopy(Mpid,gdata->M);
00280 Htot_Energy(gdata);
00281
00282
00283 ierr = PetscFPrintf(PETSC_COMM_WORLD,logfile[i],
00284 " %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e\n",
00285 gdata->time*gdata->tscale*1e9,
00286 gdata->Etot*gdata->escale*gdata->totvol/volpid[i],
00287 Mavg[0],
00288 Mavg[1],
00289 Mavg[2],
00290 sqrt(Mavg[0]*Mavg[0]+Mavg[1]*Mavg[1]+Mavg[2]*Mavg[2]),
00291 Havg[0],
00292 Havg[1],
00293 Havg[2],
00294 gdata->Edem*gdata->escale*gdata->totvol/volpid[i],
00295 #ifdef EXCH
00296 (gdata->Eexchani-gdata->Eexchonly)*gdata->escale*gdata->totvol/volpid[i],
00297 #else
00298 gdata->Eexchani*gdata->escale*gdata->totvol/volpid[i],
00299 #endif
00300 gdata->Eext*gdata->escale*gdata->totvol/volpid[i],
00301 #ifdef EXCH
00302 gdata->Eexchonly*gdata->escale*gdata->totvol/volpid[i]
00303 #else
00304 -1.0
00305 #endif
00306 );
00307 }
00308
00309
00310 VecCopy(Mbak,gdata->M);
00311 Htot_Energy(gdata);
00312
00313 MagparFunctionInfoReturn(0);
00314 }