Commit 4274526c authored by Jean-Matthieu Gallard's avatar Jean-Matthieu Gallard
Browse files

KernelGen - Use Python's Decimal instead of float in MathUtils, DGMatrices +...

KernelGen - Use Python's Decimal instead of float in MathUtils, DGMatrices + Quadrature generated with 40 decimal places instead of 15
parent b43f016f
......@@ -25,12 +25,16 @@
from .abstractModelBaseClass import AbstractModelBaseClass
from ..utils import MathsUtils #matrix operation and build functions
import time #for runtime measurements debug print if context["runtimeDebug"]
import time #for runtime measurements debug print if context["runtimeDebug"]
from decimal import Decimal
class DGMatrixModel(AbstractModelBaseClass):
def generateCode(self):
MathsUtils.setDecimalPrecision() # set Python's decimal precision to mathUtils value (default 50)
if self.context["quadratureType"] == "Gauss-Legendre":
#get GL nodes
wGPN, xGPN = MathsUtils.getGaussLegendre(self.context["nDof"])
......@@ -40,7 +44,6 @@ class DGMatrixModel(AbstractModelBaseClass):
xGPN = xGPN[::-1]
else:
raise ValueError("Quadrature type "+self.context["quadratureType"]+" not supported." )
padSize = self.context["nDofPad"] - self.context["nDof"]
self.context["nDofPad_seq"] = range(self.context["nDofPad"])
......@@ -49,8 +52,8 @@ class DGMatrixModel(AbstractModelBaseClass):
start = time.perf_counter()
# [FLCoeff 0...0]; [FRCoeff 0...0];
# right now FLCoeff, FRCoeff no pad (gives no benefit w.r.t libxsmm)
FLCoeff, _ = MathsUtils.baseFunc1d(0.0, xGPN, self.context["nDof"]) #is also F0
FRCoeff, _ = MathsUtils.baseFunc1d(1.0, xGPN, self.context["nDof"])
FLCoeff, _ = MathsUtils.baseFunc1d(Decimal("0.0"), xGPN, self.context["nDof"]) #is also F0
FRCoeff, _ = MathsUtils.baseFunc1d(Decimal("1.0"), xGPN, self.context["nDof"])
self.context["FLCoeff"] = MathsUtils.vectorPad(FLCoeff, padSize)
self.context["FRCoeff"] = MathsUtils.vectorPad(FRCoeff, padSize)
if self.context["runtimeDebug"]:
......@@ -96,9 +99,9 @@ class DGMatrixModel(AbstractModelBaseClass):
#fineGridProjector1d_T_weighted
for i in range(self.context["nDof"]):
for j in range(self.context["nDof"]):
fineGridProjector1d_0[i][j] *= wGPN[j]/wGPN[i]/3.0
fineGridProjector1d_1[i][j] *= wGPN[j]/wGPN[i]/3.0
fineGridProjector1d_2[i][j] *= wGPN[j]/wGPN[i]/3.0
fineGridProjector1d_0[i][j] *= wGPN[j]/wGPN[i]/Decimal("3.0")
fineGridProjector1d_1[i][j] *= wGPN[j]/wGPN[i]/Decimal("3.0")
fineGridProjector1d_2[i][j] *= wGPN[j]/wGPN[i]/Decimal("3.0")
self.context["fineGridProjector1d_T_weighted_0"] = MathsUtils.matrixPadAndFlatten_RowMajor(fineGridProjector1d_0,padSize)
self.context["fineGridProjector1d_T_weighted_1"] = MathsUtils.matrixPadAndFlatten_RowMajor(fineGridProjector1d_1,padSize)
self.context["fineGridProjector1d_T_weighted_2"] = MathsUtils.matrixPadAndFlatten_RowMajor(fineGridProjector1d_2,padSize)
......
......@@ -26,11 +26,15 @@ from .abstractModelBaseClass import AbstractModelBaseClass
from ..utils import MathsUtils #matrix operation and build functions
from decimal import Decimal
class QuadratureModel(AbstractModelBaseClass):
def generateCode(self):
MathsUtils.setDecimalPrecision() # set Python's decimal precision to mathUtils value (default 50)
if self.context["quadratureType"] == "Gauss-Legendre":
QuadratureWeights, QuadratureNodes = MathsUtils.getGaussLegendre(self.context["nDof"])
OtherQuadratureWeights, OtherQuadratureNodes = MathsUtils.getGaussLobatto(self.context["nDof"])
......@@ -47,7 +51,7 @@ class QuadratureModel(AbstractModelBaseClass):
self.context["w1_seq"] = range(self.context["w1Size"])
#inverse of weights1
self.context["iweights1"] = [1.0/self.context["weights1"][i] if self.context["weights1"][i] != 0. else 0. for i in self.context["w1_seq"]]
self.context["iweights1"] = [Decimal("1.0")/self.context["weights1"][i] if self.context["weights1"][i] != 0. else 0. for i in self.context["w1_seq"]]
if(self.context["nDim"] == 2):
# weightsVector is QuadratureWeights itself
......@@ -95,7 +99,7 @@ class QuadratureModel(AbstractModelBaseClass):
print("QuadratureModel: nDim not supported")
# inverse of weight3
self.context["iweights3"] = [1.0/self.context["weights3"][i] if self.context["weights3"][i] != 0. else 0. for i in self.context["w3_seq"]]
self.context["iweights3"] = [Decimal("1.0")/self.context["weights3"][i] if self.context["weights3"][i] != 0. else 0. for i in self.context["w3_seq"]]
self.context["QuadratureWeights"], self.context["QuadratureNodes"] = QuadratureWeights, QuadratureNodes
self.context["quadrature_seq"] = range(self.context["nDof"])
......
......@@ -10,6 +10,8 @@
* Released under the BSD 3 Open Source License.
* For the full license text, see LICENSE.txt
**/ #}
{# format decimal numbers with 40 decimal places, exception for 0. in the padding to avoid getting 0.00...0e+40 from Python's Decimal #}
{% macro formatDecimal(a)%}{{"{:.40e}".format(a) if a != 0.0 else "0.0"}}{% endmacro %}
#include <mm_malloc.h> //g++
......@@ -100,51 +102,51 @@ void {{codeNamespace}}::initDGMatrices() {
}
{% for i in nDofPad_seq %}
FLCoeff[{{i}}] = {{"{:.15e}".format(FLCoeff[i])}};
FLCoeff[{{i}}] = {{formatDecimal(FLCoeff[i])}};
{% endfor %}
{% for i in nDofPad_seq %}
FRCoeff[{{i}}] = {{"{:.15e}".format(FRCoeff[i])}};
FRCoeff[{{i}}] = {{formatDecimal(FRCoeff[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
dudx[{{i}}] = {{"{:.15e}".format(dudx[i])}};
dudx[{{i}}] = {{formatDecimal(dudx[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
dudx_T[{{i}}] = {{"{:.15e}".format(dudx_T[i])}};
dudx_T[{{i}}] = {{formatDecimal(dudx_T[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
iK1_T[{{i}}] = {{"{:.15e}".format(iK1_T[i])}};
iK1_T[{{i}}] = {{formatDecimal(iK1_T[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
Kxi[{{i}}] = {{"{:.15e}".format(Kxi[i])}};
Kxi[{{i}}] = {{formatDecimal(Kxi[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
Kxi_T[{{i}}] = {{"{:.15e}".format(Kxi_T[i])}};
Kxi_T[{{i}}] = {{formatDecimal(Kxi_T[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
fineGridProjector1d[0][{{i}}] = {{"{:.15e}".format(fineGridProjector1d_0[i])}};
fineGridProjector1d[0][{{i}}] = {{formatDecimal(fineGridProjector1d_0[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
fineGridProjector1d[1][{{i}}] = {{"{:.15e}".format(fineGridProjector1d_1[i])}};
fineGridProjector1d[1][{{i}}] = {{formatDecimal(fineGridProjector1d_1[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
fineGridProjector1d[2][{{i}}] = {{"{:.15e}".format(fineGridProjector1d_2[i])}};
fineGridProjector1d[2][{{i}}] = {{formatDecimal(fineGridProjector1d_2[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
fineGridProjector1d_T_weighted[0][{{i}}] = {{"{:.15e}".format(fineGridProjector1d_T_weighted_0[i])}};
fineGridProjector1d_T_weighted[0][{{i}}] = {{formatDecimal(fineGridProjector1d_T_weighted_0[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
fineGridProjector1d_T_weighted[1][{{i}}] = {{"{:.15e}".format(fineGridProjector1d_T_weighted_1[i])}};
fineGridProjector1d_T_weighted[1][{{i}}] = {{formatDecimal(fineGridProjector1d_T_weighted_1[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
fineGridProjector1d_T_weighted[2][{{i}}] = {{"{:.15e}".format(fineGridProjector1d_T_weighted_2[i])}};
fineGridProjector1d_T_weighted[2][{{i}}] = {{formatDecimal(fineGridProjector1d_T_weighted_2[i])}};
{% endfor %}
......@@ -167,51 +169,51 @@ void {{codeNamespace}}::initDGMatrices() {
}
{% for i in nDofPad_seq %}
FLCoeff_SP[{{i}}] = {{"{:.15e}".format(FLCoeff[i])}};
FLCoeff_SP[{{i}}] = {{formatDecimal(FLCoeff[i])}};
{% endfor %}
{% for i in nDofPad_seq %}
FRCoeff_SP[{{i}}] = {{"{:.15e}".format(FRCoeff[i])}};
FRCoeff_SP[{{i}}] = {{formatDecimal(FRCoeff[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
dudx_SP[{{i}}] = {{"{:.15e}".format(dudx[i])}};
dudx_SP[{{i}}] = {{formatDecimal(dudx[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
dudx_T_SP[{{i}}] = {{"{:.15e}".format(dudx_T[i])}};
dudx_T_SP[{{i}}] = {{formatDecimal(dudx_T[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
iK1_T_SP[{{i}}] = {{"{:.15e}".format(iK1_T[i])}};
iK1_T_SP[{{i}}] = {{formatDecimal(iK1_T[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
Kxi_SP[{{i}}] = {{"{:.15e}".format(Kxi[i])}};
Kxi_SP[{{i}}] = {{formatDecimal(Kxi[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
Kxi_T_SP[{{i}}] = {{"{:.15e}".format(Kxi_T[i])}};
Kxi_T_SP[{{i}}] = {{formatDecimal(Kxi_T[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
fineGridProjector1d_SP[0][{{i}}] = {{"{:.15e}".format(fineGridProjector1d_0[i])}};
fineGridProjector1d_SP[0][{{i}}] = {{formatDecimal(fineGridProjector1d_0[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
fineGridProjector1d_SP[1][{{i}}] = {{"{:.15e}".format(fineGridProjector1d_1[i])}};
fineGridProjector1d_SP[1][{{i}}] = {{formatDecimal(fineGridProjector1d_1[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
fineGridProjector1d_SP[2][{{i}}] = {{"{:.15e}".format(fineGridProjector1d_2[i])}};
fineGridProjector1d_SP[2][{{i}}] = {{formatDecimal(fineGridProjector1d_2[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
fineGridProjector1d_T_weighted_SP[0][{{i}}] = {{"{:.15e}".format(fineGridProjector1d_T_weighted_0[i])}};
fineGridProjector1d_T_weighted_SP[0][{{i}}] = {{formatDecimal(fineGridProjector1d_T_weighted_0[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
fineGridProjector1d_T_weighted_SP[1][{{i}}] = {{"{:.15e}".format(fineGridProjector1d_T_weighted_1[i])}};
fineGridProjector1d_T_weighted_SP[1][{{i}}] = {{formatDecimal(fineGridProjector1d_T_weighted_1[i])}};
{% endfor %}
{% for i in nDofPadTimesnDof_seq %}
fineGridProjector1d_T_weighted_SP[2][{{i}}] = {{"{:.15e}".format(fineGridProjector1d_T_weighted_2[i])}};
fineGridProjector1d_T_weighted_SP[2][{{i}}] = {{formatDecimal(fineGridProjector1d_T_weighted_2[i])}};
{% endfor %}
{% endif %}{# useSinglePrecision #}
......
......@@ -10,6 +10,8 @@
* Released under the BSD 3 Open Source License.
* For the full license text, see LICENSE.txt
**/ #}
{# format decimal numbers with 40 decimal places, exception for 0. in the padding to avoid getting 0.00...0e+40 from Python's Decimal #}
{% macro formatDecimal(a)%}{{"{:.40e}".format(a) if a != 0.0 else "0.0"}}{% endmacro %}
#include <mm_malloc.h> //g++
......@@ -77,31 +79,31 @@ void {{codeNamespace}}::initQuadratureNodesAndWeights() {
nodes = (double *) _mm_malloc(sizeof(double)*{{nDofPad}}, ALIGNMENT);
{% for i in w1_seq %}
weights1[{{i}}] = {{"{:.17e}".format(weights1[i])}};
weights1[{{i}}] = {{formatDecimal(weights1[i])}};
{% endfor %}
{% for i in w2_seq %}
weights2[{{i}}] = {{"{:.17e}".format(weights2[i])}};
weights2[{{i}}] = {{formatDecimal(weights2[i])}};
{% endfor %}
{% for i in w3_seq %}
weights3[{{i}}] = {{"{:.17e}".format(weights3[i])}};
weights3[{{i}}] = {{formatDecimal(weights3[i])}};
{% endfor %}
{% for i in w4_seq %}
weights4[{{i}}] = {{"{:.17e}".format(weights4[i])}};
weights4[{{i}}] = {{formatDecimal(weights4[i])}};
{% endfor %}
{% for i in w1_seq %}
iweights1[{{i}}] = {{"{:.17e}".format(iweights1[i])}};
iweights1[{{i}}] = {{formatDecimal(iweights1[i])}};
{% endfor %}
{% for i in w3_seq %}
iweights3[{{i}}] = {{"{:.17e}".format(iweights3[i])}};
iweights3[{{i}}] = {{formatDecimal(iweights3[i])}};
{% endfor %}
{% for i in quadrature_seq %}
nodes[{{i}}] = {{"{:.17e}".format(QuadratureNodes[i])}};
nodes[{{i}}] = {{formatDecimal(QuadratureNodes[i])}};
{% endfor %}
......@@ -120,31 +122,31 @@ void {{codeNamespace}}::initQuadratureNodesAndWeights() {
{% endif %}
{% for i in w1_seq %}
weights1_SP[{{i}}] = {{"{:.17e}".format(weights1[i])}};
weights1_SP[{{i}}] = {{formatDecimal(weights1[i])}};
{% endfor %}
{% for i in w2_seq %}
weights2_SP[{{i}}] = {{"{:.17e}".format(weights2[i])}};
weights2_SP[{{i}}] = {{formatDecimal(weights2[i])}};
{% endfor %}
{% for i in w3_seq %}
weights3_SP[{{i}}] = {{"{:.17e}".format(weights3[i])}};
weights3_SP[{{i}}] = {{formatDecimal(weights3[i])}};
{% endfor %}
{% for i in w4_seq %}
weights4_SP[{{i}}] = {{"{:.17e}".format(weights4[i])}};
weights4_SP[{{i}}] = {{formatDecimal(weights4[i])}};
{% endfor %}
{% for i in w1_seq %}
iweights1_SP[{{i}}] = {{"{:.17e}".format(iweights1[i])}};
iweights1_SP[{{i}}] = {{formatDecimal(iweights1[i])}};
{% endfor %}
{% for i in w3_seq %}
iweights3_SP[{{i}}] = {{"{:.17e}".format(iweights3[i])}};
iweights3_SP[{{i}}] = {{formatDecimal(iweights3[i])}};
{% endfor %}
{% for i in quadrature_seq %}
nodes_SP[{{i}}] = {{"{:.17e}".format(QuadratureNodes[i])}};
nodes_SP[{{i}}] = {{formatDecimal(QuadratureNodes[i])}};
{% endfor %}
{% endif %}{# useSinglePrecision #}
......@@ -156,15 +158,15 @@ void {{codeNamespace}}::initQuadratureNodesAndWeights() {
{% if kernelType=="limiter" %}
{% for i in uh2lob_seq %}
uh2lob[{{i}}] = {{"{:.17e}".format(uh2lob[i])}};
uh2lob[{{i}}] = {{formatDecimal(uh2lob[i])}};
{% endfor %}
{% for i in dg2fv_seq %}
dg2fv[{{i}}] = {{"{:.17e}".format(dg2fv[i])}};
dg2fv[{{i}}] = {{formatDecimal(dg2fv[i])}};
{% endfor %}
{% for i in fv2dg_seq %}
fv2dg[{{i}}] = {{"{:.17e}".format(fv2dg[i])}};
fv2dg[{{i}}] = {{formatDecimal(fv2dg[i])}};
{% endfor %}
{% endif %}{# limiter #}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment