[a0bcf1] | 1 | /** \file output.c
|
---|
| 2 | * Output of forces and energies, Visuals (density and ions for OpenDX).
|
---|
| 3 | *
|
---|
| 4 | * Herein all the functions concerning the output of data are gathered:
|
---|
| 5 | * Initialization of the FileData structure InitOutVisArray() and the files themselves
|
---|
| 6 | * InitOutputFiles(),
|
---|
| 7 | * Saving of a calculated state OutputVisSrcFiles() - 1. Psis OutputSrcPsiDensity(), 2. Ions
|
---|
| 8 | * OutSrcIons() - or retrieving Psis ReadSrcPsiDensity() and Ions ReadSrcIons(),
|
---|
| 9 | * Preparation (in case of RiemannTensor use) CalculateOutVisDensityPos(), OutVisPosRTransformPosNFRto0(),
|
---|
| 10 | * Output of visual data (for OpenDx explorer) OutputVis() - 1. OpenDX files CreateDensityOutputGeneral(),
|
---|
| 11 | * 2. electronic densitiy data OutVisDensity() (uses MPI sending density coefficients OutputOutVisDensity()
|
---|
| 12 | * and receiving and writing CombineOutVisDensity()), 3. ion data OutVisIons() -
|
---|
| 13 | * Closing of files CloseOutputFiles().
|
---|
| 14 | *
|
---|
| 15 | * There are some more routines: OutputCurrentDensity() outputs the current density for each magnetic field
|
---|
| 16 | * direction, OutputVisAllOrbital() outputs all orbitals of a certain minimsation type and
|
---|
| 17 | * TestReadnWriteSrcDensity() checks whether a certain minimisation group can be written and read again
|
---|
| 18 | * correctly.
|
---|
| 19 | *
|
---|
| 20 | * There are some helpers that open files with certain filenames, making extensive usage of all
|
---|
| 21 | * the suffixes defined in here: OpenFile(), OpenFileNo(), OpenFileNo2(), OpenFileNoNo() and
|
---|
| 22 | * OpenFileNoPost().
|
---|
| 23 | *
|
---|
| 24 | *
|
---|
| 25 | Project: ParallelCarParrinello
|
---|
| 26 | \author Jan Hamaekers, Frederik Heber
|
---|
| 27 | \date 2000, 2006
|
---|
| 28 |
|
---|
| 29 | File: output.c
|
---|
| 30 | $Id: output.c,v 1.51.2.2 2007-04-21 12:55:50 foo Exp $
|
---|
| 31 | */
|
---|
| 32 | #include<stdlib.h>
|
---|
| 33 | #include<stdio.h>
|
---|
| 34 | #include<string.h>
|
---|
| 35 | #include<math.h>
|
---|
| 36 | //#include<sys/time.h>
|
---|
| 37 | #include <time.h>
|
---|
| 38 | #include<unistd.h>
|
---|
| 39 | #include"data.h"
|
---|
| 40 | #include"density.h"
|
---|
| 41 | #include"errors.h"
|
---|
| 42 | #include"gramsch.h"
|
---|
| 43 | #include"helpers.h"
|
---|
| 44 | #include "init.h"
|
---|
| 45 | #include"myfft.h"
|
---|
| 46 | #include"mymath.h"
|
---|
| 47 | #include"output.h"
|
---|
| 48 | #include"pdbformat.h"
|
---|
| 49 | #include"perturbed.h"
|
---|
| 50 |
|
---|
| 51 |
|
---|
| 52 | /* Konvention: Rueckgabe 0 einer Funktion, bedeutet keinen Fehler (entsprechend exitcode 0) */
|
---|
| 53 | /* Oeffnet Datei P->mainname+"..." mit what*/
|
---|
| 54 |
|
---|
| 55 | /** Opens a file with FileData::mainname+suffix.
|
---|
| 56 | * Both the path under which the main programme resides and the path to the config file are
|
---|
| 57 | * tried subsequently.
|
---|
| 58 | * \param *P Problem at hand
|
---|
| 59 | * \param **file file pointer array
|
---|
| 60 | * \param *suffix suffix after mainname
|
---|
| 61 | * \param *what access parameter for the file: r, w, rw, r+, w+ ...
|
---|
| 62 | * \param verbose 1 - print status, 0 - don't print anything
|
---|
| 63 | * \return 0 - could not open file, 1 - file is open
|
---|
| 64 | */
|
---|
| 65 | int OpenFile(struct Problem *P, FILE** file, const char* suffix, const char* what, int verbose)
|
---|
| 66 | {
|
---|
| 67 | char* name; /* Zu erzeugender Dateiname */
|
---|
| 68 | name = (char*)
|
---|
| 69 | Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(suffix) + 1,"OpenFile");
|
---|
| 70 | sprintf(name, "%s%s%s", P->Files.default_path, P->Files.mainname, suffix);
|
---|
| 71 | *file = fopen(name, what);
|
---|
| 72 | if (*file == 0) {
|
---|
| 73 | fprintf(stderr,"(%i) Normal access failed: name %s, what %s\n", P->Par.me, name, what);
|
---|
| 74 | /* hier default file ausprobieren, falls nur gelesen werden soll! */
|
---|
| 75 | if(*what == 'r') {
|
---|
| 76 | name = (char*)
|
---|
| 77 | Realloc(name,strlen(P->Files.default_path) + strlen(suffix)+1,"OpenFile");
|
---|
| 78 | sprintf(name, "%s%s", P->Files.default_path,suffix);
|
---|
| 79 | *file = fopen(name,what);
|
---|
| 80 | if (*file != 0) {
|
---|
| 81 | if (verbose) fprintf(stderr,"(%i) Default file is open: %s\n",P->Par.me,name);
|
---|
[64fa9e] | 82 | Free(name, "OpenFile: name");
|
---|
[a0bcf1] | 83 | return(1);
|
---|
| 84 | }
|
---|
| 85 | }
|
---|
| 86 | fprintf(stderr,"\n(%i)Error: Cannot open neither normal nor default file: %s\n",P->Par.me,name);
|
---|
[64fa9e] | 87 | Free(name, "OpenFile: name");
|
---|
[a0bcf1] | 88 | return(0);
|
---|
| 89 | } else {
|
---|
| 90 | if(verbose) fprintf(stderr,"(%i) File is open: %s\n",P->Par.me,name);
|
---|
[64fa9e] | 91 | Free(name, "OpenFile: name");
|
---|
[a0bcf1] | 92 | return(1);
|
---|
| 93 | }
|
---|
| 94 | }
|
---|
| 95 |
|
---|
| 96 | /** Opens a file with FileData::mainname+suffix+"."+No (2 digits).
|
---|
| 97 | * Both the path under which the main programme resides and the path to the config file are
|
---|
| 98 | * tried subsequently.
|
---|
| 99 | * \param *P Problem at hand
|
---|
| 100 | * \param **file file pointer array
|
---|
| 101 | * \param *suffix suffix after mainname
|
---|
| 102 | * \param No the number with up to two digits
|
---|
| 103 | * \param *what access parameter for the file: r, w, rw, r+, w+ ...
|
---|
| 104 | * \param verbose 1 - print status, 0 - don't print anything
|
---|
| 105 | * \return 0 - could not open file, 1 - file is open
|
---|
| 106 | */
|
---|
| 107 | int OpenFileNo2(struct Problem *P, FILE** file, const char* suffix, int No, const char* what, int verbose)
|
---|
| 108 | {
|
---|
| 109 | char* name; /* Zu erzeugender Dateiname */
|
---|
| 110 | name = (char*)
|
---|
| 111 | Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(suffix) + 4,"OpenFile");
|
---|
| 112 | sprintf(name, "%s%s%s.%02i", P->Files.default_path, P->Files.mainname, suffix, No);
|
---|
| 113 | *file = fopen(name, what);
|
---|
| 114 | if (*file == 0) {
|
---|
| 115 | /* falls nur gelesen wird, auch default_path ausprobieren */
|
---|
| 116 | if (*what == 'r') {
|
---|
| 117 | name = (char*)
|
---|
| 118 | Realloc(name,strlen(P->Files.default_path) + strlen(suffix) + 4,"OpenFileNo2");
|
---|
| 119 | sprintf(name, "%s%s.%02i", P->Files.default_path, suffix, No);
|
---|
| 120 | *file = fopen(name, what);
|
---|
| 121 | if (*file != 0) {
|
---|
| 122 | if(verbose) fprintf(stderr,"(%i) Default file is open: %s\n", P->Par.me, name);
|
---|
[64fa9e] | 123 | Free(name, "OpenFileNo2: name");
|
---|
[a0bcf1] | 124 | return(1);
|
---|
| 125 | }
|
---|
| 126 | }
|
---|
| 127 | fprintf(stderr,"\n(%i)Error: Cannot open file: %s\n", P->Par.me, name);
|
---|
[64fa9e] | 128 | Free(name, "OpenFileNo2: name");
|
---|
[a0bcf1] | 129 | return(0);
|
---|
| 130 | } else {
|
---|
| 131 | if(verbose) fprintf(stderr,"(%i) File is open: %s\n",P->Par.me, name);
|
---|
[64fa9e] | 132 | Free(name, "OpenFileNo2: name");
|
---|
[a0bcf1] | 133 | return(1);
|
---|
| 134 | }
|
---|
| 135 | }
|
---|
| 136 |
|
---|
| 137 | /* Oeffnet Datei P->Files.mainname+"...".Nr(4stellig) mit what*/
|
---|
| 138 | /** Opens a file with FileData::mainname+suffix+"."+No (4 digits).
|
---|
| 139 | * Both the path under which the main programme resides and the path to the config file are
|
---|
| 140 | * tried subsequently.
|
---|
| 141 | * \param *P Problem at hand
|
---|
| 142 | * \param **file file pointer array
|
---|
| 143 | * \param *suffix suffix after mainname
|
---|
| 144 | * \param No the number with up to four digits
|
---|
| 145 | * \param *what access parameter for the file: r, w, rw, r+, w+ ...
|
---|
| 146 | * \param verbose 1 - print status, 0 - don't print anything
|
---|
| 147 | * \return 0 - could not open file, 1 - file is open
|
---|
| 148 | */
|
---|
| 149 | int OpenFileNo(struct Problem *P, FILE** file, const char* suffix, int No, const char* what, int verbose)
|
---|
| 150 | {
|
---|
| 151 | char* name; /* Zu erzeugender Dateiname */
|
---|
| 152 | name = (char*)
|
---|
| 153 | Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(suffix) + 6,"OpenFileNo");
|
---|
| 154 | sprintf(name, "%s%s%s.%04i", P->Files.default_path, P->Files.mainname, suffix, No);
|
---|
| 155 | *file = fopen(name, what);
|
---|
| 156 | if (*file == 0) {
|
---|
| 157 | /* falls nur gelesen wird, auch default_path ausprobieren */
|
---|
| 158 | if (*what == 'r') {
|
---|
| 159 | name = (char*)
|
---|
| 160 | Realloc(name,strlen(P->Files.default_path) + strlen(suffix) + 6,"OpenFileNo");
|
---|
| 161 | sprintf(name, "%s%s.%04i", P->Files.default_path, suffix, No);
|
---|
| 162 | *file = fopen(name, what);
|
---|
| 163 | if (*file != 0) {
|
---|
| 164 | if(verbose) fprintf(stderr,"(%i) Default file is open: %s\n", P->Par.me, name);
|
---|
[64fa9e] | 165 | Free(name, "OpenFileNo: name");
|
---|
[a0bcf1] | 166 | return(1);
|
---|
| 167 | }
|
---|
| 168 | }
|
---|
| 169 | fprintf(stderr,"\n(%i)Error: Cannot open file: %s\n", P->Par.me, name);
|
---|
[64fa9e] | 170 | Free(name, "OpenFileNo: name");
|
---|
[a0bcf1] | 171 | return(0);
|
---|
| 172 | } else {
|
---|
| 173 | if(verbose) fprintf(stderr,"(%i) File is open: %s\n", P->Par.me, name);
|
---|
[64fa9e] | 174 | Free(name, "OpenFileNo: name");
|
---|
[a0bcf1] | 175 | return(1);
|
---|
| 176 | }
|
---|
| 177 | }
|
---|
| 178 |
|
---|
| 179 | /* die nachfolgende Routine wird von seq und par benutzt!!! */
|
---|
| 180 | /* Oeffnet Datei P->Files.mainname+"...".No.postfix mit what*/
|
---|
| 181 |
|
---|
| 182 | /** Opens a file with FileData::mainname+suffix+"."+No(4 digits)+postfix.
|
---|
| 183 | * Only the path under which the main programme resides is tried.
|
---|
| 184 | * \param *P Problem at hand
|
---|
| 185 | * \param **file file pointer array
|
---|
| 186 | * \param *suffix suffix after mainname
|
---|
| 187 | * \param No the number with up to four digits
|
---|
| 188 | * \param *postfix post-suffix at the end of filename
|
---|
| 189 | * \param *what access parameter for the file: r, w, rw, r+, w+ ...
|
---|
| 190 | * \param verbose 1 - print status, 0 - don't print anything
|
---|
| 191 | * \return 0 - could not open file, 1 - file is open
|
---|
| 192 | */
|
---|
| 193 | int OpenFileNoPost(struct Problem *P, FILE** file, const char* suffix, int No, const char* postfix, const char* what, int verbose)
|
---|
| 194 | {
|
---|
| 195 | char* name; /* Zu erzeugender Dateiname */
|
---|
| 196 | name = (char*)
|
---|
| 197 | Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(postfix) + strlen(suffix) + 6,"OpenFileNoPost");
|
---|
| 198 | sprintf(name, "%s%s%s.%04i%s", P->Files.default_path, P->Files.mainname, suffix, No, postfix);
|
---|
| 199 | *file = fopen(name, what);
|
---|
| 200 | if (*file == 0) {
|
---|
| 201 | fprintf(stderr,"\n(%i)Error: Cannot open file: %s\n", P->Par.me, name);
|
---|
[64fa9e] | 202 | Free(name, "OpenFileNoPost: name");
|
---|
[a0bcf1] | 203 | return(0);
|
---|
| 204 | } else {
|
---|
| 205 | if(verbose) fprintf(stderr,"(%i) File is open: %s\n", P->Par.me, name);
|
---|
[64fa9e] | 206 | Free(name, "OpenFileNoPost: name");
|
---|
[a0bcf1] | 207 | return(1);
|
---|
| 208 | }
|
---|
| 209 | }
|
---|
| 210 |
|
---|
| 211 | /** Opens a file with FileData::mainname+suffix+"."+No1(4 digits)+"."+No2(4 digits).
|
---|
| 212 | * Only the path under which the main programme resides is tried.
|
---|
| 213 | * \param *P Problem at hand
|
---|
| 214 | * \param **file file pointer array
|
---|
| 215 | * \param *suffix suffix after mainname
|
---|
| 216 | * \param No1 first number with up to four digits
|
---|
| 217 | * \param No2 second number with up to four digits
|
---|
| 218 | * \param *what access parameter for the file: r, w, rw, r+, w+ ...
|
---|
| 219 | * \param verbose 1 - print status, 0 - don't print anything
|
---|
| 220 | * \return 0 - could not open file, 1 - file is open
|
---|
| 221 | */
|
---|
| 222 | int OpenFileNoNo(struct Problem *P, FILE** file, const char* suffix, int No1, int No2, const char* what, int verbose)
|
---|
| 223 | { /* zuerst suffix; No1: lfd, No2: procId */
|
---|
| 224 | char *name;
|
---|
| 225 | name = (char*)
|
---|
| 226 | Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(suffix) + 5 + 6 + 1,"OpenFileNoNo");
|
---|
| 227 | sprintf(name,"%s%s%s.%04i.%04i", P->Files.default_path, P->Files.mainname, suffix, No1, No2);
|
---|
| 228 | *file = fopen(name, what);
|
---|
| 229 | if(*file == 0) {
|
---|
| 230 | fprintf(stderr,"\n(%i)Error: Cannot open file: %s\n", P->Par.me, name);
|
---|
[64fa9e] | 231 | Free(name, "OpenFileNoNo: name");
|
---|
[a0bcf1] | 232 | return(0);
|
---|
| 233 | }
|
---|
| 234 | if(verbose) fprintf(stderr,"(%i) File is open: %s\n", P->Par.me, name);
|
---|
[64fa9e] | 235 | Free(name, "OpenFileNoNo: name");
|
---|
[a0bcf1] | 236 | return(1);
|
---|
| 237 | }
|
---|
| 238 |
|
---|
| 239 | /** Transformation of wave function of Riemann level to zeroth before output as Vis.
|
---|
| 240 | * \param *RT RiemannTensor structure, contains RiemannTensor::LatticeLevel
|
---|
| 241 | * \param *source source wave function array
|
---|
| 242 | * \param *dest destination wave function array
|
---|
| 243 | * \param NF dimensional factor (NDIM, generally 3)
|
---|
| 244 | */
|
---|
| 245 | static void OutVisPosRTransformPosNFRto0(const struct RiemannTensor *RT, fftw_real *source, fftw_real *dest, const int NF)
|
---|
| 246 | {
|
---|
| 247 | struct LatticeLevel *Lev = RT->LevR;
|
---|
| 248 | int es = Lev->NUp0[2]*NF;
|
---|
| 249 | unsigned int cpyes = sizeof(fftw_real)*es;
|
---|
| 250 | int nx=Lev->Plan0.plan->local_nx,ny=Lev->N[1],nz=Lev->N[2],Nx=Lev->NUp0[0],Ny=Lev->NUp0[1];
|
---|
| 251 | int lx,ly,z,Lx,Ly;
|
---|
| 252 | for(lx=0; lx < nx; lx++)
|
---|
| 253 | for(Lx=0; Lx < Nx; Lx++)
|
---|
| 254 | for(ly=0; ly < ny; ly++)
|
---|
| 255 | for(Ly=0; Ly < Ny; Ly++)
|
---|
| 256 | for(z=0; z < nz; z++) {
|
---|
| 257 | memcpy( &dest[es*(z+nz*(Ly+Ny*(ly+ny*(Lx+Nx*lx))))],
|
---|
| 258 | &source[es*(Ly+Ny*(Lx+Nx*(z+nz*(ly+ny*lx))))],
|
---|
| 259 | cpyes);
|
---|
| 260 | }
|
---|
| 261 | }
|
---|
| 262 |
|
---|
| 263 |
|
---|
| 264 | /** prints Norm of each wave function to screen.
|
---|
| 265 | * \param *out output destination (eg. stdout)
|
---|
| 266 | * \param *P Problem at hand
|
---|
| 267 | */
|
---|
| 268 | void OutputNorm (FILE *out, struct Problem *P) {
|
---|
| 269 | struct Lattice *Lat = &P->Lat;
|
---|
| 270 | struct Psis *Psi = &Lat->Psi;
|
---|
| 271 | struct LatticeLevel *Lev = P->R.LevS;
|
---|
| 272 | struct OnePsiElement *OnePsi, *LOnePsi;
|
---|
| 273 | int i;
|
---|
| 274 |
|
---|
| 275 | for (i=0;i<Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT;i++) {
|
---|
| 276 | OnePsi =&Psi->AllPsiStatus[i];
|
---|
| 277 | if (OnePsi->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) { // local one
|
---|
| 278 | LOnePsi = &Psi->LocalPsiStatus[OnePsi->MyLocalNo];
|
---|
| 279 | fprintf(out,"(%i) Norm of Psi %i: %e\n", P->Par.me, OnePsi->MyGlobalNo, GramSchGetNorm2(P,Lev,Lev->LPsi->LocalPsi[OnePsi->MyLocalNo]));
|
---|
| 280 | }
|
---|
| 281 | }
|
---|
| 282 | }
|
---|
| 283 |
|
---|
| 284 | /** Output of current Psis state of source level RunStruct::LevS's LevelPsi::LocalPsi to file \ref suffixsrcpsidat.
|
---|
| 285 | * In case of process 0, the doc file is written in a parsable way with minimisation type RunStruct#CurrentMin, level number
|
---|
| 286 | * LatticeLevel#LevelNo, and number of grid nodes LatticeLevel#N.
|
---|
| 287 | * The two (three) different Psis::SpinType's are discerned and in case of SpinUpDown separate data files opened.
|
---|
| 288 | *
|
---|
| 289 | * Then for every global wave function of the desired type each coefficient of the reciprocal grid is copied into
|
---|
| 290 | * Density::DensityCArray[TempDensity], fouriertransformed, copied into Density::DensityArray[TempDensity].
|
---|
| 291 | *
|
---|
| 292 | * As the file output is only handled by process 0, all other coefficient shares are sent within the Psi group to which
|
---|
| 293 | * the current global wave function belongs to process 0 of the Psi group and from there on to process 0 of all via
|
---|
| 294 | * MPI.
|
---|
| 295 | * \param *P Problem at hand
|
---|
| 296 | * \param type Current minimisation state
|
---|
| 297 | * \note This serves as a backup file, when the process is terminated and one would like to restart it at the
|
---|
| 298 | * current calculation lateron, see ReadSrcPsiDensity(). Note also that it is not necessary to specify the
|
---|
| 299 | * same number of processes on later restart, any number may be used under the condition that the number of
|
---|
| 300 | * grid nodes match and that there 2 sharing wave functions in case of SpinUpDown.
|
---|
| 301 | * \sa ReadSrcPsiDensity() - same for Reading the coefficients, TestReadnWriteSrcDensity() - checks both routines against
|
---|
| 302 | * each other
|
---|
| 303 | */
|
---|
| 304 | void OutputSrcPsiDensity(struct Problem *P, enum PsiTypeTag type)
|
---|
| 305 | {
|
---|
| 306 | int i,j,k, Index, zahl, owner;
|
---|
| 307 | struct Lattice *Lat = &P->Lat;
|
---|
| 308 | //struct RunStruct *R = &P->R;
|
---|
| 309 | struct fft_plan_3d *plan = Lat->plan;
|
---|
| 310 | struct LatticeLevel *LevS = P->R.LevS;
|
---|
| 311 | struct Psis *Psi = &Lat->Psi;
|
---|
| 312 | fftw_complex *work = (fftw_complex *)LevS->Dens->DensityArray[TempDensity];
|
---|
| 313 | double *destpsiR = (double *)LevS->Dens->DensityArray[TempDensity];
|
---|
| 314 | fftw_real *srcpsiR = (fftw_real *)LevS->Dens->DensityCArray[TempDensity];
|
---|
| 315 | fftw_complex *srcpsiC = (fftw_complex *)LevS->Dens->DensityCArray[TempDensity];
|
---|
| 316 | FILE *SrcPsiData, *SrcPsiDoc;
|
---|
| 317 | char suffixdat[255], suffixdoc[255];
|
---|
| 318 | MPI_Status status;
|
---|
| 319 | struct OnePsiElement *OnePsiA, *LOnePsiA;
|
---|
| 320 | fftw_complex *LPsiDatA;
|
---|
| 321 | int Num = 0, colorNo = 0;
|
---|
| 322 | int Sent = 0, sent = 0;
|
---|
| 323 |
|
---|
| 324 | SpeedMeasure(P,ReadnWriteTime,StartTimeDo);
|
---|
| 325 | sprintf(suffixdat, ".%.254s.L%i", P->R.MinimisationName[type], LevS->LevelNo);
|
---|
| 326 | strncpy (suffixdoc, suffixdat, 255);
|
---|
| 327 | // for the various spin cases, output the doc-file if it's process 0
|
---|
| 328 | if (!(P->Par.me_comm_ST)) {
|
---|
| 329 | switch (Lat->Psi.PsiST) {
|
---|
| 330 | case SpinDouble:
|
---|
| 331 | colorNo = 0;
|
---|
| 332 | strncat (suffixdat, suffixsrcpsidat, 255-strlen(suffixdat));
|
---|
| 333 | strncat (suffixdoc, suffixsrcpsidoc, 255-strlen(suffixdoc));
|
---|
| 334 | Num = Lat->Psi.GlobalNo[PsiMaxNoDouble];
|
---|
| 335 | break;
|
---|
| 336 | case SpinUp:
|
---|
| 337 | colorNo = 0;
|
---|
| 338 | strncat (suffixdat, suffixsrcpsiupdat, 255-strlen(suffixdat));
|
---|
| 339 | strncat (suffixdoc, suffixsrcpsiupdoc, 255-strlen(suffixdoc));
|
---|
| 340 | Num = Lat->Psi.GlobalNo[PsiMaxNoUp];
|
---|
| 341 | break;
|
---|
| 342 | case SpinDown:
|
---|
| 343 | colorNo = 1;
|
---|
| 344 | strncat (suffixdat, suffixsrcpsidowndat, 255-strlen(suffixdat));
|
---|
| 345 | strncat (suffixdoc, suffixsrcpsidowndoc, 255-strlen(suffixdoc));
|
---|
| 346 | Num = Lat->Psi.GlobalNo[PsiMaxNoDown];
|
---|
| 347 | break;
|
---|
| 348 | }
|
---|
| 349 | if (!OpenFileNo(P, &SrcPsiData, suffixdat, colorNo, "wb",P->Call.out[ReadOut])) // open SourcePsiData as write binary
|
---|
| 350 | fprintf(stderr,"(%i) Error opening file with suffix %s for writing!\n",P->Par.me, suffixdat);
|
---|
| 351 | if (!(P->Par.my_color_comm_ST_Psi)) { // if we are process 0
|
---|
| 352 | if (!OpenFile(P, &SrcPsiDoc, suffixdoc, "w",P->Call.out[ReadOut])) // open the (text) doc file
|
---|
| 353 | fprintf(stderr,"(%i) Error opening file with suffix %s for writing!\n",P->Par.me, suffixdoc);
|
---|
| 354 | fprintf(SrcPsiDoc, "Mintype\t%i\n", (int)type);
|
---|
| 355 | fprintf(SrcPsiDoc, "LevelNo\t%i\n", LevS->LevelNo);
|
---|
| 356 | fprintf(SrcPsiDoc, "GridNodes\t%i\t%i\t%i\n", LevS->N[0], LevS->N[1], LevS->N[2]);
|
---|
| 357 | fprintf(SrcPsiDoc, "PsiNo\t%i\t%i\n", Num, P->Lat.Psi.GlobalNo[PsiMaxAdd]);
|
---|
| 358 | fprintf(SrcPsiDoc, "Epsilon\t%lg\t%lg\n", P->R.RelEpsTotalEnergy, P->R.RelEpsKineticEnergy);
|
---|
| 359 | for (i = 0; i < P->Par.Max_my_color_comm_ST_Psi; i++) {
|
---|
| 360 | fprintf(SrcPsiDoc, "\t%i", Lat->Psi.RealAllLocalNo[i]); // print number of local Psis
|
---|
| 361 | }
|
---|
| 362 | fprintf(SrcPsiDoc, "\n");
|
---|
| 363 | fclose(SrcPsiDoc);
|
---|
| 364 | }
|
---|
| 365 | }
|
---|
| 366 |
|
---|
| 367 | // send/receive around and write share of coefficient array of each wave function
|
---|
| 368 | MPI_Allreduce(&sent, &Sent, 1, MPI_INT, MPI_SUM, P->Par.comm_ST); // catch all at the starter line
|
---|
| 369 | fprintf(stderr,"(%i) me (%i/%i) \t Psi (%i/%i)\t PsiT (%i/%i)\n", P->Par.me, P->Par.me_comm_ST, P->Par.Max_me_comm_ST, P->Par.me_comm_ST_Psi, P->Par.Max_me_comm_ST_Psi, P->Par.me_comm_ST_PsiT, P->Par.Max_me_comm_ST_PsiT);
|
---|
| 370 | k = -1; // k is global PsiNo counter for the desired group
|
---|
| 371 | for (j=0; j < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; j++) { // go through all wave functions (plus the extra one for each process)
|
---|
| 372 | OnePsiA = &Psi->AllPsiStatus[j]; // grab OnePsiA
|
---|
| 373 | if (OnePsiA->PsiType == type) { // only take desired minimisation group
|
---|
| 374 | k++;
|
---|
| 375 | owner = 0; // notes down in process 0 of each psi group the owner of the next coefficient share
|
---|
| 376 | //fprintf(stderr,"(%i) ST_Psi: OnePsiA %i\tP->Par.me %i\n", P->Par.me,OnePsiA->my_color_comm_ST_Psi,P->Par.my_color_comm_ST_Psi);
|
---|
| 377 | if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) { // Belongs to my Psi group?
|
---|
| 378 | LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo];
|
---|
| 379 | LPsiDatA = LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo];
|
---|
| 380 | SetArrayToDouble0((double *)srcpsiR, LevS->Dens->TotalSize*2); // zero DensityCArray[TempDensity]
|
---|
| 381 | for (i=0;i<LevS->MaxG;i++) { // for each every unique G grid vector
|
---|
| 382 | Index = LevS->GArray[i].Index;
|
---|
| 383 | srcpsiC[Index].re = LPsiDatA[i].re; // copy real value
|
---|
| 384 | srcpsiC[Index].im = LPsiDatA[i].im; // copy imaginary value
|
---|
| 385 | }
|
---|
| 386 | for (i=0; i<LevS->MaxDoubleG; i++) { // also for every doubly appearing G vector (symmetry savings)
|
---|
| 387 | srcpsiC[LevS->DoubleG[2*i+1]].re = srcpsiC[LevS->DoubleG[2*i]].re;
|
---|
| 388 | srcpsiC[LevS->DoubleG[2*i+1]].im = -srcpsiC[LevS->DoubleG[2*i]].im;
|
---|
| 389 | }
|
---|
| 390 | // do an fft transform from complex to real on these srcPsiR
|
---|
| 391 | fft_3d_complex_to_real(plan, LevS->LevelNo, FFTNF1, srcpsiC, work);
|
---|
| 392 |
|
---|
| 393 | for (i=0; i < LevS->Dens->LocalSizeR; i++)
|
---|
| 394 | destpsiR[i] = (double)srcpsiR[i];
|
---|
| 395 | } else
|
---|
| 396 | LOnePsiA = NULL;
|
---|
| 397 |
|
---|
| 398 | if (P->Par.me_comm_ST == 0) { // if we are process 0 of all, only we may access the file
|
---|
| 399 | for (i=0; i<P->Par.Max_me_comm_ST_Psi;i++) { // for each share of the coefficient in the PsiGroup
|
---|
| 400 | if (LOnePsiA == NULL) { // if it's not local, receive coefficients from correct PsiGroup (process 0 within that group)
|
---|
| 401 | if (MPI_Recv( destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, OnePsiA->my_color_comm_ST_Psi, OutputSrcPsiTag, P->Par.comm_ST_PsiT, &status ) != MPI_SUCCESS)
|
---|
| 402 | Error(SomeError, "OutputSrcPsiDensity: MPI_Recv of loaded coefficients failed!");
|
---|
| 403 | MPI_Get_count(&status, MPI_DOUBLE, &zahl);
|
---|
| 404 | if (zahl != LevS->Dens->LocalSizeR) // check number of elements
|
---|
| 405 | fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, i, zahl, LevS->Dens->LocalSizeR);
|
---|
| 406 | //else
|
---|
| 407 | //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, i);
|
---|
| 408 | } else { // if it's local ...
|
---|
| 409 | if (i != 0) { // but share of array not for us, receive from owner process within Psi group
|
---|
| 410 | if (MPI_Recv( destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, i, OutputSrcPsiTag, P->Par.comm_ST_Psi, &status ) != MPI_SUCCESS)
|
---|
| 411 | Error(SomeError, "OutputSrcPsiDensity: MPI_Recv of loaded coefficients failed!");
|
---|
| 412 | MPI_Get_count(&status, MPI_DOUBLE, &zahl);
|
---|
| 413 | if (zahl != LevS->Dens->LocalSizeR) // check number of elements
|
---|
| 414 | fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, i, zahl, LevS->Dens->LocalSizeR);
|
---|
| 415 | //else
|
---|
| 416 | //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, i);
|
---|
| 417 | } // otherwise it was our share already
|
---|
| 418 | }
|
---|
| 419 | // store the final share on disc
|
---|
| 420 | if ((zahl = fwrite(destpsiR, sizeof(double), (size_t)(LevS->Dens->LocalSizeR), SrcPsiData)) != (size_t)(LevS->Dens->LocalSizeR)) {
|
---|
| 421 | fclose(SrcPsiData);
|
---|
| 422 | //if (P->Par.me == 0)
|
---|
| 423 | fprintf(stderr, "(%i)OutputSrcPsiDensity: only %i bytes written instead of expected %i\n", P->Par.me, zahl, LevS->Dens->LocalSizeR);
|
---|
| 424 | Error(SomeError,"OutputSrcPsiDensity: fwrite Error");
|
---|
| 425 | }
|
---|
| 426 | }
|
---|
| 427 | } else { // if we are not process 0 of all, we are but a deliverer
|
---|
| 428 | if (LOnePsiA != NULL) { // send if it's local
|
---|
| 429 | if (P->Par.me_comm_ST_Psi == 0) { // if we are process 0 in the group, send final share to our process 0
|
---|
| 430 | for (owner = 0; owner < P->Par.Max_me_comm_ST_Psi; owner++) { // for all processes of our Psi group
|
---|
| 431 | if (owner != 0) { // still not our share of coefficients, receive from owner in our Psi group (increasing owner)
|
---|
| 432 | if (MPI_Recv( destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, owner, OutputSrcPsiTag, P->Par.comm_ST_Psi, &status ) != MPI_SUCCESS)
|
---|
| 433 | Error(SomeError, "OutputSrcPsiDensity: MPI_Recv of loaded coefficients failed!");
|
---|
| 434 | MPI_Get_count(&status, MPI_DOUBLE, &zahl);
|
---|
| 435 | if (zahl != LevS->Dens->LocalSizeR) // check number of elements
|
---|
| 436 | fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, owner, zahl, LevS->Dens->LocalSizeR);
|
---|
| 437 | //else
|
---|
| 438 | //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, owner);
|
---|
| 439 | } else sent++; // only count sent if it was our share
|
---|
| 440 |
|
---|
| 441 | if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, 0, OutputSrcPsiTag, P->Par.comm_ST_PsiT) != MPI_SUCCESS)
|
---|
| 442 | Error(SomeError, "OutputSrcPsiDensity: MPI_Send of loaded coefficients failed!");
|
---|
| 443 | //else
|
---|
| 444 | //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Send to process %i in PsiT group of loaded coefficients GlobalNo %i succeeded!\n", P->Par.me, 0, k);
|
---|
| 445 | }
|
---|
| 446 | } else {
|
---|
| 447 | sent++;
|
---|
| 448 | if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, 0, OutputSrcPsiTag, P->Par.comm_ST_Psi) != MPI_SUCCESS)
|
---|
| 449 | Error(SomeError, "OutputSrcPsiDensity: MPI_Send of loaded coefficients failed!");
|
---|
| 450 | //else
|
---|
| 451 | //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Send to process %i in Psi group of loaded coefficients GlobalNo %i succeeded!\n", P->Par.me, 0, k);
|
---|
| 452 | }
|
---|
| 453 | }
|
---|
| 454 | // otherwise we don't have anything to do with this
|
---|
| 455 | }
|
---|
| 456 | }
|
---|
| 457 | }
|
---|
| 458 | MPI_Allreduce(&sent, &Sent, 1, MPI_INT, MPI_SUM, P->Par.comm_ST); // catch all again at finish
|
---|
| 459 | fprintf(stderr,"(%i) Out of %i shares %i had to be sent in total, %i from this process alone.\n", P->Par.me, P->Par.Max_me_comm_ST_Psi*Psi->NoOfPsis, Sent, sent);
|
---|
| 460 | if (!(P->Par.me_comm_ST))
|
---|
| 461 | fclose(SrcPsiData);
|
---|
| 462 | SpeedMeasure(P,ReadnWriteTime,StopTimeDo);
|
---|
| 463 | }
|
---|
| 464 |
|
---|
| 465 | /** Tests whether writing and successant reading of coefficient array is working correctly.
|
---|
| 466 | * The local wave function array is written to a disc (\sa OutputSrcPsiDensity()), the by \a wavenr
|
---|
| 467 | * specified coefficient array copied to OldPsiDat and afterwards read again via ReadSrcPsiDensity().
|
---|
| 468 | * Comparison per process of each local coefficient shows incorrect read or writes.
|
---|
| 469 | * \param *P Problem at hand
|
---|
| 470 | * \param type minimisation type array to test for read and write
|
---|
| 471 | * \return 1 - successful, 0 - test failed
|
---|
| 472 | */
|
---|
| 473 | int TestReadnWriteSrcDensity(struct Problem *P, enum PsiTypeTag type)
|
---|
| 474 | {
|
---|
| 475 | struct RunStruct *R = &P->R;
|
---|
| 476 | struct LatticeLevel *LevS = R->LevS;
|
---|
| 477 | struct Lattice *Lat = &P->Lat;
|
---|
| 478 | struct Psis *Psi = &Lat->Psi;
|
---|
| 479 | int i,k;
|
---|
| 480 | fftw_complex *destpsiC, *srcpsiC;
|
---|
| 481 |
|
---|
| 482 | //fprintf(stderr,"(%i)TestReadnWriteSrcDensity\n",P->Par.me);
|
---|
| 483 | // write whole array of type to disc
|
---|
| 484 | OutputSrcPsiDensity(P,type);
|
---|
| 485 | //debug(P,"array written");
|
---|
| 486 |
|
---|
| 487 | // copy specified array to OldPsiDat
|
---|
| 488 | for (k=0; k < Lat->Psi.LocalNo; k++) // for every local wave function of type, copy coefficients
|
---|
| 489 | if (Psi->LocalPsiStatus[k].PsiType == type) { // ... yet only for given type
|
---|
| 490 | srcpsiC = LevS->LPsi->LocalPsi[k];
|
---|
| 491 | destpsiC = LevS->LPsi->OldLocalPsi[k - Psi->TypeStartIndex[type]];
|
---|
| 492 | for (i=0;i<LevS->MaxG;i++) { // for each every unique G grid vector
|
---|
| 493 | destpsiC[i].re = srcpsiC[i].re; // copy real value
|
---|
| 494 | destpsiC[i].im = srcpsiC[i].im; // copy imaginary value
|
---|
| 495 | }
|
---|
| 496 | }
|
---|
| 497 | //debug(P,"array copied");
|
---|
| 498 |
|
---|
| 499 | // read whole array again
|
---|
| 500 | if (!ReadSrcPsiDensity(P,type,0,R->LevSNo))
|
---|
| 501 | return 0;
|
---|
| 502 | //debug(P,"array read");
|
---|
| 503 |
|
---|
| 504 | // compare with copied array
|
---|
| 505 | for (k=0; k < Lat->Psi.LocalNo; k++) // for every local wave function of type, compare coefficients
|
---|
| 506 | if (Psi->LocalPsiStatus[k].PsiType == type) { // ... yet only for given type
|
---|
| 507 | srcpsiC = LevS->LPsi->LocalPsi[k];
|
---|
| 508 | destpsiC = LevS->LPsi->OldLocalPsi[k - Psi->TypeStartIndex[type]];
|
---|
| 509 | for (i=0;i<LevS->MaxG;i++) // for each every unique G grid vector
|
---|
| 510 | if ((fabs(destpsiC[i].re - srcpsiC[i].re) >= MYEPSILON) ||(fabs(destpsiC[i].im - srcpsiC[i].im) >= MYEPSILON)) {
|
---|
| 511 | fprintf(stderr,"(%i)TestReadnWriteSrcDensity: First difference at index %i - %lg+i%lg against loaded %lg+i%lg\n",P->Par.me, i, srcpsiC[i].re, srcpsiC[i].im,destpsiC[i].re,destpsiC[i].im);
|
---|
| 512 | return 0;
|
---|
| 513 | }
|
---|
| 514 | }
|
---|
| 515 | //debug(P,"array compared");
|
---|
| 516 | fprintf(stderr,"(%i)TestReadnWriteSrcDensity: OK!\n",P->Par.me);
|
---|
| 517 | return 1;
|
---|
| 518 | }
|
---|
| 519 |
|
---|
| 520 |
|
---|
| 521 | /** Read Psis state from an earlier run.
|
---|
| 522 | * The doc file is opened, mesh sizes LatticeLevel::N[], global number of Psis read and checked against the known
|
---|
| 523 | * values in the inital level RunStruct::LevS.
|
---|
| 524 | * Note, only process 0 handles the files, all read coefficients are transfered to their respective owners via MPI
|
---|
| 525 | * afterwards. Here, process 0 of a certain Psi group is used as a transferer for the coefficients of the other processes
|
---|
| 526 | * in this Psi group. He receives them all from process 0 and sends them onward accordingly. The complete set of
|
---|
| 527 | * coefficients on the real grid for one wave function in the Psi group are transformed into complex wave function
|
---|
| 528 | * coefficients by the usual fft procedure (see ChangeToLevUp()).
|
---|
| 529 | *
|
---|
| 530 | * \param *P Problem at hand
|
---|
| 531 | * \param type minimisation type to read
|
---|
| 532 | * \param test whether to just test for presence of files (1) or really read them (0)
|
---|
| 533 | * \param LevSNo level number to be read
|
---|
| 534 | * \note This is the counterpart to OutputSrcPsiDensity().
|
---|
| 535 | * \return 1 - success, 0 - failure
|
---|
| 536 | * \note It is not necessary to specify the same number of processes on later restart, any number may be used under
|
---|
| 537 | * the condition that the number of grid nodes match and that there at least 2 processes sharing wave functions
|
---|
| 538 | * in case of SpinUpDown.
|
---|
| 539 | * \sa OutputSrcPsiDensity() - same for writing the coefficients, TestReadnWriteSrcDensity() - checks both routines against
|
---|
| 540 | * each other
|
---|
| 541 | */
|
---|
| 542 | int ReadSrcPsiDensity(struct Problem *P, enum PsiTypeTag type, int test, int LevSNo)
|
---|
| 543 | {
|
---|
| 544 | int i, j, k, Index, owner;
|
---|
| 545 | struct RunStruct *R = &P->R;
|
---|
| 546 | struct Lattice *Lat = &P->Lat;
|
---|
| 547 | struct fft_plan_3d *plan = Lat->plan;
|
---|
| 548 | struct LatticeLevel *LevS = R->LevS; // keep open for LevelNo read from file
|
---|
| 549 | struct Psis *Psi = &Lat->Psi;
|
---|
| 550 | //struct Energy *E = Lat->E;
|
---|
| 551 | fftw_complex *work;
|
---|
| 552 | double *destpsiR;
|
---|
| 553 | fftw_real *srcpsiR;
|
---|
| 554 | fftw_complex *srcpsiC;
|
---|
| 555 | FILE *SrcPsiData, *SrcPsiDoc;
|
---|
| 556 | int N[NDIM], GlobalNo[2];
|
---|
| 557 | int LevelNo, readnr=0;
|
---|
| 558 | int zahl, signal = test ? 1 : 2; // 0 - ok, 1 - test failed, 2 - throw Error
|
---|
| 559 | char suffixdat[255], suffixdoc[255];
|
---|
| 560 | int read_type, Num = 0, colorNo = 0;
|
---|
| 561 | char spin[20];
|
---|
| 562 | double Eps[2];
|
---|
| 563 | MPI_Status status;
|
---|
| 564 | struct OnePsiElement *OnePsiA, *LOnePsiA;
|
---|
| 565 | int Recv=0, recv=0;
|
---|
| 566 |
|
---|
| 567 | SpeedMeasure(P,ReadnWriteTime,StartTimeDo);
|
---|
| 568 | sprintf(suffixdat, ".%.254s.L%i", P->R.MinimisationName[type], LevSNo);
|
---|
| 569 | strncpy (suffixdoc, suffixdat, 255);
|
---|
| 570 | // Depending on Psis::SpinType the source psi doc file is opened and header written
|
---|
| 571 | switch (Lat->Psi.PsiST) {
|
---|
| 572 | case SpinDouble:
|
---|
| 573 | colorNo = 0;
|
---|
| 574 | strncat (suffixdat, suffixsrcpsidat, 255-strlen(suffixdat));
|
---|
| 575 | strncat (suffixdoc, suffixsrcpsidoc, 255-strlen(suffixdoc));
|
---|
| 576 | strncpy (spin, "GlobalNoSpinDouble", 20);
|
---|
| 577 | Num = Lat->Psi.GlobalNo[PsiMaxNoDouble];
|
---|
| 578 | break;
|
---|
| 579 | case SpinUp:
|
---|
| 580 | colorNo = 0;
|
---|
| 581 | strncat (suffixdat, suffixsrcpsiupdat, 255-strlen(suffixdat));
|
---|
| 582 | strncat (suffixdoc, suffixsrcpsiupdoc, 255-strlen(suffixdoc));
|
---|
| 583 | strncpy (spin, "GlobalNoSpinUp", 20);
|
---|
| 584 | Num = Lat->Psi.GlobalNo[PsiMaxNoUp];
|
---|
| 585 | break;
|
---|
| 586 | case SpinDown:
|
---|
| 587 | colorNo = 1;
|
---|
| 588 | strncat (suffixdat, suffixsrcpsidowndat, 255-strlen(suffixdat));
|
---|
| 589 | strncat (suffixdoc, suffixsrcpsidowndoc, 255-strlen(suffixdoc));
|
---|
| 590 | strncpy (spin, "GlobalNoSpinDown", 20);
|
---|
| 591 | Num = Lat->Psi.GlobalNo[PsiMaxNoDown];
|
---|
| 592 | break;
|
---|
| 593 | }
|
---|
| 594 | // open doc file ...
|
---|
| 595 | if (!(P->Par.me_comm_ST)) {
|
---|
| 596 | if (!OpenFile(P, &SrcPsiDoc, suffixdoc, "r", test ? 0 : P->Call.out[ReadOut])) { // open doc file
|
---|
| 597 | debug(P,"ReadSrcPsiDensity: doc file pointer NULL\n");
|
---|
| 598 | if (test) {
|
---|
| 599 | if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
|
---|
| 600 | Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
|
---|
| 601 | return 0;
|
---|
| 602 | }
|
---|
| 603 | if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
|
---|
| 604 | Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
|
---|
| 605 | Error(SomeError,"ReadSrcPsiDensity: cannot open doc file!");
|
---|
| 606 | }
|
---|
| 607 | // ... and parse critical ...
|
---|
| 608 | readnr += ParseForParameter(0,SrcPsiDoc,"Mintype",0,1,1,int_type,(int *)&read_type, test ? optional : critical);
|
---|
| 609 | readnr += ParseForParameter(0,SrcPsiDoc,"LevelNo",0,1,1,int_type,&LevelNo, test ? optional : critical);
|
---|
| 610 | readnr += 3*ParseForParameter(0,SrcPsiDoc,"GridNodes",0,3,1,row_int,&N[0], test ? optional : critical);
|
---|
| 611 | readnr += 2*ParseForParameter(0,SrcPsiDoc,"PsiNo",0,2,1,row_int,&GlobalNo[0], test ? optional : critical);
|
---|
| 612 | // and optional items ...
|
---|
| 613 | if (ParseForParameter(0,SrcPsiDoc,"Epsilon",0,2,1,row_double,&Eps[0],optional))
|
---|
| 614 | if ((P->Call.ReadSrcFiles == 1) && ((Eps[1] < R->RelEpsKineticEnergy) || (Eps[0] < R->RelEpsTotalEnergy))) {
|
---|
| 615 | //fprintf(stderr,"(%i) Eps %lg %lg\tRelEps %lg %lg\n", P->Par.me, Eps[0], Eps[1], R->RelEpsTotalEnergy, R->RelEpsKineticEnergy);
|
---|
| 616 | fprintf(stderr,"(%i) NOTE: Doing minimization after source file parsing due to smaller specified epsilon stop conditions.\n",P->Par.me);
|
---|
| 617 | P->Call.ReadSrcFiles = 2; // do minimisation even after loading
|
---|
| 618 | }
|
---|
| 619 | if (readnr != 7) { // check number of items read
|
---|
| 620 | debug(P, "ReadSrcPsiDensity: too few doc items in file\n");
|
---|
| 621 | fclose(SrcPsiDoc);
|
---|
| 622 | if (test) {
|
---|
| 623 | if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
|
---|
| 624 | Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
|
---|
| 625 | return 0;
|
---|
| 626 | }
|
---|
| 627 | fprintf(stderr,"ReadSrcPsiDensity: Only %i items read!\n",readnr);
|
---|
| 628 | if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
|
---|
| 629 | Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
|
---|
| 630 | Error(SomeError, "ReadSrcPsiDensity: read error");
|
---|
| 631 | }
|
---|
| 632 | // check if levels match
|
---|
| 633 | if (LevSNo != LevelNo) {
|
---|
| 634 | debug(P,"ReadSrcPsiDensity: mismatching levels\n");
|
---|
| 635 | fclose(SrcPsiDoc);
|
---|
| 636 | if (test) {
|
---|
| 637 | if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
|
---|
| 638 | Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
|
---|
| 639 | return 0;
|
---|
| 640 | }
|
---|
| 641 | if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
|
---|
| 642 | Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
|
---|
| 643 | Error(SomeError,"ReadSrcPsiDensity: Mismatching levels!");
|
---|
| 644 | } else {
|
---|
| 645 | LevS = &P->Lat.Lev[LevelNo];
|
---|
| 646 | }
|
---|
| 647 |
|
---|
| 648 | // check if systems in memory and file match
|
---|
| 649 | if ((read_type != R->CurrentMin) || (N[0] != LevS->N[0]) || (N[1] != LevS->N[1]) || (N[2] != LevS->N[2]) || (GlobalNo[0] != Num) || (GlobalNo[1] != Lat->Psi.GlobalNo[PsiMaxAdd])) { // check read system
|
---|
| 650 | debug(P,"ReadSrcPsiDensity: srcpsi file does not fit to system\n");
|
---|
| 651 | fclose(SrcPsiDoc);
|
---|
| 652 | if (test) {
|
---|
| 653 | fprintf(stderr,"(%i) Min\t N(x,y,z)\tPsiNo+AddNo\n file: %s\t %i %i %i\t %i + %i\nsystem: %s\t %d %d %d\t %d + %d\n",P->Par.me, R->MinimisationName[read_type], N[0], N[1], N[2], GlobalNo[0], GlobalNo[1], R->MinimisationName[R->CurrentMin], LevS->N[0] , LevS->N[1], LevS->N[2], Lat->Psi.GlobalNo[PsiMaxNoDouble], Lat->Psi.GlobalNo[PsiMaxAdd]);
|
---|
| 654 | if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
|
---|
| 655 | Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
|
---|
| 656 | return 0;
|
---|
| 657 | }
|
---|
| 658 | fprintf(stderr,"ReadSrcPsiDensity: Type %i != CurrentMin %i || N[0] %i != %i || N[1] %i != %i || N[2] %i != %i || %s %i + %i != %i + %i\n", read_type, R->CurrentMin, N[0], LevS->N[0], N[1], LevS->N[1], N[2], LevS->N[2], spin, GlobalNo[0], GlobalNo[1], Num, P->Lat.Psi.GlobalNo[PsiMaxAdd]);
|
---|
| 659 | if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
|
---|
| 660 | Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
|
---|
| 661 | Error(SomeError,"ReadSrcPsiDensity: srcpsi file does not fit to system");
|
---|
| 662 | }
|
---|
| 663 | signal = 0; // everything went alright, signal ok
|
---|
| 664 | if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
|
---|
| 665 | Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
|
---|
| 666 | } else { // others wait for signal from root process
|
---|
| 667 | if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
|
---|
| 668 | Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
|
---|
| 669 | if (signal == 1)
|
---|
| 670 | return 0;
|
---|
| 671 | else if (signal == 2)
|
---|
| 672 | Error(SomeError, "ReadSrcPsiDensity: Something went utterly wrong, see root process");
|
---|
| 673 | else
|
---|
| 674 | fprintf(stderr,"(%i) ReadSrcPsiDensity: Everything went alright so far\n", P->Par.me);
|
---|
| 675 | }
|
---|
| 676 | if (MPI_Bcast(&LevelNo,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
|
---|
| 677 | Error(SomeError,"ReadSrcPsiDensity: Bcast of LevelNo failed\n");
|
---|
| 678 | LevS = &P->Lat.Lev[LevelNo];
|
---|
| 679 | //if (!test) fprintf(stderr,"(%i) LevelSNo %i\n", P->Par.me, LevS->LevelNo);
|
---|
| 680 |
|
---|
| 681 | if (!test) {
|
---|
| 682 | // set some pointers for work to follow
|
---|
| 683 | work = (fftw_complex *)LevS->Dens->DensityArray[TempDensity];
|
---|
| 684 | destpsiR = (double *)LevS->Dens->DensityArray[TempDensity];
|
---|
| 685 | srcpsiR = (fftw_real *)LevS->Dens->DensityCArray[TempDensity];
|
---|
| 686 | srcpsiC = (fftw_complex *)LevS->Dens->DensityCArray[TempDensity];
|
---|
| 687 |
|
---|
| 688 | // read share of coefficient array for each wave function and send/receive around
|
---|
| 689 | owner = 0;
|
---|
| 690 | MPI_Allreduce (&recv, &Recv, 1, MPI_INT, MPI_SUM, P->Par.comm_ST);
|
---|
| 691 | //fprintf(stderr,"(%i) me (%i/%i) \t Psi (%i/%i)\t PsiT (%i/%i)\n", P->Par.me, P->Par.me_comm_ST, P->Par.Max_me_comm_ST, P->Par.me_comm_ST_Psi, P->Par.Max_me_comm_ST_Psi, P->Par.me_comm_ST_PsiT, P->Par.Max_me_comm_ST_PsiT);
|
---|
| 692 | k = -1; // k is global PsiNo counter for the desired group
|
---|
| 693 | for (j=0; j < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; j++) { // go through all wave functions (plus the extra one for each process)
|
---|
| 694 | OnePsiA = &Psi->AllPsiStatus[j]; // grab OnePsiA
|
---|
| 695 | if (OnePsiA->PsiType == type) { // only take desired minimisation group
|
---|
| 696 | k++;
|
---|
| 697 | //fprintf(stderr,"(%i) ST_Psi: OnePsiA %i\tP->Par.me %i\n", P->Par.me,OnePsiA->my_color_comm_ST_Psi,P->Par.my_color_comm_ST_Psi);
|
---|
| 698 | if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // Belongs to my Psi group?
|
---|
| 699 | LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo];
|
---|
| 700 | else
|
---|
| 701 | LOnePsiA = NULL;
|
---|
| 702 |
|
---|
| 703 | if (P->Par.me_comm_ST == 0) { // if we are process 0 of all, we may access file
|
---|
| 704 | if (!OpenFileNo(P, &SrcPsiData, suffixdat, colorNo, "r", test ? 0 : P->Call.out[ReadOut])) {
|
---|
| 705 | Error(SomeError,"ReadSrcPsiDensity: cannot open data file");
|
---|
| 706 | }
|
---|
| 707 | for (i=P->Par.Max_me_comm_ST_Psi-1; i>=0;i--) { // load coefficients
|
---|
| 708 | fseek( SrcPsiData, (long)((k*N[0]*N[1]*N[2]+i*((long)LevS->Dens->LocalSizeR))*sizeof(double)), SEEK_SET); // seek to beginning of this process' coefficients
|
---|
| 709 | // readin
|
---|
| 710 | if ((zahl = fread(destpsiR, sizeof(double), (size_t)(LevS->Dens->LocalSizeR), SrcPsiData)) != (size_t)LevS->Dens->LocalSizeR) {
|
---|
| 711 | fclose(SrcPsiData);
|
---|
| 712 | fprintf(stderr, "(%i)ReadSrcPsiDensity: only %i bytes read instead of expected %i\n", P->Par.me, zahl, LevS->Dens->LocalSizeR);
|
---|
| 713 | Error(SomeError,"ReadSrcPsiDensity: fread Error");
|
---|
| 714 | }
|
---|
| 715 | if (LOnePsiA == NULL) { // if it's not local, send away coefficients to correct PsiGroup (process 0 within that group)
|
---|
| 716 | if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, OnePsiA->my_color_comm_ST_Psi, ReadSrcPsiTag, P->Par.comm_ST_PsiT) != MPI_SUCCESS)
|
---|
| 717 | Error(SomeError, "ReadSrcPsiDensity: MPI_Send of loaded coefficients failed!");
|
---|
| 718 | //else
|
---|
| 719 | //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Send to process %i of loaded coefficients GlobalNo %i, owner %i succeeded!\n", P->Par.me, OnePsiA->my_color_comm_ST_Psi, k, i);
|
---|
| 720 | } else { // if it's local ...
|
---|
| 721 | if (i != 0) { // but share of array not for us, send to owner process within Psi group
|
---|
| 722 | if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, i, ReadSrcPsiTag, P->Par.comm_ST_Psi) != MPI_SUCCESS)
|
---|
| 723 | Error(SomeError, "ReadSrcPsiDensity: MPI_Send within Psi group of loaded coefficients failed!");
|
---|
| 724 | //else
|
---|
| 725 | //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Send to process %i within Psi group of loaded coefficients GlobalNo %i succeeded!\n", P->Par.me, i, k);
|
---|
| 726 | } // otherwise it was our share already
|
---|
| 727 | }
|
---|
| 728 | }
|
---|
| 729 | } else {
|
---|
| 730 | if (LOnePsiA != NULL) { // receive
|
---|
| 731 | if (P->Par.me_comm_ST_Psi == 0) { // if we are process 0 in the group, receive share from process 0 of all
|
---|
| 732 | for (owner = P->Par.Max_me_comm_ST_Psi -1; owner >=0; owner--) { // for all processes of our Psi group
|
---|
| 733 | if (MPI_Recv(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, 0, ReadSrcPsiTag, P->Par.comm_ST_PsiT, &status) != MPI_SUCCESS)
|
---|
| 734 | Error(SomeError, "ReadSrcPsiDensity: MPI_Recv of loaded coefficients failed!");
|
---|
| 735 | MPI_Get_count(&status, MPI_DOUBLE, &zahl);
|
---|
| 736 | if (zahl != LevS->Dens->LocalSizeR) // check number of elements
|
---|
| 737 | fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Recv from process 0 of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, owner, zahl, LevS->Dens->LocalSizeR);
|
---|
| 738 | //else
|
---|
| 739 | //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Recv from process 0 of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, owner);
|
---|
| 740 |
|
---|
| 741 | if (owner != 0) { // not our share of coefficients, send to owner in our Psi group
|
---|
| 742 | if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, owner, ReadSrcPsiTag, P->Par.comm_ST_Psi) != MPI_SUCCESS)
|
---|
| 743 | Error(SomeError, "ReadSrcPsiDensity: MPI_Send within Psi group of loaded coefficients failed!");
|
---|
| 744 | //else
|
---|
| 745 | //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Send to process %i within Psi group of loaded coefficients GlobalNo %i succeeded!\n", P->Par.me, owner, k);
|
---|
| 746 | } else recv++;
|
---|
| 747 | }
|
---|
| 748 | // otherwise it's our share!
|
---|
| 749 | } else { // our share within Psi Group not belonging to process 0 of all
|
---|
| 750 | recv++;
|
---|
| 751 | if (MPI_Recv(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, 0, ReadSrcPsiTag, P->Par.comm_ST_Psi, &status) != MPI_SUCCESS)
|
---|
| 752 | Error(SomeError, "ReadSrcPsiDensity: MPI_Recv of loaded coefficients failed!");
|
---|
| 753 | MPI_Get_count(&status, MPI_DOUBLE, &zahl);
|
---|
| 754 | if (zahl != LevS->Dens->LocalSizeR) // check number of elements
|
---|
| 755 | fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, P->Par.me_comm_ST_Psi, zahl, LevS->Dens->LocalSizeR);
|
---|
| 756 | //else
|
---|
| 757 | //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, P->Par.me_comm_ST_Psi);
|
---|
| 758 | }
|
---|
| 759 | }
|
---|
| 760 | // otherwise we don't have anything to do with this
|
---|
| 761 | }
|
---|
| 762 |
|
---|
| 763 | if (LOnePsiA != NULL) {
|
---|
| 764 | SetArrayToDouble0((double *)srcpsiR, LevS->Dens->TotalSize*2);
|
---|
| 765 | for (i=0; i < LevS->Dens->LocalSizeR; i++) // copy dest to src
|
---|
| 766 | srcpsiR[i] = (fftw_real)destpsiR[i];
|
---|
| 767 |
|
---|
| 768 | fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, srcpsiC, work); // fft transform
|
---|
| 769 | //if (P->Call.out[PsiOut])
|
---|
| 770 | //fprintf(stderr,"(%i) LevSNo %i\t LocalPsi %p\n", P->Par.me, LevS->LevelNo, LevS->LPsi->LocalPsi[LOnePsiA->MyLocalNo]);
|
---|
| 771 | for (i=0;i<LevS->MaxG;i++) { // and copy into wave functions coefficients
|
---|
| 772 | Index = LevS->GArray[i].Index;
|
---|
| 773 | LevS->LPsi->LocalPsi[LOnePsiA->MyLocalNo][i].re = srcpsiC[Index].re/LevS->MaxN;
|
---|
| 774 | LevS->LPsi->LocalPsi[LOnePsiA->MyLocalNo][i].im = srcpsiC[Index].im/LevS->MaxN;
|
---|
| 775 | }
|
---|
| 776 | }
|
---|
| 777 | if ((P->Par.me_comm_ST == 0) && (SrcPsiData != NULL)) fclose(SrcPsiData);
|
---|
| 778 | }
|
---|
| 779 | }
|
---|
| 780 | MPI_Allreduce (&recv, &Recv, 1, MPI_INT, MPI_SUM, P->Par.comm_ST);
|
---|
| 781 | fprintf(stderr,"(%i) Out of %i shares %i had to be received in total, %i from this process alone.\n", P->Par.me, P->Par.Max_me_comm_ST_Psi*Psi->NoOfPsis, Recv, recv);
|
---|
| 782 | SpeedMeasure(P,ReadnWriteTime,StopTimeDo);
|
---|
| 783 | }
|
---|
| 784 | return 1; // everything went well till the end
|
---|
| 785 | }
|
---|
| 786 |
|
---|
| 787 | /** Creates the density \ref suffixdensdoc and \ref suffixdensdx files for OpenDx.
|
---|
| 788 | * Opens \ref suffixdensdoc, fills (pos&data file name, byte order, max mesh points, matrix alignment, steps)
|
---|
| 789 | * and closes it.
|
---|
| 790 | * Opens \ref suffixdensdx, then for every FileData::OutVisStep a describing structure for DX is written and
|
---|
| 791 | * the file closed again.
|
---|
| 792 | * \param *P Problem at hand
|
---|
| 793 | * \param me my process number in communicator Psi (0 - do nothing, else - do)
|
---|
| 794 | */
|
---|
| 795 | static void CreateDensityOutputGeneral(struct Problem *P, const int me)
|
---|
| 796 | {
|
---|
| 797 | FILE *DensityDoc, *DensityDx;
|
---|
| 798 | char *posname, *datname;
|
---|
| 799 | struct LatticeLevel *Lev = &P->Lat.Lev[STANDARTLEVEL];
|
---|
| 800 | unsigned int i, MaxPoints, N[NDIM];
|
---|
| 801 | double *RB = P->Lat.RealBasis;
|
---|
| 802 | if (me) return;
|
---|
| 803 | N[0] = Lev->N[0]*Lev->NUp[0];
|
---|
| 804 | N[1] = Lev->N[1]*Lev->NUp[1];
|
---|
| 805 | N[2] = Lev->N[2]*Lev->NUp[2];
|
---|
| 806 | MaxPoints = (N[0]+1)*(N[1]+1)*(N[2]+1);
|
---|
| 807 | posname = (char*)
|
---|
| 808 | Malloc(strlen(P->Files.mainname) + strlen(suffixdenspos) + 1,"OpenFile");
|
---|
| 809 | sprintf(posname, "%s%s", P->Files.mainname, suffixdenspos);
|
---|
| 810 | datname = (char*)
|
---|
| 811 | Malloc(strlen(P->Files.mainname) + strlen(suffixdensdat) + 1,"OpenFile");
|
---|
| 812 | sprintf(datname, "%s%s", P->Files.mainname, suffixdensdat);
|
---|
| 813 | // write doc file
|
---|
| 814 | OpenFile(P, &DensityDoc, suffixdensdoc, "w",P->Call.out[ReadOut]);
|
---|
| 815 | fprintf(DensityDoc,"DensityPositions file = %s.####\n", posname);
|
---|
| 816 | fprintf(DensityDoc,"DensityData file = %s.####\n", datname);
|
---|
| 817 | fprintf(DensityDoc,"format = ieee float (Bytes %lu) %s\n",(unsigned long) sizeof(float),msb);
|
---|
| 818 | fprintf(DensityDoc,"points = %i\n", MaxPoints);
|
---|
| 819 | fprintf(DensityDoc,"majority = row\n");
|
---|
| 820 | fprintf(DensityDoc,"TimeSeries = %i\n",P->Files.OutVisStep+1);
|
---|
| 821 | fclose(DensityDoc);
|
---|
| 822 | // write DX file
|
---|
| 823 | OpenFile(P, &DensityDx, suffixdensdx, "w",P->Call.out[ReadOut]);
|
---|
| 824 | for (i=0; i < (unsigned int)P->Files.OutVisStep+1; i++) { // for every OutVis step
|
---|
| 825 | if (i==0) {
|
---|
| 826 | fprintf(DensityDx,"object \"gridcon\" class gridconnections counts %i %i %i\n\n",(N[0]+1),(N[1]+1),(N[2]+1));
|
---|
| 827 | if (P->Files.OutputPosType[i] != active)
|
---|
| 828 | fprintf(DensityDx, "object \"posdens\" class gridpositions counts %i %i %i\norigin 0.0 0.0 0.0\ndelta %f %f %f\ndelta %f %f %f\ndelta %f %f %f\n\n",
|
---|
| 829 | (N[0]+1),(N[1]+1),(N[2]+1),
|
---|
| 830 | (float)(RB[0]/N[0]),(float)(RB[1]/N[0]),(float)RB[2]/N[0],
|
---|
| 831 | (float)(RB[3]/N[1]),(float)(RB[4]/N[1]),(float)RB[5]/N[1],
|
---|
| 832 | (float)(RB[6]/N[2]),(float)(RB[7]/N[2]),(float)RB[8]/N[2]);
|
---|
| 833 | }
|
---|
| 834 | if (P->Files.OutputPosType[i] == active) {
|
---|
| 835 | fprintf(DensityDx,
|
---|
| 836 | "object \"pos.%04u\" class array type float rank 1 shape 3 items %i %s binary\n",i,MaxPoints,msb);
|
---|
| 837 | fprintf(DensityDx,"data file %s.%04u,0\n",posname,i);
|
---|
| 838 | }
|
---|
| 839 | fprintf(DensityDx,"attribute \"dep\" string \"positions\"\n");
|
---|
| 840 | fprintf(DensityDx,"# %lu - %lu Bytes\n\n",MaxPoints*i*(unsigned long)sizeof(float)*NDIM,MaxPoints*(i+1)*(unsigned long)sizeof(float)*NDIM-1);
|
---|
| 841 |
|
---|
| 842 | fprintf(DensityDx,"object \"dat.%04u\" class array type float rank 0 items %i %s binary\n",i,MaxPoints,msb);
|
---|
| 843 | fprintf(DensityDx,"data file %s.%04u,0\n",datname,i);
|
---|
| 844 | fprintf(DensityDx,"attribute \"dep\" string \"positions\"\n");
|
---|
| 845 | fprintf(DensityDx,"# %lu - %lu Bytes\n\n",MaxPoints*i*(unsigned long)sizeof(float),MaxPoints*(i+1)*(unsigned long)sizeof(float)-1);
|
---|
| 846 |
|
---|
| 847 | fprintf(DensityDx,"object \"obj.%04u\" class field\n",i);
|
---|
| 848 | if (P->Files.OutputPosType[i] == active)
|
---|
| 849 | fprintf(DensityDx,"component \"positions\" \"pos.%04i\"\n",i);
|
---|
| 850 | if (P->Files.OutputPosType[i] != active)
|
---|
| 851 | fprintf(DensityDx,"component \"positions\" \"posdens\"\n");
|
---|
| 852 | fprintf(DensityDx,"component \"connections\" \"gridcon\"\n");
|
---|
| 853 | fprintf(DensityDx,"component \"data\" \"dat.%04i\"\n",i);
|
---|
| 854 | }
|
---|
| 855 | fprintf(DensityDx,"\nobject \"series\" class series\n");
|
---|
| 856 | for (i=0; i < (unsigned int)P->Files.OutVisStep+1; i++)
|
---|
| 857 | fprintf(DensityDx,"member %i \"obj.%04u\" position %f\n",i,i,(float)i);
|
---|
| 858 | fprintf(DensityDx,"end\n");
|
---|
| 859 |
|
---|
| 860 | fclose(DensityDx);
|
---|
[64fa9e] | 861 | Free(posname, "CreateDensityOutputGeneral: posname");
|
---|
| 862 | Free(datname, "CreateDensityOutputGeneral: datname");
|
---|
[a0bcf1] | 863 | }
|
---|
| 864 |
|
---|
| 865 | /** Calculates the OutVis density of the RiemannTensor level.
|
---|
| 866 | * The usual pattern arises when a density is fftransformed:
|
---|
| 867 | * -# over all grid vectors up to MaxG
|
---|
| 868 | * -# over all doubly grid vectors up to MaxDoubleG
|
---|
| 869 | * -# call to fft_3d_complex_to_real()
|
---|
| 870 | *
|
---|
| 871 | * In this case here followed by call to OutVisPosRTransformPosNFRto0() and finally FileData::work
|
---|
| 872 | * is copied to FileData::PosR.
|
---|
| 873 | * \param *Lat Lattice structure, containing Lattice::plan and LatticeLevel
|
---|
| 874 | * \param *F FileData structure, containing FileData::PosC, FileData::PosR, FileData::work, FileData::Totalsize, FileData::LocalSizeR
|
---|
| 875 | */
|
---|
| 876 | static void CalculateOutVisDensityPos(struct Lattice *Lat, struct FileData *F/*, const double FactorC_R, const double FactorR_C*/)
|
---|
| 877 | {
|
---|
| 878 | struct fft_plan_3d *plan = Lat->plan;
|
---|
| 879 | struct RiemannTensor *RT = &Lat->RT;
|
---|
| 880 | struct LatticeLevel *LevR = RT->LevR;
|
---|
| 881 | fftw_complex *destC = F->PosC;
|
---|
| 882 | fftw_real *destR = F->PosR;
|
---|
| 883 | fftw_complex *work = F->work;
|
---|
| 884 | fftw_real *workR = (fftw_real*)work;
|
---|
| 885 | fftw_complex *PosFactor = F->PosFactor;
|
---|
| 886 | fftw_complex *posfac, *destpos, *destRCS, *destRCD;
|
---|
| 887 | fftw_complex *coeff;
|
---|
| 888 | fftw_complex source;
|
---|
| 889 | int i, Index, pos, n;
|
---|
| 890 | int NF = NDIM, MaxNUp = F->MaxNUp, TotalSize = F->TotalSize, LocalSizeR = F->LocalSizeR;
|
---|
| 891 | SetArrayToDouble0((double *)destC, TotalSize*2);
|
---|
| 892 | for (i=0;i < LevR->MaxG;i++) {
|
---|
| 893 | Index = LevR->GArray[i].Index;
|
---|
| 894 | posfac = &PosFactor[MaxNUp*i];
|
---|
| 895 | destpos = &destC[MaxNUp*Index*NF];
|
---|
| 896 | coeff = &RT->Coeff[i];
|
---|
| 897 | for (pos=0; pos < MaxNUp; pos++) {
|
---|
| 898 | for (n=0; n < NF; n++) {
|
---|
| 899 | source.re = coeff[n].re;
|
---|
| 900 | source.im = coeff[n].im;
|
---|
| 901 | destpos[n+NF*pos].re = source.re*posfac[pos].re-source.im*posfac[pos].im;
|
---|
| 902 | destpos[n+NF*pos].im = source.re*posfac[pos].im+source.im*posfac[pos].re;
|
---|
| 903 | }
|
---|
| 904 | }
|
---|
| 905 | }
|
---|
| 906 | for (i=0; i < LevR->MaxDoubleG; i++) {
|
---|
| 907 | destRCS = &destC[LevR->DoubleG[2*i]*MaxNUp*NF];
|
---|
| 908 | destRCD = &destC[LevR->DoubleG[2*i+1]*MaxNUp*NF];
|
---|
| 909 | for (pos=0; pos < MaxNUp; pos++) {
|
---|
| 910 | for (n=0; n < NF; n++) {
|
---|
| 911 | destRCD[n+NF*pos].re = destRCS[n+NF*pos].re;
|
---|
| 912 | destRCD[n+NF*pos].im = -destRCS[n+NF*pos].im;
|
---|
| 913 | }
|
---|
| 914 | }
|
---|
| 915 | }
|
---|
| 916 | fft_3d_complex_to_real(plan, LevR->LevelNo, FFTNFRVecUp0, destC, work);
|
---|
| 917 | OutVisPosRTransformPosNFRto0(RT, destR, workR, NF);
|
---|
| 918 | memcpy(destR,workR,sizeof(fftw_real)*LocalSizeR);
|
---|
| 919 | }
|
---|
| 920 |
|
---|
| 921 | /** Prepare Density::DensityArray for output.
|
---|
| 922 | * Into FileData::work subsequently each node (all z, all y, all x) is written as \f$\log(1+x)\f$,
|
---|
| 923 | * where x is Density::DensityArray[TotalDensity]. In the end result is send to process 0 (yet not
|
---|
| 924 | * received here, see CombineOutVisArray()). In case of RiemannTensor use, some more complex calculations
|
---|
| 925 | * are made: FileData::PosR is used, the coefficient offset'ed to the current node and the log taken there.
|
---|
| 926 | * \param *P Problem at hand
|
---|
| 927 | * \param myPE this ranks process in the Psi communcator ParallelSimulationData::me_comm_ST_Psi
|
---|
| 928 | * \param *srcdens Pointer to DensityArray which is to be displayed
|
---|
| 929 | */
|
---|
| 930 | static void OutputOutVisDensity(struct Problem *P, const int myPE, fftw_real *srcdens)
|
---|
| 931 | {
|
---|
| 932 | int N[NDIM], n[NDIM], pos[NDIM];
|
---|
| 933 | int destpos = 0;
|
---|
| 934 | double fac[NDIM], posd[NDIM];
|
---|
| 935 | float posf[NDIM+1];
|
---|
| 936 | struct Lattice *Lat = &P->Lat;
|
---|
| 937 | struct LatticeLevel *Lev0 = &Lat->Lev[0];
|
---|
| 938 | fftw_real *srcpos = P->Files.PosR;
|
---|
| 939 | //fftw_real *srcdens = Lev0->Dens->DensityArray[ActualDensity]; //[TotalDensity]; trick to display single density
|
---|
| 940 | float *dest = (float *)P->Files.work;
|
---|
| 941 | int Nx = Lev0->Plan0.plan->N[0];
|
---|
| 942 | int i;
|
---|
| 943 | double min, max;
|
---|
| 944 |
|
---|
| 945 | N[0] = Lev0->Plan0.plan->local_nx;
|
---|
| 946 | N[1] = Lev0->Plan0.plan->N[1];
|
---|
| 947 | N[2] = Lev0->Plan0.plan->N[2];
|
---|
| 948 |
|
---|
| 949 | max = min = srcdens[0];
|
---|
| 950 | for (i=1;i<P->R.Lev0->Dens->LocalSizeR;i++) {
|
---|
| 951 | if (srcdens[i] < min) min = srcdens[i];
|
---|
| 952 | if (srcdens[i] > max) max = srcdens[i];
|
---|
| 953 | }
|
---|
| 954 | if (P->Call.out[PsiOut]) fprintf(stderr,"(%i)OutputOutVisDensity: min %e\tmax %e\n",P->Par.me, min, max);
|
---|
| 955 |
|
---|
| 956 | // go through all nodes
|
---|
| 957 | for (n[0]=0; n[0] < N[0]; n[0]++) {
|
---|
| 958 | pos[0] = (n[0] == N[0] ? 0 : n[0]);
|
---|
| 959 | for (n[1]=0; n[1] <= N[1]; n[1]++) {
|
---|
| 960 | pos[1] = (n[1] == N[1] ? 0 : n[1]);
|
---|
| 961 | for (n[2]=0; n[2] <= N[2]; n[2]++) {
|
---|
| 962 | pos[2] = (n[2] == N[2] ? 0 : n[2]);
|
---|
| 963 | // depending on RiemannTensor use, fill FileData::work
|
---|
| 964 | switch (Lat->RT.ActualUse) {
|
---|
| 965 | case inactive:
|
---|
| 966 | case standby:
|
---|
| 967 | if ((srcdens[pos[2]+N[2]*(pos[1]+N[1]*pos[0])]) > 0.)
|
---|
| 968 | dest[destpos] = log(1.0+(srcdens[pos[2]+N[2]*(pos[1]+N[1]*pos[0])]));
|
---|
| 969 | else
|
---|
| 970 | dest[destpos] = 0.;
|
---|
| 971 | destpos++;
|
---|
| 972 | break;
|
---|
| 973 | case active:
|
---|
| 974 | posf[0] = srcpos[0+NDIM*(pos[2]+N[2]*(pos[1]+N[1]*pos[0]))];
|
---|
| 975 | posf[1] = srcpos[1+NDIM*(pos[2]+N[2]*(pos[1]+N[1]*pos[0]))];
|
---|
| 976 | posf[2] = srcpos[2+NDIM*(pos[2]+N[2]*(pos[1]+N[1]*pos[0]))];
|
---|
| 977 | fac[0] = ((n[0]+N[0]*myPE)/(double)Nx);
|
---|
| 978 | fac[1] = (n[1]/(double)N[1]);
|
---|
| 979 | fac[2] = (n[2]/(double)N[2]);
|
---|
| 980 | RMat33Vec3(posd, Lat->RealBasis, fac);
|
---|
| 981 | posf[0] += posd[0];
|
---|
| 982 | posf[1] += posd[1];
|
---|
| 983 | posf[2] += posd[2];
|
---|
| 984 | if ((srcdens[pos[2]+N[2]*(pos[1]+N[1]*pos[0])]) > 0.)
|
---|
| 985 | posf[3] = log(1.0+(srcdens[pos[2]+N[2]*(pos[1]+N[1]*pos[0])]));
|
---|
| 986 | else
|
---|
| 987 | posf[3] = 0.;
|
---|
| 988 | dest[destpos+0] = posf[0];
|
---|
| 989 | dest[destpos+1] = posf[1];
|
---|
| 990 | dest[destpos+2] = posf[2];
|
---|
| 991 | dest[destpos+3] = posf[3];
|
---|
| 992 | destpos += 4;
|
---|
| 993 | break;
|
---|
| 994 | }
|
---|
| 995 | }
|
---|
| 996 | }
|
---|
| 997 | }
|
---|
| 998 | if (myPE) MPI_Send(dest, destpos, MPI_FLOAT, 0, OutputDensTag, P->Par.comm_ST_Psi);
|
---|
| 999 | }
|
---|
| 1000 |
|
---|
| 1001 | /** Combines prepared electronic Psis density and output to file.
|
---|
| 1002 | * If we are process 0, open file suffixdensdat (only when RiemannTensor is used) and suffixdenspos, receive
|
---|
| 1003 | * FileData::work logarithmic coefficients sent by the other processes in OutputOutVisDensity(), go through all
|
---|
| 1004 | * nodes and save the coefficient to file - again depending on RiemannTensor use - followed by FileData::PosTemp
|
---|
| 1005 | * (for y and z nodes), close file(s).
|
---|
| 1006 | * \param *P Problem at hand
|
---|
| 1007 | * \param me this ranks process in the Psi communcator ParallelSimulationData::me_comm_ST_Psi
|
---|
| 1008 | * \param Maxme number of processes in this Psi communcator ParallelSimulationData::Max_me_comm_ST_Psi
|
---|
| 1009 | */
|
---|
| 1010 | static void CombineOutVisDensity(struct Problem *P, const int me, const int Maxme)
|
---|
| 1011 | {
|
---|
| 1012 | int i,n[NDIM], N[NDIM];
|
---|
| 1013 | float posf[NDIM+1];
|
---|
| 1014 | float *source = (float *)P->Files.work;
|
---|
| 1015 | double posd[NDIM], fac[NDIM];
|
---|
| 1016 | float *Temp = (float *)P->Files.PosTemp;
|
---|
| 1017 | struct Lattice *Lat = &P->Lat;
|
---|
| 1018 | struct LatticeLevel *Lev0 = &Lat->Lev[0];
|
---|
| 1019 | float step = P->Files.OutVisStep;
|
---|
| 1020 | int No=0, destpos;
|
---|
| 1021 | FILE *DensityData, *DensityPos;
|
---|
| 1022 | int Nx = Lev0->Plan0.plan->N[0]+1;
|
---|
| 1023 | MPI_Status status;
|
---|
| 1024 | if (me) return; // if we are process 0!
|
---|
| 1025 | N[0] = Lev0->Plan0.plan->local_nx;
|
---|
| 1026 | N[1] = Lev0->Plan0.plan->N[1]+1;
|
---|
| 1027 | N[2] = Lev0->Plan0.plan->N[2]+1;
|
---|
| 1028 | // Open respective file depending on RiemannTensor use
|
---|
| 1029 | switch (Lat->RT.ActualUse) {
|
---|
| 1030 | case active:
|
---|
| 1031 | OpenFileNo(P, &DensityPos, suffixdenspos, (int)step, "wb",P->Call.out[ReadOut]);
|
---|
| 1032 | case inactive:
|
---|
| 1033 | case standby:
|
---|
| 1034 | OpenFileNo(P, &DensityData, suffixdensdat, (int)step, "wb",P->Call.out[ReadOut]);
|
---|
| 1035 | break;
|
---|
| 1036 | }
|
---|
| 1037 | // for all processes in the communicator
|
---|
| 1038 | for (i=0; i< Maxme; i++) {
|
---|
| 1039 | if (i) { // if process != 0, receive from this process
|
---|
| 1040 | /* MPI_Probe( i, OutputDensTag, P->Par.comm_ST_Psi, &status );*/
|
---|
| 1041 | switch (Lat->RT.ActualUse) {
|
---|
| 1042 | case inactive:
|
---|
| 1043 | case standby:
|
---|
| 1044 | MPI_Recv( source, N[0]*N[1]*N[2], MPI_FLOAT, i, OutputDensTag, P->Par.comm_ST_Psi, &status );
|
---|
| 1045 | break;
|
---|
| 1046 | case active:
|
---|
| 1047 | MPI_Recv( source, N[0]*N[1]*N[2]*4, MPI_FLOAT, i, OutputDensTag, P->Par.comm_ST_Psi, &status );
|
---|
| 1048 | break;
|
---|
| 1049 | }
|
---|
| 1050 | }
|
---|
| 1051 | destpos = 0;
|
---|
| 1052 | // go through all nodes and save the coefficient to file DensityData, depending on RiemannTensor
|
---|
| 1053 | for (n[0]=0; n[0] < N[0]; n[0]++) {
|
---|
| 1054 | for (n[1]=0; n[1] < N[1]; n[1]++) {
|
---|
| 1055 | for (n[2]=0; n[2] < N[2]; n[2]++) {
|
---|
| 1056 | switch (Lat->RT.ActualUse) {
|
---|
| 1057 | case inactive:
|
---|
| 1058 | case standby:
|
---|
| 1059 | posf[3] = source[destpos];
|
---|
| 1060 | destpos++;
|
---|
| 1061 | (void)fwrite(&posf[3], sizeof(float), (size_t)(1), DensityData);
|
---|
| 1062 | No++;
|
---|
| 1063 | if (i==0 && n[0] == 0)
|
---|
| 1064 | Temp[(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[3];
|
---|
| 1065 | break;
|
---|
| 1066 | case active:
|
---|
| 1067 | posf[0] = source[destpos+0];
|
---|
| 1068 | posf[1] = source[destpos+1];
|
---|
| 1069 | posf[2] = source[destpos+2];
|
---|
| 1070 | posf[3] = source[destpos+3];
|
---|
| 1071 | destpos += 4;
|
---|
| 1072 | (void)fwrite(posf, sizeof(float), (size_t)(NDIM), DensityPos);
|
---|
| 1073 | (void)fwrite(&posf[3], sizeof(float), (size_t)(1), DensityData);
|
---|
| 1074 | No++;
|
---|
| 1075 | if (i==0 && n[0] == 0) {
|
---|
| 1076 | fac[0] = ((n[0]+N[0]*i)/(double)(Nx-1));
|
---|
| 1077 | fac[1] = (n[1]/(double)(N[1]-1));
|
---|
| 1078 | fac[2] = (n[2]/(double)(N[2]-1));
|
---|
| 1079 | RMat33Vec3(posd, Lat->RealBasis, fac);
|
---|
| 1080 | posf[0] -= posd[0];
|
---|
| 1081 | posf[1] -= posd[1];
|
---|
| 1082 | posf[2] -= posd[2];
|
---|
| 1083 | fac[0] = ((Nx-1)/(double)(Nx-1));
|
---|
| 1084 | fac[1] = (n[1]/(double)(N[1]-1));
|
---|
| 1085 | fac[2] = (n[2]/(double)(N[2]-1));
|
---|
| 1086 | RMat33Vec3(posd, Lat->RealBasis, fac);
|
---|
| 1087 | posf[0] += posd[0];
|
---|
| 1088 | posf[1] += posd[1];
|
---|
| 1089 | posf[2] += posd[2];
|
---|
| 1090 | Temp[0+(NDIM+1)*(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[0];
|
---|
| 1091 | Temp[1+(NDIM+1)*(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[1];
|
---|
| 1092 | Temp[2+(NDIM+1)*(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[2];
|
---|
| 1093 | Temp[3+(NDIM+1)*(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[3];
|
---|
| 1094 | }
|
---|
| 1095 | break;
|
---|
| 1096 | }
|
---|
| 1097 | }
|
---|
| 1098 | }
|
---|
| 1099 | }
|
---|
| 1100 | }
|
---|
| 1101 | n[0] = N[0];
|
---|
| 1102 | for (n[1]=0; n[1] < N[1]; n[1]++) {
|
---|
| 1103 | for (n[2]=0; n[2] < N[2]; n[2]++) {
|
---|
| 1104 | switch (Lat->RT.ActualUse) {
|
---|
| 1105 | case inactive:
|
---|
| 1106 | case standby:
|
---|
| 1107 | (void)fwrite(&Temp[n[2]+N[2]*(n[1])], sizeof(float), (size_t)(1), DensityData);
|
---|
| 1108 | No++;
|
---|
| 1109 | break;
|
---|
| 1110 | case active:
|
---|
| 1111 | (void)fwrite(&Temp[(NDIM+1)*(n[2]+N[2]*(n[1]))], sizeof(float), (size_t)(NDIM), DensityPos);
|
---|
| 1112 | (void)fwrite(&Temp[(NDIM+1)*(n[2]+N[2]*(n[1]))+3], sizeof(float), (size_t)(1), DensityData);
|
---|
| 1113 | No++;
|
---|
| 1114 | break;
|
---|
| 1115 | }
|
---|
| 1116 | }
|
---|
| 1117 | }
|
---|
| 1118 | if (No != Nx*N[1]*N[2]) Error(SomeError,"CombineOutVisDensity: No != points");
|
---|
| 1119 | switch (Lat->RT.ActualUse) {
|
---|
| 1120 | case active:
|
---|
| 1121 | fclose(DensityPos);
|
---|
| 1122 | case inactive:
|
---|
| 1123 | case standby:
|
---|
| 1124 | fclose(DensityData);
|
---|
| 1125 | break;
|
---|
| 1126 | }
|
---|
| 1127 | }
|
---|
| 1128 |
|
---|
| 1129 | /** Main output electronic Psis density for OpenDX.
|
---|
| 1130 | * If FileData::MeOutVis is set, calls OutputOutVisDensity() followed by CombineOutVisDensity().
|
---|
| 1131 | * Beforehand CalculateOutVisDensityPos() is called if RiemannTensor is used.
|
---|
| 1132 | * \param *P Problem at hand
|
---|
| 1133 | * \param *src_dens Pointer to DensityArray which is to be displayed
|
---|
| 1134 | */
|
---|
| 1135 | static void OutVisDensity(struct Problem *P, fftw_real *src_dens)
|
---|
| 1136 | {
|
---|
| 1137 | if (!P->Files.MeOutVis) return;
|
---|
| 1138 | if (P->Lat.RT.ActualUse == active) CalculateOutVisDensityPos(&P->Lat, &P->Files/*, P->Lat.FactorDensityR, P->Lat.FactorDensityC*/);
|
---|
| 1139 | OutputOutVisDensity(P, P->Par.me_comm_ST_Psi, src_dens);
|
---|
| 1140 | /* Achtung hier: P->Files.work (RT->TempC, Dens->DensityCArray[TempDensity]) fuer myPE == 0 nicht veraendern !!! */
|
---|
| 1141 | CombineOutVisDensity(P, P->Par.me_comm_ST_Psi, P->Par.Max_me_comm_ST_Psi);
|
---|
| 1142 | }
|
---|
| 1143 |
|
---|
| 1144 | /** Opening and Initializing of output measurement files.
|
---|
| 1145 | * If this is process 0, opens and writes top line of FileData::ForcesFile, FileData::EnergyFile.
|
---|
| 1146 | * and sets FileData::MeOutVis and FileData::MeOutMes (if output desired) to 1, otherwise 0.
|
---|
| 1147 | * \param *P Problem at hand
|
---|
| 1148 | */
|
---|
| 1149 | void InitOutputFiles(struct Problem *P)
|
---|
| 1150 | {
|
---|
| 1151 | struct FileData *F = &P->Files;
|
---|
| 1152 | F->ForcesFile = NULL;
|
---|
| 1153 | F->EnergyFile = NULL;
|
---|
| 1154 | F->HamiltonianFile = NULL;
|
---|
| 1155 | F->MinimisationFile = NULL;
|
---|
| 1156 | F->SpreadFile = NULL;
|
---|
| 1157 | // process 0 ?
|
---|
| 1158 | F->MeOutVis = ((P->Par.my_color_comm_ST == 0 && P->Par.my_color_comm_ST_Psi == 0 && F->DoOutVis) ? 1 : 0);
|
---|
| 1159 | F->MeOutCurr = ((P->Par.my_color_comm_ST == 0 && P->Par.my_color_comm_ST_Psi == 0 && F->DoOutCurr) ? 1 : 0);
|
---|
| 1160 | F->MeOutMes = ((P->Par.me == 0 && F->DoOutMes) ? 1 : 0);
|
---|
| 1161 | OpenFile(P, &F->HamiltonianFile, suffixhamiltonianall, "w",P->Call.out[ReadOut]);
|
---|
| 1162 |
|
---|
| 1163 | if (!F->MeOutMes) return;
|
---|
| 1164 | OpenFile(P, &F->ForcesFile, suffixforcesall, "w",P->Call.out[ReadOut]);
|
---|
| 1165 | if (F->ForcesFile == NULL) fprintf(stderr,"Error opening ForcesFile\n");
|
---|
| 1166 | // write header of forces file
|
---|
| 1167 | fprintf(F->ForcesFile, "%s\t%s\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\n",
|
---|
| 1168 | "Type", "No",
|
---|
| 1169 | "Pos0", "Pos1", "Pos2",
|
---|
| 1170 | "Total0", "Total1", "Total2",
|
---|
| 1171 | "Local0", "Local1", "Local2",
|
---|
| 1172 | "NLocal0", "NLocal1", "NLocal2",
|
---|
| 1173 | "Ewald0", "Ewald1", "Ewald2");
|
---|
| 1174 | OpenFile(P, &F->EnergyFile, suffixenergyall, "w",P->Call.out[ReadOut]);
|
---|
| 1175 | if (F->EnergyFile == NULL) fprintf(stderr,"Error opening EnergyFile\n");
|
---|
| 1176 | // write header of energy file
|
---|
| 1177 | if (P->R.DoUnOccupied) {
|
---|
| 1178 | fprintf(F->EnergyFile, "%s\t\t%s\t\t%s\t%s\t\t%s\t%s\t\t%s\t%s\t%s\t\t%s\t\t%s\t%s\t\t%s\t\t%s\t\t%s\n",
|
---|
| 1179 | "Time",
|
---|
| 1180 | "Total",
|
---|
| 1181 | "Total+Gap",
|
---|
| 1182 | "Kinetic", "NonLocal",
|
---|
| 1183 | "GapPsi",
|
---|
| 1184 | "Correlation", "Exchange",
|
---|
| 1185 | "Pseudo", "Hartree",
|
---|
| 1186 | "GapDensity",
|
---|
| 1187 | "-Gauss",
|
---|
| 1188 | "Ewald",
|
---|
| 1189 | "IonKin",
|
---|
| 1190 | "ETotal");
|
---|
| 1191 | } else {
|
---|
| 1192 | fprintf(F->EnergyFile, "%s\t\t%s\t\t%s\t\t%s\t%s\t%s\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\n",
|
---|
| 1193 | "Time",
|
---|
| 1194 | "Total",
|
---|
| 1195 | "Kinetic", "NonLocal",
|
---|
| 1196 | "Correlation", "Exchange",
|
---|
| 1197 | "Pseudo", "Hartree",
|
---|
| 1198 | "-Gauss",
|
---|
| 1199 | "Ewald",
|
---|
| 1200 | "IonKin",
|
---|
| 1201 | "ETotal");
|
---|
| 1202 | }
|
---|
| 1203 | OpenFile(P, &F->MinimisationFile, suffixminall, "w",P->Call.out[ReadOut]);
|
---|
| 1204 | if (F->MinimisationFile == NULL) fprintf(stderr,"Error opening MinimsationFile\n");
|
---|
| 1205 | // write header of minimsation file
|
---|
| 1206 | fprintf(F->MinimisationFile, "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n","Step", "Psi", "PsiStep", "TE", "ATE", "delta", "dEdt0", "ddEdtt0");
|
---|
| 1207 | }
|
---|
| 1208 |
|
---|
| 1209 | /** Closes all output measurement files.
|
---|
| 1210 | * Closes FileData::ForcesFile and FileData::EnergyFile
|
---|
| 1211 | * \param *P Problem at hand
|
---|
| 1212 | */
|
---|
| 1213 | void CloseOutputFiles(struct Problem *P)
|
---|
| 1214 | {
|
---|
| 1215 | struct FileData *F = &P->Files;
|
---|
| 1216 | if (!F->MeOutMes) return;
|
---|
| 1217 | if (F->ForcesFile != NULL) fclose(F->ForcesFile); // only they've been opened (thus not pointing to NULL)
|
---|
| 1218 | if (F->EnergyFile != NULL) fclose(F->EnergyFile);
|
---|
| 1219 | if (F->HamiltonianFile != NULL) fclose(F->HamiltonianFile);
|
---|
| 1220 | if (F->MinimisationFile != NULL) fclose(F->MinimisationFile);
|
---|
| 1221 | if (F->SpreadFile != NULL) fclose(F->SpreadFile);
|
---|
| 1222 | }
|
---|
| 1223 |
|
---|
| 1224 | /** Initialization of Problem::FileData structure for visual output.
|
---|
| 1225 | * If this is process 0 (and OutVis desired), allocate memory for FileData::PosTemp, FileData::work,
|
---|
| 1226 | * set the entries of FileData all to their corresponding values from RiemannTensor,
|
---|
| 1227 | * FileData::OutVisStep to zero.
|
---|
| 1228 | * \param *P Problem at hand
|
---|
| 1229 | */
|
---|
| 1230 | void InitOutVisArray(struct Problem *P)
|
---|
| 1231 | {
|
---|
| 1232 | struct FileData *F = &P->Files;
|
---|
| 1233 | struct Lattice *Lat = &P->Lat;
|
---|
| 1234 | struct RiemannTensor *RT = &Lat->RT;
|
---|
| 1235 | struct LatticeLevel *Lev0 = &Lat->Lev[0];
|
---|
| 1236 | F->OutputPosType = NULL;
|
---|
| 1237 | if (!F->MeOutVis) return;
|
---|
| 1238 | F->OutVisStep = 0;
|
---|
| 1239 | F->PosTemp = (fftw_complex *)
|
---|
| 1240 | Malloc(sizeof(float)*(Lev0->Plan0.plan->N[1]+1)*(Lev0->Plan0.plan->N[2]+1)*
|
---|
| 1241 | ((Lat->RT.Use != UseRT ? 0 : NDIM)+1), "InitOutVisArray: TempC");
|
---|
| 1242 | F->work = (fftw_complex *)Lev0->Dens->DensityCArray[TempDensity];
|
---|
| 1243 | if (Lat->RT.Use != UseRT) return;
|
---|
| 1244 | F->TotalSize = RT->TotalSize[RTAIRT]/NDIM;
|
---|
| 1245 | F->LocalSizeR = RT->LocalSizeR[RTAiGcg];
|
---|
| 1246 | F->LocalSizeC = RT->LocalSizeC[RTAiGcg];
|
---|
| 1247 | F->MaxNUp = RT->MaxNUp[RTPFRto0];
|
---|
| 1248 | F->PosC = RT->DensityC[RTAiGcg];
|
---|
| 1249 | F->PosR = (fftw_real *)F->PosC;
|
---|
| 1250 | F->work = RT->TempC;
|
---|
| 1251 | F->PosFactor = RT->PosFactor[RTPFRto0];
|
---|
| 1252 | }
|
---|
| 1253 |
|
---|
| 1254 | static const char suffixionfor[] = ".ions.force"; //!< Suffix for ion forces file
|
---|
| 1255 | static const char suffixionZ[] = ".ions.datZ"; //!< Suffix for ion datZ file
|
---|
| 1256 | static const char suffixionpos[] = ".ions.pos"; //!< Suffix for ion positions file
|
---|
| 1257 | static const char suffixiondx[] = ".ions.dx"; //!< Suffix for ions dx file
|
---|
| 1258 | static const char suffixiondoc[] = ".ions.doc"; //!< Suffix for ions doc file
|
---|
| 1259 | static const char suffixsrciondoc[] = ".srcion.doc"; //!< Suffix for state ions doc file
|
---|
| 1260 | static const char suffixsrciondat[] = ".srcion.data"; //!< Suffix for state ions data file
|
---|
| 1261 |
|
---|
| 1262 | /** Output current ions state.
|
---|
| 1263 | * If this is process0, open file suffixsrciondoc for writing, output Ions::Max_Types and
|
---|
| 1264 | * Ions::Max_IonsOfType of each type - each in a new line - closes it.
|
---|
| 1265 | * Then opens suffixsrciondat for binary writing, outputs Lattice:RealBasis vectors and
|
---|
| 1266 | * position IonType::R and speed IonType::U, closes it.
|
---|
| 1267 | * \param *P Problem at hand
|
---|
| 1268 | * \note This is the ionic counterpart to the elecontric OutputSrcPsiDensity(), storing a so far made
|
---|
| 1269 | * calculation to file.
|
---|
| 1270 | */
|
---|
| 1271 | static void OutSrcIons(struct Problem *P)
|
---|
| 1272 | {
|
---|
| 1273 | struct Ions *I = &P->Ion;
|
---|
| 1274 | double *U, *pos;
|
---|
| 1275 | double data[2*NDIM];
|
---|
| 1276 | int is,ia,i;
|
---|
| 1277 | FILE *SrcIonDoc, *SrcIonData;
|
---|
| 1278 |
|
---|
| 1279 | if (!(P->Par.me == 0)) return;
|
---|
| 1280 |
|
---|
| 1281 | // output of ion types and numbers per type
|
---|
| 1282 | OpenFile(P, &SrcIonDoc, suffixsrciondoc, "w",P->Call.out[ReadOut]);
|
---|
| 1283 | fprintf(SrcIonDoc,"%i\n", I->Max_Types);
|
---|
| 1284 | for (is=0; is < I->Max_Types; is++)
|
---|
| 1285 | fprintf(SrcIonDoc,"%i\n", I->I[is].Max_IonsOfType);
|
---|
| 1286 | fclose(SrcIonDoc);
|
---|
| 1287 |
|
---|
| 1288 | OpenFile(P, &SrcIonData, suffixsrciondat, "wb",P->Call.out[ReadOut]);
|
---|
| 1289 | (void)fwrite(P->Lat.RealBasis, sizeof(double), (size_t)(NDIM_NDIM), SrcIonData);
|
---|
| 1290 | for (is=0; is < I->Max_Types; is++) {
|
---|
| 1291 | for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
|
---|
| 1292 | U = &I->I[is].U[NDIM*ia];
|
---|
| 1293 | pos = &I->I[is].R[NDIM*ia];
|
---|
| 1294 | for (i=0;i<NDIM;i++) {
|
---|
| 1295 | data[i] = pos[i];
|
---|
| 1296 | data[i+NDIM] = U[i];
|
---|
| 1297 | }
|
---|
| 1298 | (void)fwrite(&data, sizeof(double), (size_t)(2*NDIM), SrcIonData);
|
---|
| 1299 | }
|
---|
| 1300 | }
|
---|
| 1301 | fclose(SrcIonData);
|
---|
| 1302 | }
|
---|
| 1303 |
|
---|
| 1304 | /** Read ions state from a file.
|
---|
| 1305 | * Reads the suffixsrciondoc file and checks it against the current state in Ions regarding
|
---|
| 1306 | * IonType::MaxTypes and IonType::Max_IonsOfType, closes it.
|
---|
| 1307 | * Afterwards, opens suffixsrciondat for binary reading, retrieves the basis checking it against
|
---|
| 1308 | * Problem::Lattice::RealBasis. If ok, reads position IonType::R and speed IonType::U, closes it.
|
---|
| 1309 | * \param *P Problem at hand
|
---|
| 1310 | * \return 1 - successful, 0 - failure
|
---|
| 1311 | * \note This is the ionic counterpart to the electronic ReadSrcPsiDensity(), see also OutSrcIons().
|
---|
| 1312 | */
|
---|
| 1313 | int ReadSrcIons(struct Problem *P)
|
---|
| 1314 | {
|
---|
| 1315 | struct Ions *I = &P->Ion;
|
---|
| 1316 | double *U, *pos;
|
---|
| 1317 | double data[2*NDIM];
|
---|
| 1318 | int is,ia,i;
|
---|
| 1319 | int Max_Types;
|
---|
| 1320 | int *Max_IonsOfType = NULL;
|
---|
| 1321 | double RealBasis[NDIM_NDIM];
|
---|
| 1322 | FILE *SrcIonDoc, *SrcIonData;
|
---|
| 1323 | // read the doc file and check
|
---|
| 1324 | if (OpenFile(P, &SrcIonDoc, suffixsrciondoc, "r",P->Call.out[ReadOut])) {
|
---|
| 1325 | if (fscanf(SrcIonDoc,"%i", &Max_Types) != 1)
|
---|
| 1326 | //Error(SomeError, "ReadSrcIons: read error");
|
---|
| 1327 | return 0;
|
---|
| 1328 | if (Max_Types != I->Max_Types)
|
---|
| 1329 | //Error(SomeError, "ReadSrcIons: srcion file does not fit to system, MaxTypes");
|
---|
| 1330 | return 0;
|
---|
| 1331 | Max_IonsOfType = (int *) Malloc(Max_Types*sizeof(int), "ReadSrcIons: Max_IonsOfType");
|
---|
| 1332 | for (is=0; is < Max_Types; is++) {
|
---|
| 1333 | if (fscanf(SrcIonDoc,"%i", &Max_IonsOfType[is]) != 1)
|
---|
| 1334 | //Error(SomeError, "ReadSrcIons: read error");
|
---|
| 1335 | return 0;
|
---|
| 1336 | if (Max_IonsOfType[is] != I->I[is].Max_IonsOfType)
|
---|
| 1337 | //Error(SomeError, "ReadSrcIons: srcion file does not fit to system, Max_IonsOfType");
|
---|
| 1338 | return 0;
|
---|
| 1339 | }
|
---|
| 1340 | fclose(SrcIonDoc);
|
---|
| 1341 | // read basis, then positions and speeds of ions
|
---|
| 1342 | if (OpenFile(P, &SrcIonData, suffixsrciondat, "rb",P->Call.out[ReadOut])) {
|
---|
| 1343 | if (fread(RealBasis, sizeof(double), (size_t)(NDIM_NDIM), SrcIonData) != NDIM_NDIM)
|
---|
| 1344 | //Error(SomeError, "ReadSrcIons: read error");
|
---|
| 1345 | return 0;
|
---|
| 1346 | for (i=0; i < NDIM_NDIM; i++)
|
---|
| 1347 | if (RealBasis[i] != P->Lat.RealBasis[i])
|
---|
| 1348 | //Error(SomeError, "ReadSrcIons: srcion file does not fit to system, RealBasis");
|
---|
| 1349 | return 0;
|
---|
| 1350 | for (is=0; is < I->Max_Types; is++) {
|
---|
| 1351 | for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
|
---|
| 1352 | if (fread(&data, sizeof(double), (size_t)(2*NDIM), SrcIonData) != 2*NDIM)
|
---|
| 1353 | //Error(SomeError, "ReadSrcIons: read error");
|
---|
| 1354 | return 0;
|
---|
| 1355 | U = &I->I[is].U[NDIM*ia];
|
---|
| 1356 | pos = &I->I[is].R[NDIM*ia];
|
---|
| 1357 | for (i=0;i<NDIM;i++) {
|
---|
| 1358 | pos[i] = data[i];
|
---|
| 1359 | U[i] = data[i+NDIM];
|
---|
| 1360 | }
|
---|
| 1361 | }
|
---|
| 1362 | }
|
---|
| 1363 | fclose(SrcIonData);
|
---|
| 1364 | }
|
---|
[64fa9e] | 1365 | if (Max_IonsOfType != NULL) Free(Max_IonsOfType, "Max_IonsOfType");
|
---|
[a0bcf1] | 1366 | }
|
---|
| 1367 | return 1;
|
---|
| 1368 | }
|
---|
| 1369 |
|
---|
| 1370 | /** Output of ion doc, dx, forces and positions file for OpenDX.
|
---|
| 1371 | * If this is process 0,
|
---|
| 1372 | * open, fill and close IonDoc file suffixiondoc,
|
---|
| 1373 | * open, fill for each FileData::OutVisStep and close IonDX file suffixiondx
|
---|
| 1374 | * for every
|
---|
| 1375 | * open suffixionfor, suffixionpos (and suffixionZ in case of only FileData::OutVisStep), fill
|
---|
| 1376 | * them with the ion forces IonType::FIon and positions IonType::R of each type and each ion per type,
|
---|
| 1377 | * close them all.
|
---|
| 1378 | * \param *P Problem at hand
|
---|
| 1379 | */
|
---|
| 1380 | static void OutVisIons(struct Problem *P)
|
---|
| 1381 | {
|
---|
| 1382 | struct Ions *I = &P->Ion;
|
---|
| 1383 | struct FileData *F = &P->Files;
|
---|
| 1384 | int i,is,ia;
|
---|
| 1385 | double *fion, *pos;
|
---|
| 1386 | float data[6]; // holds temporarily twice NDIM values as write buffer
|
---|
| 1387 | int Z;
|
---|
| 1388 | char *datnamef, *datnameZ, *posname;
|
---|
| 1389 | FILE *IonsDataF, *IonsDataZ, *IonsPos, *IonsDoc, *IonsDx;
|
---|
| 1390 | if (!P->Files.MeOutVis && P->Par.me == 0) return;
|
---|
| 1391 | // generate file names
|
---|
| 1392 | datnamef = (char*)
|
---|
| 1393 | malloc(strlen(P->Files.mainname)+strlen(suffixionfor) + 1);
|
---|
| 1394 | sprintf(datnamef, "%s%s", P->Files.mainname, suffixionfor);
|
---|
| 1395 | datnameZ = (char*)
|
---|
| 1396 | malloc(strlen(P->Files.mainname)+strlen(suffixionZ) + 1);
|
---|
| 1397 | sprintf(datnameZ, "%s%s", P->Files.mainname, suffixionZ);
|
---|
| 1398 | posname = (char*)
|
---|
| 1399 | malloc(strlen(P->Files.mainname)+strlen(suffixionpos) + 1);
|
---|
| 1400 | sprintf(posname, "%s%s", P->Files.mainname, suffixionpos);
|
---|
| 1401 | // open, fill and close doc file
|
---|
| 1402 | if (OpenFile(P, &IonsDoc, suffixiondoc, "w",P->Call.out[ReadOut])) {
|
---|
| 1403 | fprintf(IonsDoc,"IonsPos file = %s.####\n", posname);
|
---|
| 1404 | fprintf(IonsDoc,"IonsForce file = %s.####\n", datnamef);
|
---|
| 1405 | fprintf(IonsDoc,"format = ieee float (Bytes %lu) %s = Force(3)\n",(unsigned long) sizeof(float),msb);
|
---|
| 1406 | fprintf(IonsDoc,"IonsZ file = %s.####\n", datnameZ);
|
---|
| 1407 | fprintf(IonsDoc,"format = int (Bytes %lu) %s = Z(1)\n",(unsigned long) sizeof(int),msb);
|
---|
| 1408 | fprintf(IonsDoc,"points = %i\n", I->Max_TotalIons);
|
---|
| 1409 | fprintf(IonsDoc,"majority = row\n");
|
---|
| 1410 | fprintf(IonsDoc,"TimeSeries = %i\n",F->OutVisStep+1);
|
---|
| 1411 | fclose(IonsDoc);
|
---|
| 1412 | }
|
---|
| 1413 | // open dx file and fill it with each output step, close it
|
---|
| 1414 | if (OpenFile(P, &IonsDx, suffixiondx, "w",P->Call.out[ReadOut])) {
|
---|
| 1415 | for (i=0; i < F->OutVisStep+1; i++) {
|
---|
| 1416 | if (i==0) {
|
---|
| 1417 | /* fprintf(IonsDx,"object \"ioncon\" class array type int rank 1 shape 2 items 0 data follows\nattribute \"element type\" string \"lines\"\nattribute \"ref\" string \"positions\"\n\n"); */
|
---|
| 1418 | fprintf(IonsDx,"object \"iondatZ\" class array type int rank 0 items %i %s binary\n",I->Max_TotalIons,msb);
|
---|
| 1419 | fprintf(IonsDx,"data file %s,0\n",datnameZ);
|
---|
| 1420 | fprintf(IonsDx,"attribute \"dep\" string \"positions\"\n\n");
|
---|
| 1421 | }
|
---|
| 1422 |
|
---|
| 1423 | fprintf(IonsDx,"object \"ionpos.%04i\" class array type float rank 1 shape 3 items %i %s binary\n",i,I->Max_TotalIons,msb);
|
---|
| 1424 | fprintf(IonsDx,"data file %s.%04i,0\n\n",posname,i);
|
---|
| 1425 |
|
---|
| 1426 | fprintf(IonsDx,"object \"iondatF.%04i\" class array type float rank 1 shape 3 items %i %s binary\n",i,I->Max_TotalIons,msb);
|
---|
| 1427 | fprintf(IonsDx,"data file %s.%04i,0\n",datnamef,i);
|
---|
| 1428 | fprintf(IonsDx,"attribute \"dep\" string \"positions\"\n\n");
|
---|
| 1429 |
|
---|
| 1430 | fprintf(IonsDx,"object \"ionobjF.%04i\" class field\n",i);
|
---|
| 1431 | fprintf(IonsDx,"component \"positions\" \"ionpos.%04i\"\n",i);
|
---|
| 1432 | fprintf(IonsDx,"component \"data\" \"iondatF.%04i\"\n",i);
|
---|
| 1433 | /*fprintf(IonsDx,"component \"connections\" \"ioncon\"\n\n");*/
|
---|
| 1434 |
|
---|
| 1435 | fprintf(IonsDx,"object \"ionobjZ.%04i\" class field\n",i);
|
---|
| 1436 | fprintf(IonsDx,"component \"positions\" \"ionpos.%04i\"\n",i);
|
---|
| 1437 | fprintf(IonsDx,"component \"data\" \"iondatZ\"\n");
|
---|
| 1438 | /* fprintf(IonsDx,"component \"connections\" \"ioncon\"\n\n");*/
|
---|
| 1439 | }
|
---|
| 1440 | fprintf(IonsDx,"\nobject \"ionseriesF\" class series\n");
|
---|
| 1441 | for (i=0; i < F->OutVisStep+1; i++)
|
---|
| 1442 | fprintf(IonsDx,"member %i \"ionobjF.%04i\" position %f\n",i,i,(float)i);
|
---|
| 1443 | fprintf(IonsDx,"\nobject \"ionseriesZ\" class series\n");
|
---|
| 1444 | for (i=0; i < F->OutVisStep+1; i++)
|
---|
| 1445 | fprintf(IonsDx,"member %i \"ionobjZ.%04i\" position %f\n",i,i,(float)i);
|
---|
| 1446 | fprintf(IonsDx,"end\n");
|
---|
| 1447 | fclose(IonsDx);
|
---|
| 1448 | }
|
---|
| 1449 | free(datnamef);
|
---|
| 1450 | free(datnameZ);
|
---|
| 1451 | free(posname);
|
---|
| 1452 | // open IonForces, IonZ and IonPosition file, write forces respectively positions for each ion of each type, close them
|
---|
| 1453 | if (OpenFileNo(P, &IonsDataF, suffixionfor, F->OutVisStep, "wb",P->Call.out[ReadOut])) {
|
---|
| 1454 | if (F->OutVisStep == 0)
|
---|
| 1455 | OpenFile(P, &IonsDataZ, suffixionZ, "wb",P->Call.out[ReadOut]);
|
---|
| 1456 | if (OpenFileNo(P, &IonsPos, suffixionpos, F->OutVisStep, "wb",P->Call.out[ReadOut])) {
|
---|
| 1457 | for (is=0; is < I->Max_Types; is++) {
|
---|
| 1458 | for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
|
---|
| 1459 | fion = &I->I[is].FIon[NDIM*ia];
|
---|
| 1460 | pos = &I->I[is].R[NDIM*ia];
|
---|
| 1461 | for (i=0;i<3;i++) {
|
---|
| 1462 | data[i+3] = fion[i];
|
---|
| 1463 | data[i] = pos[i];
|
---|
| 1464 | }
|
---|
| 1465 | Z = I->I[is].Z;
|
---|
| 1466 | if (fwrite(&data[0],sizeof(float), (size_t)(3),IonsPos) != 3)
|
---|
| 1467 | Error(FileOpenParams, "Error writing ion positions!");
|
---|
| 1468 | if (F->OutVisStep == 0) (void)fwrite(&Z,sizeof(int), (size_t)(1),IonsDataZ);
|
---|
| 1469 | if (fwrite(&data[3],sizeof(float), (size_t)(3),IonsDataF) != 3)
|
---|
| 1470 | Error(FileOpenParams, "Error writing ionic forces!");
|
---|
| 1471 | }
|
---|
| 1472 | }
|
---|
| 1473 | fclose(IonsPos);
|
---|
| 1474 | }
|
---|
| 1475 | if (F->OutVisStep == 0)
|
---|
| 1476 | fclose(IonsDataZ);
|
---|
| 1477 | fclose(IonsDataF);
|
---|
| 1478 | }
|
---|
| 1479 | }
|
---|
| 1480 |
|
---|
| 1481 | /** Output electronic density and ion state files with so far made calculations.
|
---|
| 1482 | * If CallOptions::WriteSrcFiles is set, OutputSrcPsiDensity() and OutSrcIons() are called.
|
---|
| 1483 | * \param *P Problem at hand
|
---|
| 1484 | * \param type which PsiTypeTag should be put to file
|
---|
| 1485 | */
|
---|
| 1486 | void OutputVisSrcFiles(struct Problem *P, enum PsiTypeTag type)
|
---|
| 1487 | {
|
---|
| 1488 | if (P->Call.WriteSrcFiles) {
|
---|
| 1489 | if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Write srcpsi to disk\n", P->Par.me);
|
---|
| 1490 | OutputSrcPsiDensity(P, type);
|
---|
| 1491 | if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Write srcion to disk\n", P->Par.me);
|
---|
| 1492 | OutSrcIons(P);
|
---|
| 1493 | }
|
---|
| 1494 | // if (!P->Files.MeOutVis) return;
|
---|
| 1495 | // P->Files.OutputPosType =
|
---|
| 1496 | // Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+1),"OutputVis");
|
---|
| 1497 | // P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse;
|
---|
| 1498 | // CreateDensityOutputGeneral(P, P->Par.me_comm_ST_Psi);
|
---|
| 1499 | // OutVisDensity(P);
|
---|
| 1500 | // OutVisIons(P);
|
---|
| 1501 | // if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Written OutVisStep %i to disk\n", P->Par.me, P->Files.OutVisStep);
|
---|
| 1502 | // /* P->Files.OutVisStep++; Genau ebend nicht hochzaehlen - wird immer ueberschrieben */
|
---|
| 1503 | }
|
---|
| 1504 |
|
---|
| 1505 | /** Main output total electronic density and ion data for OpenDX.
|
---|
| 1506 | * Calls subsequently preparing CreateDensityOutputGeneral(), then output of electronic
|
---|
| 1507 | * densities OutVisDensity() and ion data OutVisIons(), increasing finally FileData::OutVisStep.
|
---|
| 1508 | * \param *P Problem at hand
|
---|
| 1509 | * \param *srcdens Pointer to DensityArray which is to be displayed
|
---|
| 1510 | * \note Output is made only RunStruct::OutVisStep steps and if FileData::MeOutVis is set.
|
---|
| 1511 | */
|
---|
| 1512 | void OutputVis(struct Problem *P, fftw_real *srcdens)
|
---|
| 1513 | {
|
---|
| 1514 | if (!P->Files.MeOutVis) return;
|
---|
| 1515 | P->Files.OutputPosType = (enum ModeType *) Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+1),"OutputVis");
|
---|
| 1516 | P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse;
|
---|
| 1517 |
|
---|
| 1518 | CreateDensityOutputGeneral(P, P->Par.me_comm_ST_Psi);
|
---|
| 1519 | OutVisDensity(P, srcdens);
|
---|
| 1520 | OutVisIons(P);
|
---|
| 1521 | if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Written OutVisStep %i to disk\n", P->Par.me, P->Files.OutVisStep);
|
---|
| 1522 | P->Files.OutVisStep++;
|
---|
| 1523 | }
|
---|
| 1524 |
|
---|
| 1525 | /** Output of each current density component for OpenDX.
|
---|
| 1526 | * \param *P Problem at hand
|
---|
| 1527 | * \note Output is made only RunStruct::OutVisStep steps and if FileData::MeOutVis is set.
|
---|
| 1528 | */
|
---|
| 1529 | void OutputCurrentDensity(struct Problem *P)
|
---|
| 1530 | {
|
---|
| 1531 | int index, i, r;
|
---|
| 1532 | fftw_real *density = P->R.Lev0->Dens->DensityArray[ActualDensity];
|
---|
| 1533 | fftw_real *CurrentDensity[NDIM*NDIM];
|
---|
| 1534 | if (!P->Files.MeOutCurr) return;
|
---|
| 1535 |
|
---|
| 1536 | P->Files.OutputPosType = (enum ModeType *) Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+(1)*NDIM),"OutputVis");
|
---|
| 1537 | for (i=0;i<(1)*NDIM;i++)
|
---|
| 1538 | P->Files.OutputPosType[P->Files.OutVisStep+i] = P->Lat.RT.ActualUse;
|
---|
| 1539 | fprintf(stderr,"(%i) OutVisStep %i, OutputPosType %p\n",P->Par.me, P->Files.OutVisStep, P->Files.OutputPosType);
|
---|
| 1540 |
|
---|
| 1541 | // due to preprocessor values we can't put the following stuff into a loop
|
---|
| 1542 | CurrentDensity[0] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity0];
|
---|
| 1543 | CurrentDensity[1] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity1];
|
---|
| 1544 | CurrentDensity[2] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity2];
|
---|
| 1545 | CurrentDensity[3] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity3];
|
---|
| 1546 | CurrentDensity[4] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity4];
|
---|
| 1547 | CurrentDensity[5] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity5];
|
---|
| 1548 | CurrentDensity[6] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity6];
|
---|
| 1549 | CurrentDensity[7] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity7];
|
---|
| 1550 | CurrentDensity[8] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity8];
|
---|
| 1551 |
|
---|
| 1552 | // output current density, not vector component
|
---|
| 1553 | for (index=0;index<NDIM;index++) {
|
---|
| 1554 | // evaluate euclidian norm for given B component
|
---|
| 1555 | SetArrayToDouble0((double *)density,P->R.Lev0->Dens->TotalSize*2); // reset
|
---|
| 1556 | for (r=0;r<P->R.Lev0->Dens->LocalSizeR;r++) {
|
---|
| 1557 | for (i=0;i<NDIM;i++)
|
---|
| 1558 | density[r] += CurrentDensity[i + index*NDIM][r]*CurrentDensity[i + index*NDIM][r];
|
---|
| 1559 | density[r] = sqrt(density[r]);
|
---|
| 1560 | }
|
---|
| 1561 | // output
|
---|
| 1562 | CreateDensityOutputGeneral(P, P->Par.me_comm_ST_Psi);
|
---|
| 1563 | OutVisDensity(P, density);
|
---|
| 1564 | OutVisIons(P);
|
---|
| 1565 | if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Written OutVisStep %i to disk\n", P->Par.me, P->Files.OutVisStep);
|
---|
| 1566 | P->Files.OutVisStep++;
|
---|
| 1567 | }
|
---|
| 1568 | }
|
---|
| 1569 |
|
---|
| 1570 | /** Output each orbital in a certain step order for OpenDX.
|
---|
| 1571 | * \param *P Problem at hand
|
---|
| 1572 | * \param offset from which step do we start
|
---|
| 1573 | * \param increment by which increment do we advance step-wise
|
---|
| 1574 | * \param type Only PsiTypeTag orbitals are displayed
|
---|
| 1575 | */
|
---|
| 1576 | void OutputVisAllOrbital(struct Problem *P, int offset, int increment, enum PsiTypeTag type) {
|
---|
| 1577 | struct Lattice *Lat = &P->Lat;
|
---|
| 1578 | struct Psis *Psi = &Lat->Psi;
|
---|
| 1579 | struct RunStruct *R = &P->R;
|
---|
| 1580 | struct LatticeLevel *LevS = R->LevS;
|
---|
| 1581 | struct LatticeLevel *Lev0 = R->Lev0;
|
---|
| 1582 | struct Density *Dens0 = Lev0->Dens;
|
---|
| 1583 | struct OnePsiElement *OnePsiA, *LOnePsiA;
|
---|
| 1584 | MPI_Status status;
|
---|
| 1585 | int ElementSize = (sizeof(fftw_complex) / sizeof(double));
|
---|
| 1586 | fftw_complex *LPsiDatA;
|
---|
| 1587 | int i, p, RecvSource;
|
---|
| 1588 | int Num = Psi->NoOfPsis;
|
---|
| 1589 |
|
---|
| 1590 | if (!P->Files.MeOutVis) return;
|
---|
| 1591 |
|
---|
| 1592 | P->Files.OutputPosType = (enum ModeType *) Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+(Num)),"OutputVis");
|
---|
| 1593 |
|
---|
| 1594 | P->Files.OutVisStep += offset;
|
---|
| 1595 | P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse;
|
---|
| 1596 | for (i=0; i < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; i++) { // go through all wave functions (plus the extra one for each process)
|
---|
| 1597 | OnePsiA = &Psi->AllPsiStatus[i]; // grab the desired OnePsiA
|
---|
| 1598 | if (OnePsiA->PsiType == type) { // drop if extra one
|
---|
| 1599 | if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // local?
|
---|
| 1600 | LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo];
|
---|
| 1601 | else
|
---|
| 1602 | LOnePsiA = NULL;
|
---|
| 1603 | if (LOnePsiA == NULL) { // if it's not local ... receive it from respective process into TempPsi
|
---|
| 1604 | RecvSource = OnePsiA->my_color_comm_ST_Psi;
|
---|
| 1605 | MPI_Recv( LevS->LPsi->TempPsi, LevS->MaxG*ElementSize, MPI_DOUBLE, RecvSource, WannierTag1, P->Par.comm_ST_PsiT, &status );
|
---|
| 1606 | LPsiDatA=LevS->LPsi->TempPsi;
|
---|
| 1607 | } else { // .. otherwise send it to all other processes (Max_me... - 1)
|
---|
| 1608 | for (p=0;p<P->Par.Max_me_comm_ST_PsiT;p++)
|
---|
| 1609 | if (p != OnePsiA->my_color_comm_ST_Psi)
|
---|
| 1610 | MPI_Send( LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo], LevS->MaxG*ElementSize, MPI_DOUBLE, p, WannierTag1, P->Par.comm_ST_PsiT);
|
---|
| 1611 | LPsiDatA=LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo];
|
---|
| 1612 | } // LPsiDatA is now set to the coefficients of OnePsi either stored or MPI_Received
|
---|
| 1613 |
|
---|
| 1614 | if (P->Files.MeOutVis) {
|
---|
| 1615 | P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse;
|
---|
| 1616 | // recalculate density for the specific wave function ...
|
---|
| 1617 | CalculateOneDensityR(Lat, LevS, Dens0, LPsiDatA, Dens0->DensityArray[ActualDensity], R->FactorDensityR, 0);
|
---|
| 1618 | // ... and output (wherein ActualDensity is used instead of TotalDensity)
|
---|
| 1619 | //OutputVis(P);
|
---|
| 1620 | CreateDensityOutputGeneral(P, P->Par.me_comm_ST_Psi);
|
---|
| 1621 | OutVisDensity(P, Dens0->DensityArray[ActualDensity]);
|
---|
| 1622 | OutVisIons(P);
|
---|
| 1623 | if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Written OutVisStep %i to disk\n", P->Par.me, P->Files.OutVisStep);
|
---|
| 1624 | P->Files.OutVisStep+=increment;
|
---|
| 1625 | P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse;
|
---|
| 1626 | }
|
---|
| 1627 | }
|
---|
| 1628 | }
|
---|
| 1629 | }
|
---|
| 1630 |
|
---|
| 1631 | /** Read source files containing stored calculations.
|
---|
| 1632 | * Calls ReadSrcPsiDensity() and ReadSrcIons().
|
---|
| 1633 | * \param *P Problem at hand
|
---|
| 1634 | */
|
---|
| 1635 | void ReadSrcFiles(struct Problem *P)
|
---|
| 1636 | {
|
---|
| 1637 | if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)ReadSrcPsiDensity\n", P->Par.me);
|
---|
| 1638 | ReadSrcPsiDensity(P, Occupied, 0, P->R.LevSNo);
|
---|
| 1639 | ReadSrcPsiDensity(P, UnOccupied, 0, P->R.LevSNo);
|
---|
| 1640 | if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)ReadSrcIons\n", P->Par.me);
|
---|
| 1641 | ReadSrcIons(P);
|
---|
| 1642 | }
|
---|
| 1643 |
|
---|
| 1644 | /** Plots a cut plane of the real density of one wave function.
|
---|
| 1645 | * \param *P Problem at hand
|
---|
| 1646 | * \param index index of axis (vector orthogonal to plane)
|
---|
| 1647 | * \param node node specifying where to cut at the given axis
|
---|
| 1648 | * \param wavenr global number of wave function
|
---|
| 1649 | * \param *density density array to plot
|
---|
| 1650 | * \sa PlotVectorPlane() - very similar
|
---|
| 1651 | */
|
---|
| 1652 | void PlotSrcPlane(struct Problem *P, int index, double node, int wavenr, fftw_real *density)
|
---|
| 1653 | {
|
---|
| 1654 | struct RunStruct *R = &P->R;
|
---|
| 1655 | struct Lattice *Lat = &P->Lat;
|
---|
| 1656 | struct LatticeLevel *Lev0 = R->Lev0;
|
---|
| 1657 | const int myPE = P->Par.me_comm_ST_Psi;
|
---|
| 1658 | char filename[255], spin[12];
|
---|
| 1659 | char *suchpointer;
|
---|
| 1660 | FILE *PlotFile = NULL;
|
---|
| 1661 | time_t seconds;
|
---|
| 1662 |
|
---|
| 1663 |
|
---|
| 1664 | // open file
|
---|
| 1665 | time(&seconds); // get current time
|
---|
| 1666 |
|
---|
| 1667 | switch (Lat->Psi.PsiST) {
|
---|
| 1668 | case SpinDouble:
|
---|
| 1669 | sprintf(&filename[0], ".psi%i_cut%i.csv", wavenr, index);
|
---|
| 1670 | strncat(spin,"SpinDouble",12);
|
---|
| 1671 | break;
|
---|
| 1672 | case SpinUp:
|
---|
| 1673 | sprintf(&filename[0], ".psi%i_cut%i_up.csv", wavenr, index);
|
---|
| 1674 | strncat(spin,"SpinUp",12);
|
---|
| 1675 | break;
|
---|
| 1676 | case SpinDown:
|
---|
| 1677 | sprintf(&filename[0], ".psi%i_cut%i_down.csv", wavenr, index);
|
---|
| 1678 | strncat(spin,"SpinDown",12);
|
---|
| 1679 | break;
|
---|
| 1680 | }
|
---|
| 1681 |
|
---|
| 1682 | if (!myPE) { // only process 0 writes to file
|
---|
| 1683 | OpenFile(P, &PlotFile, filename, "w", P->Call.out[ReadOut]);
|
---|
| 1684 | strcpy(filename, ctime(&seconds));
|
---|
| 1685 | suchpointer = strchr(filename, '\n');
|
---|
| 1686 | if (suchpointer != NULL)
|
---|
| 1687 | *suchpointer = '\0';
|
---|
| 1688 | if (PlotFile != NULL) {
|
---|
| 1689 | fprintf(PlotFile,"# Psi %i, type %s (real density) plot of plane perpendicular to direction e_%i at node %lg, seed %i, config %s, run on %s, #cpus %i", wavenr, spin, index, node, R->Seed, P->Files.default_path, filename, P->Par.Max_me_comm_ST_Psi);
|
---|
| 1690 | fprintf(PlotFile,"\n");
|
---|
| 1691 | } else { Error(SomeError, "PlotSrcPlane: Opening Plot File"); }
|
---|
| 1692 | }
|
---|
| 1693 |
|
---|
| 1694 | // plot density
|
---|
| 1695 | PlotRealDensity(P, Lev0, PlotFile, index, node, density, density);
|
---|
| 1696 |
|
---|
| 1697 | if (PlotFile != NULL) {
|
---|
| 1698 | // close file
|
---|
| 1699 | fclose(PlotFile);
|
---|
| 1700 | }
|
---|
| 1701 | }
|
---|
| 1702 |
|
---|
| 1703 | /** plots a cut plane of a given 3d real density.
|
---|
| 1704 | * \param *P Problem at hand, contains pointer to Lattice structure
|
---|
| 1705 | * \param *Lev LatticeLevel of the real density
|
---|
| 1706 | * \param PlotFile file pointer (already open and valid)
|
---|
| 1707 | * \param index index of lattice axis
|
---|
| 1708 | * \param n_orth position on lattice axis where to cut
|
---|
| 1709 | * \param *density1 first real density array
|
---|
| 1710 | * \param *density2 second real density array (point to \a *density1 if not needed)
|
---|
| 1711 | */
|
---|
| 1712 | void PlotRealDensity(struct Problem *P, struct LatticeLevel *Lev, FILE *PlotFile, int index, double n_orth, fftw_real *density1, fftw_real *density2)
|
---|
| 1713 | {
|
---|
| 1714 | struct Lattice *Lat = &P->Lat;
|
---|
| 1715 | int n[NDIM], n0;
|
---|
| 1716 | int N[NDIM];
|
---|
| 1717 | N[0] = Lev->Plan0.plan->N[0];
|
---|
| 1718 | N[1] = Lev->Plan0.plan->N[1];
|
---|
| 1719 | N[2] = Lev->Plan0.plan->N[2];
|
---|
| 1720 | const int N0 = Lev->Plan0.plan->local_nx;
|
---|
| 1721 | const int myPE = P->Par.me_comm_ST_Psi;
|
---|
| 1722 | double fac[NDIM], x[NDIM];
|
---|
| 1723 | int i0, i = 0;
|
---|
| 1724 | int PE, zahl;
|
---|
| 1725 | double *buffer;
|
---|
| 1726 | int node = ceil(n_orth/Lat->RealBasisQ[index]*N[index]);
|
---|
| 1727 | MPI_Status status;
|
---|
| 1728 | int sizes[P->Par.Max_me_comm_ST_Psi], c0, c1;
|
---|
| 1729 |
|
---|
| 1730 | switch (index) {
|
---|
| 1731 | case 0:
|
---|
| 1732 | zahl = 4*N[1]*N[2];
|
---|
| 1733 | break;
|
---|
| 1734 | case 1:
|
---|
| 1735 | zahl = 4*N0*N[2];
|
---|
| 1736 | break;
|
---|
| 1737 | case 2:
|
---|
| 1738 | zahl = 4*N0*N[1];
|
---|
| 1739 | break;
|
---|
| 1740 | }
|
---|
| 1741 | buffer = Malloc(sizeof(double)*zahl,"PlotRealDensity: buffer");
|
---|
| 1742 |
|
---|
| 1743 | c0 = cross(index,0);
|
---|
| 1744 | c1 = cross(index,1);
|
---|
| 1745 | // then for every point on the grid in real space ...
|
---|
| 1746 | for (n0=0;n0<N0;n0++) // only local points on x axis
|
---|
| 1747 | for (n[1]=0;n[1]<N[1];n[1]++)
|
---|
| 1748 | for (n[2]=0;n[2]<N[2];n[2]++) {
|
---|
| 1749 | n[0]=n0 + N0*myPE; // global relative coordinate: due to partitoning of x-axis in PEPGamma>1 case
|
---|
| 1750 | if (n[index] == node) { // only on the correct plane orthogonal to desired axis and at desired node ...
|
---|
| 1751 | fac[0] = (double)n[0]/(double)N[0];
|
---|
| 1752 | fac[1] = (double)n[1]/(double)N[1];
|
---|
| 1753 | fac[2] = (double)n[2]/(double)N[2];
|
---|
| 1754 | RMat33Vec3(x, Lat->RealBasis, fac); // relative coordinate times basis matrix gives absolute ones
|
---|
| 1755 | i0 = n[2]+N[2]*(n[1]+N[1]*n0); // index to local density array
|
---|
| 1756 |
|
---|
| 1757 | buffer[i++] = x[c0]; // fill buffer
|
---|
| 1758 | buffer[i++] = x[c1];
|
---|
| 1759 | buffer[i++] = density1[i0];
|
---|
| 1760 | buffer[i++] = density2[i0];
|
---|
| 1761 | if (i > zahl) Error(SomeError, "PlotRealDensity: buffer too small!");
|
---|
| 1762 | }
|
---|
| 1763 | }
|
---|
| 1764 | // exchange sizes of each buffer
|
---|
| 1765 | MPI_Allgather(&i, 1, MPI_INT, sizes, 1, MPI_INT, P->Par.comm_ST_Psi);
|
---|
| 1766 | if (myPE == 0) {
|
---|
| 1767 | for (PE=0; PE < P->Par.Max_me_comm_ST_Psi; PE++) {
|
---|
| 1768 | if (PE != 0) {
|
---|
| 1769 | // receive them
|
---|
| 1770 | if (MPI_Recv(buffer, sizes[PE], MPI_DOUBLE, PE, PlotRealDensityTag, P->Par.comm_ST_Psi, &status) != MPI_SUCCESS)
|
---|
| 1771 | Error(SomeError, "PlotRealDensity: MPI_Recv failure!");
|
---|
| 1772 | MPI_Get_count(&status, MPI_DOUBLE, &zahl);
|
---|
| 1773 | if (zahl != sizes[PE])
|
---|
| 1774 | Error(SomeError, "PlotRealDensity: received unexpected amount of elements!");
|
---|
| 1775 | }
|
---|
| 1776 | //write them: local one (still in buffer) and received ones
|
---|
| 1777 | for (i0 = 0; i0 < sizes[PE];) {
|
---|
| 1778 | fprintf(PlotFile,"%e", buffer[(i0)++]);
|
---|
| 1779 | if ((i0 % 4) == 0) {
|
---|
| 1780 | fprintf(PlotFile,"\n");
|
---|
| 1781 | } else {
|
---|
| 1782 | fprintf(PlotFile,"\t");
|
---|
| 1783 | }
|
---|
| 1784 | }
|
---|
| 1785 | }
|
---|
| 1786 | } else { // send them
|
---|
| 1787 | if (MPI_Send(buffer, i, MPI_DOUBLE, 0, PlotRealDensityTag, P->Par.comm_ST_Psi) != MPI_SUCCESS)
|
---|
| 1788 | Error(SomeError, "PlotRealDensity: MPI_Send failure!");
|
---|
| 1789 | }
|
---|
[64fa9e] | 1790 | Free(buffer, "PlotRealDensity: buffer");
|
---|
[a0bcf1] | 1791 | }
|
---|