IR_PrintVisualizationSWE.scala 13.6 KB
Newer Older
1
2
package exastencils.applications.swe.ir

3
import scala.collection.immutable.ListMap
4
5
import scala.collection.mutable.ListBuffer

6
import exastencils.base.ir._
7
import exastencils.base.ir.IR_ImplicitConversion._
8
import exastencils.baseExt.ir._
9
10
import exastencils.config.Knowledge
import exastencils.core.Duplicate
11
import exastencils.domain.ir.IR_IV_IsValidForDomain
12
import exastencils.field.ir._
13
import exastencils.grid.ir.IR_AtNode
14
import exastencils.grid.ir.IR_VF_NodePositionAsVec
15
import exastencils.io.ir._
16
17
import exastencils.logger.Logger
import exastencils.util.ir.IR_Print
18
import exastencils.visualization.ir.postprocessing.IR_PrintVisualizationTriangles
19

20
trait IR_PrintVisualizationSWE extends IR_PrintVisualizationTriangles {
21
22
  def numDimsGrid = 2

23
24
  def numCells_x : Int = etaDiscLower0.layout.layoutsPerDim(0).numInnerLayers
  def numCells_y : Int = etaDiscLower0.layout.layoutsPerDim(1).numInnerLayers
25
  def numCells_z = 1
26
  def numCellsPerFrag : IR_Expression = 2 * numCells_x * numCells_y * numCells_z
27

28
  override def dimsPositionsFrag : ListBuffer[IR_Expression] = if (Knowledge.swe_nodalReductionPrint) ListBuffer(numCells_x + 1, numCells_y + 1) else ListBuffer(6, numCells_x, numCells_y)
29

30
31
32
  // passed as arguments to function
  def discFieldCollection : ListBuffer[ListBuffer[IR_Field]]
  def nodalFieldCollection : ListBuffer[IR_Field]
33

34
35
  def nodalFields : ListMap[String, IR_Field] = ListMap(
    nodalFieldCollection.map(field => field.name -> field) : _*)
36
  def discFields : ListMap[String, ListBuffer[IR_Field]] = ListMap(
37
    discFieldCollection.map(discField => getBasenameDiscField(discField) -> discField) : _*)
38
  def discFieldsReduced : ListMap[String, IR_IV_TemporaryBuffer] = discFields.map { discField =>
39
    discField._1 -> IR_IV_TemporaryBuffer(discField._2.head.resolveBaseDatatype, IR_AtNode, discField._1, domainIndex, blockwise = true, dimsPositionsFrag)
40
  }
41

42
  def etaDiscLower0 : IR_Field = IR_FieldCollection.getByIdentifier("etaDiscLower0", level).get
43
  def someCellField : IR_Field = etaDiscLower0
44

45
46
  def fields : ListMap[String, ListBuffer[IR_Field]] = nodalFields.map(field => field._1 -> ListBuffer(field._2)) ++ discFields
  def fieldnames : ListBuffer[String] = fields.keys.to[ListBuffer]
47
  def numFields : Int = fieldnames.length
48

49
  def nodePosVecAsDataBuffers(accessIndices : Option[ListBuffer[IR_Index]], datasets : Option[ListBuffer[IR_Expression]]) : ListBuffer[IR_DataBuffer] = {
50
    (0 until numDimsGrid).map(dim =>
51
      IR_DataBuffer(IR_VF_NodePositionAsVec.find(level).associatedField, accessIndices, if (datasets.isDefined) Some(datasets.get(dim)) else None, dim, 0, canonicalFileLayout)
52
    ).to[ListBuffer]
53
  }
54

55
56
  // get the common prefix of a disc field and use as name (e.g. etaDiscLower0, etaDiscLower1, ... -> etaDisc)
  def getBasenameDiscField(discField : ListBuffer[IR_Field]) : String = {
57
58
59
60
61
62
    val basename = discField.map(_.name).reduce((a, b) => (a zip b).takeWhile(Function.tupled(_ == _)).map(_._1).mkString)
    if (basename.isEmpty) {
      Logger.error("\"IR_PrintVisualizationSWE:\" Could not extract a common name from disc field components. Components do not belong to the same disc field.")
    }

    basename
63
  }
64

65
66
67
68
69
70
71
72
  /*
  Store disc components interleaved in buffer
    - in exodus, the individual datasets of each disc field component cannot be joined together as in xdmf
      - this causes the disc components to be written in an interleaved fashion
      - to reduce the number of interleaved accesses, the disc components are joined together in a buffer before writing
    - can be enabled for comparison in Xdmf printers via "swe_interleaveDiscComponentsPrint" flag
  */
  def discFieldBuffers : ListMap[String, IR_IV_TemporaryBuffer] = discFields.map { discField =>
73
    discField._1 -> IR_IV_TemporaryBuffer(discField._2.head.resolveBaseDatatype, IR_AtNode, discField._1, domainIndex, blockwise = true, dimsPositionsFrag)
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
  }

