m2c/MaterialVolumeOutput.cpp at main · kevinwgy/m2c · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
/************************************************************************
 * Copyright © 2020 The Multiphysics Modeling and Computation (M2C) Lab
 * <kevin.wgy@gmail.com> <kevinw3@vt.edu>
 ************************************************************************/

#include <MaterialVolumeOutput.h>
#include <IoData.h>

//-------------------------------------------------------------------------

MaterialVolumeOutput::MaterialVolumeOutput(MPI_Comm &comm_, IoData &iod,
                                           SpaceVariable3D& cell_volume_) :
                      comm(comm_), cell_volume(cell_volume_)
{
  frequency = iod.output.materialVolumes.frequency;
  frequency_dt = iod.output.materialVolumes.frequency_dt;

  last_snapshot_time = -1.0;

  numMaterials = iod.eqs.materials.dataMap.size() + 1; //an extra one for "ghost/inactive"

  // open file (all processors will open the same file, but only proc 0 will write to it)
  int spn = strlen(iod.output.prefix) + 1;
  if(iod.output.materialVolumes.filename[0] != 0) {
    char *full_fname = new char[spn + strlen(iod.output.materialVolumes.filename)];
    snprintf(full_fname, spn + strlen(iod.output.materialVolumes.filename),
             "%s%s", iod.output.prefix, iod.output.materialVolumes.filename);
    file = fopen(full_fname, "w");

    if(!file) {
      print_error("*** Error: Cannot open file '%s' for output.\n", full_fname);
      exit_mpi();
    }

    print(file, "## Number of physical materials: %d (0 - %d); ID for ghost/inactive cells: %d.\n",
                numMaterials-1, numMaterials-2, numMaterials-1);
    print(file, "## Total volume of computational domain: %e.\n",
                (iod.mesh.xmax-iod.mesh.x0)*(iod.mesh.ymax-iod.mesh.y0)*
                (iod.mesh.zmax-iod.mesh.z0));
    print(file, "## The last column is the sum over all material subdomains "
                "(including ghost/inactive).\n");
    mpi_barrier();
    fflush(file);
    delete [] full_fname;
  }
  else
    file = NULL;

}

//-------------------------------------------------------------------------

MaterialVolumeOutput::~MaterialVolumeOutput()
{
  if(file)
    fclose(file);
}

//-------------------------------------------------------------------------

void
MaterialVolumeOutput::WriteSolution(double time, double dt, int time_step, SpaceVariable3D& ID,
                                    bool force_write)
{

  if(file == NULL)
    return;

  if(!isTimeToWrite(time,dt,time_step,frequency_dt,frequency,last_snapshot_time,force_write))
    return;

  double vol[numMaterials];
  ComputeMaterialVolumes(ID, vol);

  print(file, "%8d  %16.8e  ", time_step, time);
  double sum = 0.0;
  for(int i=0; i<numMaterials; i++) {
    print(file, "%16.8e  ", vol[i]);
    sum += vol[i];
  }
  print(file, "%16.8e\n", sum);
  mpi_barrier();
  fflush(file);

  last_snapshot_time = time;

}

//-------------------------------------------------------------------------

void
MaterialVolumeOutput::ComputeMaterialVolumes(SpaceVariable3D& ID, double* vol)
{

  double*** id   = (double***)ID.GetDataPointer();
  double*** cell = (double***)cell_volume.GetDataPointer();

  int i0, j0, k0, imax, jmax, kmax;
  ID.GetCornerIndices(&i0, &j0, &k0, &imax, &jmax, &kmax);

  //initialize volumes to 0
  for(int i=0; i<numMaterials; i++)
    vol[i] = 0.0;

  // loop through the true domain (excluding the ghost layer)
  int myid;
  for(int k=k0; k<kmax; k++) {
    for(int j=j0; j<jmax; j++) {
      for(int i=i0; i<imax; i++) {
        myid = id[k][j][i];
        if(myid<0 || myid>=numMaterials) {
          fprintf(stdout,"*** Error: Detected an unrecognized material id (%d)\n",
                  myid);
          exit(-1);
        }
        vol[myid] += cell[k][j][i];
      }
    }
  }

  //MPI communication
  MPI_Allreduce(MPI_IN_PLACE, vol, numMaterials, MPI_DOUBLE, MPI_SUM, comm);

  ID.RestoreDataPointerToLocalVector();
  cell_volume.RestoreDataPointerToLocalVector();

}

//-------------------------------------------------------------------------