Merge pull request #78 from kjetil-lye/master

Adding entropy functional

alsfvm/include/alsfvm/functional/LogEntropy.hpp

@@ -0,0 +1,87 @@
+/* Copyright (c) 2018 ETH Zurich, Kjetil Olsen Lye
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+#include "alsfvm/functional/Functional.hpp"
+
+namespace alsfvm {
+namespace functional {
+
+//!
+//! Computes the thermodynamic entropy, see
+//!
+//! Fjordholm, U. S., Mishra, S., & Tadmor, E. (2012). Arbitrarily high-order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws, 50(2), 544–573.
+//!
+//!  for a definition (2.16), E.
+//!
+//!  Note that we compute the spatial integral of the entropy at all points.
+//!
+//! Specifically, we set
+//!
+//! \f[ s = \log(p)-\gamma\log(p)\f]
+//!
+//! where \f$\gamma>0\f$ is the gas constant. We set the entropy $E$ to be
+//!
+//! \f[E=\frac{-\rho s}{\gamma-1}\f]
+//!
+//! This functional computes
+//!
+//! \f[\int_D E \;dx\f]
+//!
+//!
+//!
+//!
+class LogEntropy : public Functional {
+public:
+
+    //! Uses no parameter
+    LogEntropy(const Parameters& parameters);
+
+    //! Computes the operator value on the given input data
+    //!
+    //! @note In order to support time integration, the result should be
+    //!       added to conservedVolumeOut and extraVolumeOut, not overriding
+    //!       it.
+    //!
+    //! @param[out] conservedVolumeOut at the end, should have the contribution
+    //!             of the functional for the conservedVariables
+    //!
+    //!
+    //! @param[in] conservedVolumeIn the state of the conserved variables
+    //!
+    //! @param[in] weight the current weight to be applied to the functional. Ie, the functional should compute
+    //!                   \code{.cpp}
+    //!                   conservedVolumeOut += weight + f(conservedVolumeIn)
+    //!                   \endcode
+    //!
+    //! @param[in] grid the grid to use
+    //!
+    virtual void operator()(volume::Volume& conservedVolumeOut,
+        const volume::Volume& conservedVolumeIn,
+        const real weight,
+        const grid::Grid& grid
+    ) override;
+
+    //! Returns grid.getDimensions()
+    virtual ivec3 getFunctionalSize(const grid::Grid& grid) const override;
+
+
+private:
+
+    const real gamma;
+
+};
+} // namespace functional
+} // namespace alsfvm

alsfvm/include/alsfvm/functional/LogEntropyCUDA.hpp

@@ -0,0 +1,94 @@
+/* Copyright (c) 2018 ETH Zurich, Kjetil Olsen Lye
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+#include "alsfvm/functional/Functional.hpp"
+#include "alsfvm/cuda/CudaMemory.hpp"
+
+
+namespace alsfvm {
+namespace functional {
+
+//!
+//! Computes the thermodynamic entropy, see
+//!
+//! Fjordholm, U. S., Mishra, S., & Tadmor, E. (2012). Arbitrarily high-order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws, 50(2), 544–573.
+//!
+//!  for a definition (2.16), E.
+//!
+//!  Note that we compute the spatial integral of the entropy at all points.
+//!
+//! Specifically, we set
+//!
+//! \f[ s = \log(p)-\gamma\log(p)\f]
+//!
+//! where \f$\gamma>0\f$ is the gas constant. We set the entropy $E$ to be
+//!
+//! \f[E=\frac{-\rho s}{\gamma-1}\f]
+//!
+//! This functional computes
+//!
+//! \f[\int_D E \;dx\f]
+//!
+//!
+//!
+//!
+class LogEntropyCUDA : public Functional {
+public:
+
+    //! Uses no parameter
+    LogEntropyCUDA(const Parameters& parameters);
+
+    //! Computes the operator value on the given input data
+    //!
+    //! @note In order to support time integration, the result should be
+    //!       added to conservedVolumeOut and extraVolumeOut, not overriding
+    //!       it.
+    //!
+    //! @param[out] conservedVolumeOut at the end, should have the contribution
+    //!             of the functional for the conservedVariables
+    //!
+    //!
+    //! @param[in] conservedVolumeIn the state of the conserved variables
+    //!
+    //! @param[in] weight the current weight to be applied to the functional. Ie, the functional should compute
+    //!                   \code{.cpp}
+    //!                   conservedVolumeOut += weight + f(conservedVolumeIn)
+    //!                   \endcode
+    //!
+    //! @param[in] grid the grid to use
+    //!
+    virtual void operator()(volume::Volume& conservedVolumeOut,
+        const volume::Volume& conservedVolumeIn,
+        const real weight,
+        const grid::Grid& grid
+    ) override;
+
+    //! Returns grid.getDimensions()
+    virtual ivec3 getFunctionalSize(const grid::Grid& grid) const override;
+
+
+private:
+
+    const real gamma;
+
+    alsfvm::shared_ptr<cuda::CudaMemory<real>> buffer = nullptr;
+    alsfvm::shared_ptr<cuda::CudaMemory<real>> bufferOut = nullptr;
+    alsfvm::shared_ptr<cuda::CudaMemory<real>> temporaryReductionMemory = nullptr;
+    size_t temporaryReductionMemoryStorageSizeBytes = 0;
+
+};
+} // namespace functional
+} // namespace alsfvm

alsfvm/src/config/SimulatorSetup.cpp

@@ -81,6 +81,22 @@ vec3<T> parseVector(const std::string& vectorAsString) {
     return vec3<T>(a, b, c);
 
 }
+
+template<class ParameterClass>
+void readAllEquationParameters(const boost::property_tree::ptree&
+    configuration, ParameterClass& parameters) {
+    auto fvmNode = configuration.get_child("fvm");
+
+    if (fvmNode.find("equationParameters") != fvmNode.not_found()) {
+
+        for (const auto& node : fvmNode.get_child("equationParameters")) {
+            const auto name = node.first;
+            const auto value = node.second.data();
+
+            parameters.addStringParameter(name, value);
+        }
+    }
+}
 }
 
 std::pair<alsfvm::shared_ptr<simulator::Simulator>,
@@ -489,7 +505,8 @@ void SimulatorSetup::readEquationParameters(const SimulatorSetup::ptree&
     }
 }
 
-alsfvm::shared_ptr<diffusion::DiffusionOperator> SimulatorSetup::createDiffusion(
+alsfvm::shared_ptr<diffusion::DiffusionOperator>
+SimulatorSetup::createDiffusion(
     const SimulatorSetup::ptree& configuration,
     const grid::Grid& grid,
     const simulator::SimulatorParameters& simulatorParameters,
@@ -562,6 +579,8 @@ std::vector<io::WriterPointer> SimulatorSetup::createFunctionals(
                     io::Parameters(functional.second.get_child("writer")));
             functional::Functional::Parameters parameters(functional.second);
 
+            readAllEquationParameters(configuration, parameters);
+
             auto functionalPointer = functionalFactory.makeFunctional(this->readPlatform(
                         configuration), name, parameters);
 

alsfvm/src/functional/Identity.cpp

@@ -28,14 +28,7 @@ void Identity::operator()(volume::Volume& conservedVolumeOut,
     const real weight,
     const grid::Grid& grid) {
 
-    const auto lengths = grid.getCellLengths();
-
-    if (lengths.z > 1 || lengths.y == 1) {
-        THROW("For now, Legendre polynomials only support 2d, givne dimensions " <<
-            lengths);
-    }
-
-
+
     const auto ghostCells = conservedVolumeIn.getNumberOfGhostCells();
 
 

alsfvm/src/functional/LogEntropy.cpp

@@ -0,0 +1,86 @@
+#include "alsfvm/functional/LogEntropy.hpp"
+#include "alsfvm/volume/volume_foreach.hpp"
+#include "alsfvm/functional/register_functional.hpp"
+#include <iostream>
+namespace alsfvm {
+namespace functional {
+
+LogEntropy::LogEntropy(const Functional::Parameters& parameters)
+    : gamma(parameters.getDouble("gamma")) {
+
+
+}
+
+void LogEntropy::operator()(volume::Volume& conservedVolumeOut,
+    const volume::Volume& conservedVolumeIn, const real weight,
+    const grid::Grid& grid) {
+
+    const auto lengths = grid.getCellLengths();
+
+    const real dxdydz = lengths.x * lengths.y * lengths.z;
+
+
+    real integral = 0.0;
+
+    const auto& densityView = conservedVolumeIn.getScalarMemoryArea(
+            "rho")->getView();
+
+    const auto& energyView = conservedVolumeIn.getScalarMemoryArea(
+            "E")->getView();
+
+    const size_t numberOfComponents = grid.getActiveDimension();
+
+    std::vector<const real*> momentumPointers;
+
+
+    momentumPointers.push_back(
+        conservedVolumeIn.getScalarMemoryArea("mx")->getPointer());
+
+    if (numberOfComponents > 1) {
+        momentumPointers.push_back(
+            conservedVolumeIn.getScalarMemoryArea("my")->getPointer());
+    }
+
+    if (numberOfComponents > 2) {
+        momentumPointers.push_back(
+            conservedVolumeIn.getScalarMemoryArea("mz")->getPointer());
+    }
+
+
+
+
+
+
+    volume::for_each_midpoint(conservedVolumeIn, grid, [&](real, real, real,
+    size_t i) {
+
+        const real density = densityView.at(i);
+        const real energy = energyView.at(i);
+
+        double momentumSquared = 0.0;
+
+        for (size_t component = 0; component < numberOfComponents; ++component) {
+            const auto momentum = momentumPointers[component][i];
+            momentumSquared += momentum * momentum;
+        }
+
+        const real pressure = (gamma - 1) * (energy - 1.0 / (2 * density) *
+                momentumSquared);
+
+
+        const real s = log(pressure) - gamma * log(density);
+        const real E = (-density * s) / (gamma - 1);
+
+        integral += E * dxdydz;
+    });
+
+    conservedVolumeOut.getScalarMemoryArea("E")->getPointer()[0] += weight
+        * integral;
+}
+
+ivec3 LogEntropy::getFunctionalSize(const grid::Grid&) const {
+    return {1, 1, 1};
+}
+REGISTER_FUNCTIONAL(cpu, log_entropy, LogEntropy)
+}
+}