  def setupNonReducedDiscData() : ListBuffer[IR_Statement] = {
    var stmts : ListBuffer[IR_Statement] = ListBuffer()

    // join disc components
    discFieldBuffers.values.foreach { tmpBuf =>
      val indexRangeCells = IR_ExpressionIndexRange(
        IR_ExpressionIndex((0 until numDimsGrid).toArray.map(dim => etaDiscLower0.layout.idxById("IB", dim) - Duplicate(etaDiscLower0.referenceOffset(dim)) : IR_Expression)),
        IR_ExpressionIndex((0 until numDimsGrid).toArray.map(dim => etaDiscLower0.layout.idxById("IE", dim) - Duplicate(etaDiscLower0.referenceOffset(dim)) : IR_Expression)))

      val numAccessesPerCell = nodeOffsets.length // for non-reduced SWE "6"
      val offset = IR_LoopOverFragments.defIt * dimsPositionsFrag.reduce(_ * _) + numAccessesPerCell * indexRangeCells.linearizeIndex(IR_LoopOverDimensions.defIt(numDimsGrid))

      stmts += tmpBuf.allocateMemory
      stmts += IR_LoopOverFragments(
        IR_IfCondition(IR_IV_IsValidForDomain(domainIndex),
          IR_LoopOverDimensions(numDimsGrid, indexRangeCells,
            (0 until numAccessesPerCell).to[ListBuffer].map(idx => {
              IR_Assignment(
Richard Angersbach's avatar
Richard Angersbach committed
94
                IR_ArrayAccess(tmpBuf, offset + idx),
95
                IR_FieldAccess(discFields(tmpBuf.name)(idx), IR_IV_ActiveSlot(someCellField), IR_LoopOverDimensions.defIt(numDimsGrid))) : IR_Statement
96
97
98
99
100
101
            }))))
    }

    stmts
  }

102
  // access pattern dependent on reduction mode for blockstructured meshes
103
  def accessIndices : Option[ListBuffer[IR_Index]] = if (Knowledge.swe_nodalReductionPrint)
104
105
106
107
    None
  else
    Some(nodeOffsets.map(_.toExpressionIndex))

108
  def nodalAccess(field : IR_Field) = if (Knowledge.swe_nodalReductionPrint)
109
    IR_RegularAccessPattern(IR_AccessFieldFunction(field, IR_IV_ActiveSlot(field)))
110
  else
111
    IR_SWEAccessPattern(IR_AccessFieldFunction(field, IR_IV_ActiveSlot(field)), accessIndices)
112

113
  // glue logic for disc fields to be mapped to data buffers
114
  def discFieldsToDatabuffers(discField : ListBuffer[IR_Field]) : ListBuffer[IR_DataBuffer] = ???
115

116
  // calculating an average over a disc field to reduce the amount of data written to file
117
  def setupReducedData : ListBuffer[IR_Statement] = {
118
119
    var stmts : ListBuffer[IR_Statement] = ListBuffer()

120
121
122
123
    // allocate buffers before calculation
    discFieldsReduced.values.foreach { tmpBuf =>
      stmts += tmpBuf.allocateMemory
    }
124
125
126

    // compute averages and store into buffer before writing
    var stmtsFragLoop : ListBuffer[IR_Statement] = ListBuffer()
127
    discFieldsReduced.values.foreach { tmpBuf =>
128
      stmtsFragLoop += IR_Comment(s"\n Nodal data reduction for ${ tmpBuf.name } \n")
129
      stmtsFragLoop ++= reducedCellPrint(IR_VariableAccess(tmpBuf.name, tmpBuf.resolveDatatype()), discFields(tmpBuf.name))
130
131
132
    }

    stmts += IR_LoopOverFragments(
133
      IR_IfCondition(IR_IV_IsValidForDomain(domainIndex),
134
135
136
137
138
139
140
        stmtsFragLoop
      )
    )

    stmts
  }

141
  // nodal data reduction
142
  def reducedCellPrint(buf : IR_VariableAccess, discField : ListBuffer[IR_Field], indentation : Option[IR_StringConstant] = None) : ListBuffer[IR_Statement] = {
143

144
145
    val low : ListBuffer[IR_Field] = discField.take(3)
    val upp : ListBuffer[IR_Field] = discField.takeRight(3)
146

147
    if (discField.length != 6) {
148
      Logger.error("Wrong usage of \"addReducedNodePrint\" in IR_PrintVisualizationSWE.")
149
150
    }

151
    def getIdxNodalLoop(idx : IR_ExpressionIndex) : IR_Expression = IR_ExpressionIndexRange(
152
153
154
155
      IR_ExpressionIndex((0 until numDimsGrid).toArray.map(dim => etaDiscLower0.layout.idxById("IB", dim) - Duplicate(etaDiscLower0.referenceOffset(dim)) : IR_Expression)),
      IR_ExpressionIndex((0 until numDimsGrid).toArray.map(dim => 1 + etaDiscLower0.layout.idxById("IE", dim) - Duplicate(etaDiscLower0.referenceOffset(dim)) : IR_Expression))
    ).linearizeIndex(idx)

156
    def storeOperation(toStore : IR_Expression, idx : IR_ExpressionIndex) : IR_Statement = buf.datatype match {
157
      case IR_SpecialDatatype("std::ofstream") => IR_Print(buf,
158
        (if (indentation.isDefined) indentation.get :: Nil else Nil) :+ toStore :+ IR_Print.newline : _*)
159
      case IR_PointerDatatype(_)               => IR_Assignment(IR_ArrayAccess(buf, numPointsPerFrag * IR_LoopOverFragments.defIt + getIdxNodalLoop(idx)), toStore)
160
    }
161
162
163
164
165
166
167
168
169

    def constructLoopForDim(dim : Int, offStart : Int, offEnd : Int, loopStmts : IR_Statement*) = IR_ForLoop(
      IR_VariableDeclaration(IR_LoopOverDimensions.defItForDim(dim), etaDiscLower0.layout.idxById("IB", dim) + offStart - Duplicate(etaDiscLower0.referenceOffset(dim))),
      IR_LoopOverDimensions.defItForDim(dim) < etaDiscLower0.layout.idxById("IE", dim) - offEnd - Duplicate(etaDiscLower0.referenceOffset(dim)),
      IR_PreIncrement(IR_LoopOverDimensions.defItForDim(dim)),
      loopStmts.to[ListBuffer])

    var statementsFragLoop : ListBuffer[IR_Statement] = ListBuffer()

170
171
    // expression indices for the edges
    val baseIdx_edgeLower = IR_ExpressionIndex(IR_LoopOverDimensions.defItForDim(0), 0)
172
    val baseIdx_edgeLeft = IR_ExpressionIndex(0, IR_LoopOverDimensions.defItForDim(1))
173
174
175
176
177
178
179
180
181
    val baseIdx_edgeRight = IR_ExpressionIndex(someCellField.layout.layoutsPerDim(0).numInnerLayers, IR_LoopOverDimensions.defItForDim(1))
    val baseIdx_edgeUpper = IR_ExpressionIndex(IR_LoopOverDimensions.defItForDim(0), someCellField.layout.layoutsPerDim(1).numInnerLayers)

    // expression indices for the corners
    val baseIdx_cornerLL = IR_ExpressionIndex(0, 0)
    val baseIdx_cornerLR = IR_ExpressionIndex(low(1).layout.layoutsPerDim(0).numInnerLayers, 0)
    val baseIdx_cornerUL = IR_ExpressionIndex(0, low(2).layout.layoutsPerDim(1).numInnerLayers)
    val baseIdx_cornerUR = IR_ExpressionIndex((0 until numDimsGrid).map(d => someCellField.layout.layoutsPerDim(d).numInnerLayers).toArray)

182
    statementsFragLoop += IR_Comment("lower left corner: vertex contained by 1 triangle")
183
184
185
    statementsFragLoop += storeOperation(
      IR_FieldAccess(low(0), IR_IV_ActiveSlot(low(0)), baseIdx_cornerLL),
      baseIdx_cornerLL)
186
187
188

    statementsFragLoop += IR_Comment("lowermost row w/o corners: vertex contained by 3 triangles")
    statementsFragLoop += constructLoopForDim(dim = 0, offStart = 1, offEnd = 0,
189
190
191
192
193
      storeOperation(
        (IR_FieldAccess(low(0), IR_IV_ActiveSlot(low(0)), baseIdx_edgeLower) +
          IR_FieldAccess(low(1), IR_IV_ActiveSlot(low(1)), baseIdx_edgeLower + IR_ConstIndex(-1, 0)) +
          IR_FieldAccess(upp(2), IR_IV_ActiveSlot(upp(2)), baseIdx_edgeLower + IR_ConstIndex(-1, 0))) / 3.0,
        baseIdx_edgeLower))
194
195

    statementsFragLoop += IR_Comment("lower right corner: vertex contained by 2 triangles")
196
197
198
199
    statementsFragLoop += storeOperation(
      (IR_FieldAccess(low(1), IR_IV_ActiveSlot(low(1)), baseIdx_cornerLR + IR_ConstIndex(-1, 0)) +
        IR_FieldAccess(upp(2), IR_IV_ActiveSlot(upp(2)), baseIdx_cornerLR + IR_ConstIndex(-1, 0))) / 2.0,
      baseIdx_cornerLR)
200
201
202
203

    statementsFragLoop += IR_Comment("inner rows")
    statementsFragLoop += constructLoopForDim(dim = 1, offStart = 1, offEnd = 0,
      IR_Comment("leftmost column: vertex contained by 3 triangles"),
204
205
206
207
208
      storeOperation(
        (IR_FieldAccess(low(0), IR_IV_ActiveSlot(low(0)), baseIdx_edgeLeft) +
          IR_FieldAccess(low(2), IR_IV_ActiveSlot(low(2)), baseIdx_edgeLeft + IR_ConstIndex(0, -1)) +
          IR_FieldAccess(upp(1), IR_IV_ActiveSlot(upp(1)), baseIdx_edgeLeft + IR_ConstIndex(0, -1))) / 3.0,
        baseIdx_edgeLeft),
209
210
      IR_Comment("inner points: vertex contained by 6 triangles"),
      constructLoopForDim(dim = 0, offStart = 1, offEnd = 0,
211
        storeOperation(
212
213
214
215
216
217
          (IR_FieldAccess(low(0), IR_IV_ActiveSlot(low(0)), IR_LoopOverDimensions.defIt(numDimsGrid)) +
            IR_FieldAccess(low(1), IR_IV_ActiveSlot(low(1)), IR_LoopOverDimensions.defIt(numDimsGrid) + IR_ConstIndex(-1, 0)) +
            IR_FieldAccess(upp(2), IR_IV_ActiveSlot(upp(2)), IR_LoopOverDimensions.defIt(numDimsGrid) + IR_ConstIndex(-1, 0)) +
            IR_FieldAccess(low(2), IR_IV_ActiveSlot(low(2)), IR_LoopOverDimensions.defIt(numDimsGrid) + IR_ConstIndex(0, -1)) +
            IR_FieldAccess(upp(1), IR_IV_ActiveSlot(upp(1)), IR_LoopOverDimensions.defIt(numDimsGrid) + IR_ConstIndex(0, -1)) +
            IR_FieldAccess(upp(0), IR_IV_ActiveSlot(upp(0)), IR_LoopOverDimensions.defIt(numDimsGrid) + IR_ConstIndex(-1, -1))) / 6.0,
218
          IR_LoopOverDimensions.defIt(numDimsGrid))),
219
      IR_Comment("rightmost column: vertex contained by 3 triangles"),
220
221
222
223
224
      storeOperation(
        (IR_FieldAccess(upp(0), IR_IV_ActiveSlot(upp(0)), baseIdx_edgeRight + IR_ConstIndex(-1, -1)) +
          IR_FieldAccess(low(1), IR_IV_ActiveSlot(low(1)), baseIdx_edgeRight + IR_ConstIndex(-1, 0)) +
          IR_FieldAccess(upp(2), IR_IV_ActiveSlot(upp(2)), baseIdx_edgeRight + IR_ConstIndex(-1, 0))) / 3.0,
        baseIdx_edgeRight)
225
226
227
    )

    statementsFragLoop += IR_Comment("upper left corner: vertex contained by 2 triangles")
228
229
230
231
    statementsFragLoop += storeOperation(
      (IR_FieldAccess(low(2), IR_IV_ActiveSlot(low(2)), baseIdx_cornerUL + IR_ConstIndex(0, -1)) +
        IR_FieldAccess(upp(1), IR_IV_ActiveSlot(upp(1)), baseIdx_cornerUL + IR_ConstIndex(0, -1))) / 2.0,
      baseIdx_cornerUL)
232
233
234

    statementsFragLoop += IR_Comment("uppermost row w/o corners: vertex contained by 3 triangles")
    statementsFragLoop += constructLoopForDim(dim = 0, offStart = 1, offEnd = 0,
235
236
237
238
239
      storeOperation(
        (IR_FieldAccess(upp(0), IR_IV_ActiveSlot(upp(0)), baseIdx_edgeUpper + IR_ConstIndex(-1, -1)) +
          IR_FieldAccess(low(2), IR_IV_ActiveSlot(low(2)), baseIdx_edgeUpper + IR_ConstIndex(0, -1)) +
          IR_FieldAccess(upp(1), IR_IV_ActiveSlot(upp(1)), baseIdx_edgeUpper + IR_ConstIndex(0, -1))) / 3.0,
        baseIdx_edgeUpper))
240
241

    statementsFragLoop += IR_Comment("upper right corner: vertex contained by 1 triangle")
242
243
244
    statementsFragLoop += storeOperation(
      IR_FieldAccess(upp(0), IR_IV_ActiveSlot(upp(0)), baseIdx_cornerUR + IR_ConstIndex(-1, -1)),
      baseIdx_cornerUR)
245
246
247
248

    statementsFragLoop
  }
}