writelog_pid.c

Go to the documentation of this file.
00001 /*
00002     This file is part of magpar.
00003 
00004     Copyright (C) 2002-2010 Werner Scholz
00005 
00006     www:   http://www.magpar.net/
00007     email: magpar(at)magpar.net
00008 
00009     magpar is free software; you can redistribute it and/or modify
00010     it under the terms of the GNU General Public License as published by
00011     the Free Software Foundation; either version 2 of the License, or
00012     (at your option) any later version.
00013 
00014     magpar is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00017     GNU General Public License for more details.
00018 
00019     You should have received a copy of the GNU General Public License
00020     along with magpar; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
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     /* can use NORM_1 because all volumes are positive */
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     PetscReal v;
00137     VecStrideNorm(maskpidone[i],0,NORM_1,&v);
00138     PetscPrintf(PETSC_COMM_WORLD,"pid: %i  one: %g\n",i,v);
00139 */
00140   }
00141 
00142   if (rank) {
00143     MagparFunctionLogReturn(0);
00144   }
00145 
00146   for (int i=0;i<gdata->n_prop;i++) {
00147     /* skip if volume is zero */
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     /*      1    2    3    4    5    6    7    8    9   10   11   12   13 */
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     /*      1    2    3    4    5    6    7    8    9   10   11   12   13 */
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     /*      1    2    3    4    5    6    7    8    9   10 */
00209     ierr = PetscFPrintf(PETSC_COMM_WORLD,logfile[i],
00210       "#.%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n",
00211       /* "time",      */ "..........(ns)",
00212       /* "Etot",      */ "(J/m^3)",
00213       /* "Mavg_x",    */ "(1)",
00214       /* "Mavg_y",    */ "(1)",
00215       /* "Mavg_z",    */ "(1)",
00216       /* "|Mavg|",    */ "(1)",
00217       /* "Havg_x",    */ "(T)",
00218       /* "Havg_y",    */ "(T)",
00219       /* "Havg_z",    */ "(T)",
00220       /* "Edem",      */ "(J/m^3)",
00221       /* "Eexchani",  */ "(J/m^3)",
00222       /* "Eext",      */ "(J/m^3)",
00223 #ifdef EXCH
00224       /* "Eexchonly", */ "(J/m^3)"
00225 #else
00226       /* "void",      */ "-"
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     /* skip if volume is zero */
00255     if (volpid[i]<D_EPS) continue;
00256 
00257     ierr = VecPointwiseMult(Mpid,Mbak,maskpidvol[i]);CHKERRQ(ierr);
00258 
00259     /* calculate Mavg (ignoring Js) */
00260     PetscReal Mavg[ND];
00261     for (int j=0;j<ND;j++) {
00262       /* cannot use VecStrideNorm with NORM_1 because values are pos. and neg.! */
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     /* calculate Havg */
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 /*          1    2    3    4    5    6    7    8    9   10   11   12   13 */
00283     ierr = PetscFPrintf(PETSC_COMM_WORLD,logfile[i],
00284       "  %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e\n",
00285       /* "time",      "(ns)",    */ gdata->time*gdata->tscale*1e9,
00286       /* "Etot",      "(J/m^3)", */ gdata->Etot*gdata->escale*gdata->totvol/volpid[i],
00287       /* "Mavg_x",    "(1)",     */ Mavg[0],
00288       /* "Mavg_y",    "(1)",     */ Mavg[1],
00289       /* "Mavg_z",    "(1)",     */ Mavg[2],
00290       /* "|Mavg|",    "(1)",     */ sqrt(Mavg[0]*Mavg[0]+Mavg[1]*Mavg[1]+Mavg[2]*Mavg[2]),
00291       /* "Havg_x",    "(1)",     */ Havg[0],
00292       /* "Havg_y",    "(1)",     */ Havg[1],
00293       /* "Havg_z",    "(1)",     */ Havg[2],
00294       /* "Edem",      "(J/m^3)", */ gdata->Edem*gdata->escale*gdata->totvol/volpid[i],
00295 #ifdef EXCH
00296       /* "Eanionly",  "(J/m^3)", */ (gdata->Eexchani-gdata->Eexchonly)*gdata->escale*gdata->totvol/volpid[i],
00297 #else
00298       /* "Eexchani",  "(J/m^3)", */ gdata->Eexchani*gdata->escale*gdata->totvol/volpid[i],
00299 #endif
00300       /* "Eext",      "(J/m^3)", */ gdata->Eext*gdata->escale*gdata->totvol/volpid[i],
00301 #ifdef EXCH
00302       /* "Eexchonly", "(J/m^3)", */ gdata->Eexchonly*gdata->escale*gdata->totvol/volpid[i]
00303 #else
00304       /* "void",      "",        */ -1.0
00305 #endif
00306     );
00307   }
00308 
00309   /* recalculate energy with correct M (otherwise wrong results at end of simulation) */
00310   VecCopy(Mbak,gdata->M);
00311   Htot_Energy(gdata);
00312 
00313   MagparFunctionInfoReturn(0);
00314 }

magpar - Parallel Finite Element Micromagnetics Package
Copyright (C) 2002-2010 Werner Scholz