Merged Runge-Kutta

povanberg · Mar 20, 2019 · 33ba871 · 33ba871
1 parent 023d160
commit 33ba871
Show file tree

Hide file tree

Showing 6 changed files with 205 additions and 29 deletions.
diff --git a/doc/config/config.conf b/doc/config/config.conf
@@ -27,7 +27,7 @@ elementType=Lagrange
 
 # Time integration method:
 # ["Euler1", "Euler2", "Runge-Kutta"...]
-timeIntMethod=Euler1
+timeIntMethod=Runge-Kutta
 
 # Boundary condition:
 # physicalGroupName = Dirichelet 0 or Neumann 42 (<name> <value>)

diff --git a/include/solver.h b/include/solver.h
@@ -6,7 +6,9 @@
 
 namespace solver {
 
+    void dtfu(Mesh &mesh, Config config, std::vector<double> &u, std::vector<double> &a, double beta, int N);
     void solveForwardEuler(Mesh &mesh, Config config);
+    void solveRungeKutta(Mesh &mesh, Config config);
 }
 
 #endif
diff --git a/include/utils.h b/include/utils.h
@@ -26,11 +26,13 @@ namespace lapack {
     double dot(double *A, double *b, int N);
 
     // Matrix/vector product:  y := alpha*A*x + y
-    void linEq(double *A, double *X, double *Y, double &alpha, int &N);
+    void linEq(double *A, double *X, double *Y, double &alpha, double beta, int &N);
 
     void minus(double *A, double *B, int N);
 
     void plus(double *A, double *B, int N);
+
+    void plusTimes(double *A, double *B, double c, int N);
 }
 
 namespace eigen {

diff --git a/src/dgalerkin.cpp b/src/dgalerkin.cpp
@@ -31,7 +31,12 @@ int main(int argc, char **argv)
     // Discontinuous Galerkin simulation.
     Mesh mesh(msh_name, config);
 
-    solver::solveForwardEuler(mesh, config);
+    if(config.timeIntMethod == "Euler1")
+        solver::solveForwardEuler(mesh, config);
+    else if(config.timeIntMethod == "Runge-Kutta")
+        solver::solveRungeKutta(mesh, config);
+    else
+        solver::solveForwardEuler(mesh, config);
 
     //gmsh::fltk::run();
     gmsh::finalize();

diff --git a/src/solver.cpp b/src/solver.cpp
@@ -14,6 +14,22 @@ namespace solver {
         //return 1;
     }
 
+    void dtfu(Mesh &mesh, Config config, std::vector<double> &u_next, std::vector<double> &a, double beta, int N){
+
+        // Element variables
+        double elFlux[mesh.getElNumNodes()];
+        double elStiffvector[mesh.getElNumNodes()];
+
+        mesh.precomputeFlux(a.data(), u_next.data());
+        for(int el=0; el<mesh.getElNum(); ++el){
+
+            mesh.getElFlux(el, elFlux);
+            mesh.getElStiffVector(el, a.data(), u_next.data(), elStiffvector);
+            lapack::minus(elStiffvector, elFlux, N);
+            lapack::linEq(&mesh.elMassMatrix(el), &elStiffvector[0], &u_next[el*N], config.timeStep, beta, N);
+        }
+    }
+
     void solveForwardEuler(Mesh &mesh, Config config) {
 
         // Solution vectors
@@ -43,10 +59,6 @@ namespace solver {
         mesh.precomputeMassMatrix();
         mesh.setNumFlux(config.flux, a.data(), config.fluxCoeff);
 
-        // Element variables
-        double elFlux[mesh.getElNumNodes()];
-        double elStiffvector[mesh.getElNumNodes()];
-
         int step = 0;
         for(double t=config.timeStart, tDisplay =0; t<config.timeEnd; t+=config.timeStep, tDisplay+=config.timeStep, ++step) {
 
@@ -65,21 +77,148 @@ namespace solver {
                                     "s] Step number : " + std::to_string(step));
             }
 
-            mesh.precomputeFlux(a.data(), u.data());
+            dtfu(mesh, config, u_next, a, 1.0, N);
+        }
+
+        gmsh::view::write(viewTag, "data.msh");
+    }
+
+    void solveRungeKutta(Mesh &mesh, Config config) {
+
+        // Solution vectors
+        int N = mesh.getElNumNodes();
+        std::vector<int> nodeTags(&mesh.elNodeTag(0), &mesh.elNodeTag(0) + mesh.getNumNodes());
+        std::vector<double> u(mesh.getNumNodes());
+        std::vector<double> u_next(mesh.getNumNodes());
+
+        // Save results initialization
+        int viewTag = gmsh::view::add("Results");
+        std::vector<std::string> names;
+        gmsh::model::list(names);
+        std::vector<std::vector<double>> solution(nodeTags.size(), std::vector<double>(1));
+
+        // Set Initial conditions
+        std::vector<double> coord;
+        std::vector<double> parametricCoord;
+        for(int n=0; n<mesh.getNumNodes(); n++) {
+            gmsh::model::mesh::getNode(nodeTags[n], coord, parametricCoord);
+            u_next[n] = f(coord);
+        }
 
-            for(int el=0; el<mesh.getElNum(); ++el) {
+        // Convection vector
+        std::vector<double> a = {3, 0, 0};
 
-                mesh.getElFlux(el, elFlux);
-                mesh.getElStiffVector(el, a.data(), u.data(), elStiffvector);
+        // The mass matrix is invariant along iteration
+        mesh.precomputeMassMatrix();
+        mesh.setNumFlux(config.flux, a.data(), config.fluxCoeff);
+
+        int step = 0;
+        for(double t=config.timeStart, tDisplay =0; t<config.timeEnd; t+=config.timeStep, tDisplay+=config.timeStep, ++step) {
+
+            // Note that Neumann BCs are directly incorporated in flux calculation
+            mesh.enforceDiricheletBCs(u_next.data());
 
-                //display::print(std::vector<double>(&mesh.elMassMatrix(el),&mesh.elMassMatrix(el)+N*N), N, true);
-                //std::cout << "----" << std::endl;
+            u = u_next;
+
+            // Savings
+            if(step==0 || tDisplay>=config.timeRate){
+                tDisplay = 0;
+                for (int i = 0; i < nodeTags.size(); ++i)
+                    solution[i][0] = u[i];
+                gmsh::view::addModelData(viewTag, step, names[0], "NodeData", nodeTags, solution, t, 1);
+                gmsh::logger::write("[" + std::to_string(t) + "/" + std::to_string(config.timeEnd) +
+                                    "s] Step number : " + std::to_string(step));
+            }
 
-                lapack::minus(elStiffvector, elFlux, N);
-                lapack::linEq(&mesh.elMassMatrix(el), &elStiffvector[0], &u_next[el*N], config.timeStep, N);
+            // k vectors
+            std::vector<double> k1 = u_next;
+            std::vector<double> k2 = u_next;
+            std::vector<double> k3 = u_next;
+            std::vector<double> k4 = u_next;
+
+            dtfu(mesh, config, k1, a, 0, N); // k1
+            lapack::plusTimes(k2.data(), k1.data(), 0.5, k2.size()); // k2
+            dtfu(mesh, config, k2, a, 0, N);
+            lapack::plusTimes(k3.data(), k2.data(), 0.5, k3.size()); // k3
+            dtfu(mesh, config, k3, a, 0, N);
+            lapack::plusTimes(k4.data(), k3.data(), 1.0, k4.size()); // k4
+            dtfu(mesh, config, k4, a, 0, N);
+
+            for(int i=0; i<u_next.size(); ++i){
+            u_next[i] += (k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
             }
         }
 
         gmsh::view::write(viewTag, "data.msh");
     }
+/*
+    void solveRungeKutta(Mesh &mesh, Config config) {
+
+        // Solution vectors
+        int N = mesh.getElNumNodes();
+        std::vector<int> nodeTags(&mesh.elNodeTag(0), &mesh.elNodeTag(0) + mesh.getNumNodes());
+        std::vector<double> u(mesh.getNumNodes());
+        std::vector<double> u_next(mesh.getNumNodes());
+
+        // Save results initialization
+        int viewTag = gmsh::view::add("Results");
+        std::vector<std::string> names;
+        gmsh::model::list(names);
+        std::vector<std::vector<double>> solution(nodeTags.size(), std::vector<double>(1));
+
+        // Set Initial conditions
+        std::vector<double> coord;
+        std::vector<double> parametricCoord;
+        for(int n=0; n<mesh.getNumNodes(); n++) {
+            gmsh::model::mesh::getNode(nodeTags[n], coord, parametricCoord);
+            u_next[n] = f(coord);
+        }
+
+        // Convection vector
+        std::vector<double> a = {1, 1, 0};
+
+        // The mass matrix is invariant along iteration
+        mesh.precomputeMassMatrix();
+
+        int step = 0;
+        for(double t=config.timeStart, tDisplay =0; t<config.timeEnd; t+=config.timeStep, tDisplay+=config.timeStep, ++step) {
+
+            // Note that Neumann BCs are directly incorporated in flux calculation
+            mesh.enforceDiricheletBCs(u_next.data());
+
+            u = u_next;
+
+            // Savings
+            if(step==0 || tDisplay>=config.timeRate){
+                tDisplay = 0;
+                for (int i = 0; i < nodeTags.size(); ++i)
+                    solution[i][0] = u[i];
+                gmsh::view::addModelData(viewTag, step, names[0], "NodeData", nodeTags, solution, t, 1);
+                gmsh::logger::write("[" + std::to_string(t) + "/" + std::to_string(config.timeEnd) +
+                                    "s] Step number : " + std::to_string(step));
+            }
+
+            // k vectors
+            std::vector<double> k1 = u_next;
+            std::vector<double> k2 = u_next;
+            std::vector<double> k3 = u_next;
+            std::vector<double> k4 = u_next;
+
+            dtfu(mesh, config, k1, a, 0, N); // k1
+            lapack::plusTimes(k2.data(), k1.data(), 0.5, k2.size()); // k2
+            dtfu(mesh, config, k2, a, 0, N);
+            lapack::plusTimes(k3.data(), k2.data(), 0.5, k3.size()); // k3
+            dtfu(mesh, config, k3, a, 0, N);
+            lapack::plusTimes(k4.data(), k3.data(), 1.0, k4.size()); // k4
+            dtfu(mesh, config, k4, a, 0, N);
+
+            for(int i=0; i<u_next.size(); ++i){
+                //u_next[i] += (k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
+                u_next[i] += k1[i];
+            }
+        }
+
+        gmsh::view::write(viewTag, "data.msh");
+
+    }*/
 }
diff --git a/src/utils.cpp b/src/utils.cpp
@@ -1,6 +1,7 @@
 #include <iostream>
 #include <iomanip>
 #include <Eigen/Dense>
+#include <iomanip>
 
 extern "C" {
     // LU decomoposition of a general matrix
@@ -63,11 +64,10 @@ namespace lapack {
         for(int i=0; i<N; i++){A[i] /= norm;};
     }
 
-    // Matrix/vector product:  y := alpha*A*x + y
-    void linEq(double *A, double *X, double *Y, double &alpha, int &N){
+    // Matrix/vector product:  y := alpha*A*x + beta*y
+    void linEq(double *A, double *X, double *Y, double &alpha, double beta, int &N){
         char TRANS = 'T';
         int INC = 1;
-        double beta=1;
         dgemv_(TRANS, N, N, alpha, A, N, X, INC, beta, Y, INC);
     }
 
@@ -78,11 +78,20 @@ namespace lapack {
     }
 
     void minus(double *A, double *B, int N) {
-        std::transform(A, A + N, B, A, std::minus<double>());
+        for(int i=0; i<N; ++i)
+            A[i] -= B[i];
+        //std::transform(A, A + N, B, A, std::minus<double>());
     }
 
     void plus(double *A, double *B, int N) {
-        std::transform(A, A + N, B, A, std::plus<double>());
+        for(int i=0; i<N; ++i)
+            A[i] += B[i];
+        //std::transform(A, A + N, B, A, std::plus<double>());
+    }
+
+    void plusTimes(double *A, double *B, double c, int N) {
+        for(int i=0; i<N; ++i)
+            A[i] += B[i]*c;
     }
 }
 
@@ -98,14 +107,33 @@ namespace eigen {
         B_eigen = A_eigen.lu().solve(B_eigen);
     }
 
-    void inverse(double *A, int &N) {
-        Eigen::Map<Eigen::MatrixXd> A_eigen(A, N, N);
-        Eigen::JacobiSVD<Eigen::MatrixXd> svd(A_eigen);
-        double cond = svd.singularValues()(0)
-                      / svd.singularValues()(svd.singularValues().size()-1);
-        std::cout << "cond " << cond << std::endl;
-        A_eigen = A_eigen.inverse();
-    }
-
 }
 
+namespace display {
+    template<typename Container>
+    void print(const Container& cont, int row = 1, bool colMajor=false) {
+
+        if(colMajor){
+            for(int rowIt=0; rowIt<row; ++rowIt){
+                int colIt = 0;
+                for (auto const& x : cont) {
+                    if(colIt%row == rowIt) {
+                        std::cout << std::setprecision(4) << std::left << std::setw(10) << x << " ";
+                    }
+                    colIt++;
+                }
+                std::cout << std::endl;
+            }
+        }
+        else {
+            int colIt = 0;
+            for (auto const& x : cont) {
+                std::cout << std::setprecision(4) << std::left << std::setw(10) << x << " ";
+                colIt++;
+                if(colIt%row == 0)
+                    std::cout << std::endl;
+            }
+        }
+
+    }
+}