Major commit taking out numerical flux

Physics and numerical fluxes are now computed outside of the DG discretization. Partial implementation of vector-valued problems by adding the std::array structures Diffusive terms require a function that will evaluate the solution numerical flux and the viscous numerical flux. Current implementation of SIPG is working to recover optimal convergence orders of the solution but not of the output functional Suboptimal output error is perceived for diffusion (maybe, I get 2p) and convection-diffusion (definitely) TODO: Find the adjoint consistent boundary condition for the diffusion problem and convection-diffusion problem TODO: Find out why the linear solver is soooo much slower for the diffusion problems Start 1: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION 1/12 Test #1: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION .............. Passed 0.53 sec Start 2: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION 2/12 Test #2: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION .............. Passed 5.84 sec Start 3: 3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION 3/12 Test #3: 3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION .............. Passed 47.98 sec Start 4: 1D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION 4/12 Test #4: 1D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION .............. Passed 0.64 sec Start 5: 2D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION 5/12 Test #5: 2D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION .............. Passed 3.23 sec Start 6: 3D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION 6/12 Test #6: 3D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION .............. Passed 237.84 sec Start 7: 1D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION 7/12 Test #7: 1D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ... Passed 0.75 sec Start 8: 2D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION 8/12 Test #8: 2D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ... Passed 3.45 sec Start 9: 3D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION 9/12 Test #9: 3D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ... Passed 266.84 sec Start 10: 1D_numerical_flux_conservation 10/12 Test #10: 1D_numerical_flux_conservation ........................... Passed 0.27 sec Start 11: 2D_numerical_flux_conservation 11/12 Test #11: 2D_numerical_flux_conservation ........................... Passed 0.25 sec Start 12: 3D_numerical_flux_conservation 12/12 Test #12: 3D_numerical_flux_conservation ........................... Passed 0.24 sec 100% tests passed, 0 tests failed out of 12 Total Test time (real) = 567.87 sec
dougshidong · May 18, 2019 · 593089b · 593089b
1 parent 67e7119
commit 593089b
Show file tree

Hide file tree

Showing 23 changed files with 1,000 additions and 289 deletions.
diff --git a/TODO b/TODO
@@ -2,15 +2,11 @@ Look into member variable naming for fe_values
     This one just bit me.
     In face_boundary I used the member fe_values instead of passed fe_values_face by mistake.
 
-Take physics out of the DiscontinuousGalerkin class.
+Find the adjoint consistent boundary condition for the diffusion problem and convection-diffusion problem
 
-Move boundary conditions to Physics class.
-
-Compute auxiliary variable $ \sigma = \nabla u $ independently.
+Find out why the linear solver is soooo much slower for the diffusion problems. Maybe just need pre-conditioner.
 
-Take numerical flux out of the DiscontinuousGalerkin class.
-
-Adjoint consistency, output integral 2p.
+Move boundary conditions to Physics class.
 
 Vector-valued problems
     - Figure out global solution ordering

diff --git a/input_file.prm b/input_file.prm
@@ -0,0 +1,43 @@
+# Listing of Parameters
+# ---------------------
+# Number of dimensions
+set dimension = 1
+
+# The kind of solver for the linear system. Choices are
+# <advection|convection_diffusion>.
+set pde_type  = advection
+
+
+subsection ODE solver
+  # Maximum nonlinear solver iterations
+  set nonlinear_max_iterations            = 3
+
+  # Nonlinear solver residual tolerance
+  set nonlinear_steady_residual_tolerance = 1e-13
+
+  # Print every print_iteration_modulo iterations of the nonlinear solver
+  set print_iteration_modulo              = 1
+
+  # Explicit or implicit solverChoices are <explicit|implicit>.
+  set solver_type                         = implicit
+end
+
+
+subsection manufactured solution convergence study
+  # Last degree used for convergence study
+  set degree_end   = 4
+
+  # Starting degree for convergence study
+  set degree_start = 1
+
+  # Initial grid of size (initial_grid_size)^dim
+  set initial_grid_size = 2
+
+  # Number of grids in grid study
+  set number_of_grids   = 5
+
+  # Multiplier on grid size. nth-grid will be of size
+  # (initial_grid^grid_progression)^dim
+  set grid_progression  = 1.5
+end
+
diff --git a/run_this_input.py b/run_this_input.py
@@ -0,0 +1,28 @@
+#!/usr/bin/env python3
+import os
+import sys
+import re
+
+fname = 'input_file.prm'
+
+f = open(fname)
+
+regexp = re.compile(r'set dimension.*?([1-3.-]+)')
+
+dim = -1
+with open(fname) as f:
+    for line in f:
+        match = regexp.match(line)
+        if match:
+            dim = int(match.group(1))
+            break
+if dim == -1:
+    sys.exit("No valid 'set dimension = [1-3]' line found")
+
+
+
+bin_path = './bin/PHiLiP_'+str(dim)+'D'
+
+command = bin_path + ' -p ' + fname
+print("Running: " + command)
+os.system(command)
diff --git a/src/assemble_implicit.cpp b/src/assemble_implicit.cpp
diff --git a/src/dg.cpp b/src/dg.cpp
@@ -70,6 +70,8 @@ namespace PHiLiP
         pde_physics = PhysicsFactory<dim, nstate, ADtype >::create_Physics(parameters->pde_type);
         conv_num_flux = NumericalFluxFactory<dim, nstate, ADtype>::create_convective_numerical_flux
             (parameters->conv_num_flux_type, pde_physics);
+        diss_num_flux = NumericalFluxFactory<dim, nstate, ADtype>::create_dissipative_numerical_flux
+            (parameters->diss_num_flux_type, pde_physics);
     }
     // Destructor
     template <int dim, int nstate, typename real>
@@ -99,6 +101,7 @@ namespace PHiLiP
     void DiscontinuousGalerkin<dim,nstate,real>::delete_fe_values ()
     {
         delete conv_num_flux;
+        delete diss_num_flux;
         delete pde_physics;
 
         if (fe_values_cell          != NULL) delete fe_values_cell;
@@ -538,6 +541,44 @@ namespace PHiLiP
       data_out.write_gnuplot(gnuplot_output);
     }
 
+    template <int dim, int nstate, typename real>
+    void DiscontinuousGalerkin<dim,nstate,real>::evaluate_inverse_mass_matrices ()
+    {
+        // Invert and store mass matrix
+        // Using Gauss-Jordan since it's deal.II's default invert function
+        // Could store Cholesky decomposition for more efficient pre-processing
+        const int n_quad_pts      = quadrature.size();
+        const int n_dofs_per_cell = fe.dofs_per_cell;
+
+        typename DoFHandler<dim>::active_cell_iterator
+           cell = dof_handler.begin_active(),
+           endc = dof_handler.end();
+
+        inv_mass_matrix.resize(triangulation->n_active_cells(),
+                               FullMatrix<real>(n_dofs_per_cell));
+        FullMatrix<real> mass_matrix(n_dofs_per_cell);
+        for (; cell!=endc; ++cell) {
+
+            int cell_index = cell->index();
+            fe_values_cell->reinit(cell);
+
+            for (int itest=0; itest<n_dofs_per_cell; ++itest) {
+                for (int itrial=itest; itrial<n_dofs_per_cell; ++itrial) {
+                    mass_matrix[itest][itrial] = 0.0;
+                    for (int iquad=0; iquad<n_quad_pts; ++iquad) {
+                        mass_matrix[itest][itrial] +=
+                            fe_values_cell->shape_value(itest,iquad)
+                            * fe_values_cell->shape_value(itrial,iquad)
+                            * fe_values_cell->JxW(iquad);
+                    }
+                    mass_matrix[itrial][itest] = mass_matrix[itest][itrial];
+                }
+            }
+            inv_mass_matrix[cell_index].invert(mass_matrix);
+        }
+    }
+
+
 
 
     template class DiscontinuousGalerkin <PHILIP_DIM, 1, double>;

diff --git a/src/dg.h b/src/dg.h
@@ -44,6 +44,7 @@ namespace PHiLiP
 
         void allocate_system ();
         void assemble_system ();
+        void evaluate_inverse_mass_matrices ();
         double get_residual_l2norm ();
 
         void delete_fe_values ();
@@ -84,6 +85,7 @@ namespace PHiLiP
         /// Contains the physics of the PDE
         Physics<dim, nstate, Sacado::Fad::DFad<real> >  *pde_physics;
         NumericalFluxConvective<dim, nstate, Sacado::Fad::DFad<real> > *conv_num_flux;
+        NumericalFluxDissipative<dim, nstate, Sacado::Fad::DFad<real> > *diss_num_flux;
 
 
         void allocate_system_explicit ();
@@ -141,6 +143,8 @@ namespace PHiLiP
             Vector<real>          &current_cell_rhs,
             Vector<real>          &neighbor_cell_rhs);
 
+        std::vector<FullMatrix<real>> inv_mass_matrix;
+
         // QGauss is Gauss-Legendre quadrature nodes
         const QGauss<dim>   quadrature;
         const QGauss<dim-1> face_quadrature;

diff --git a/src/grid_study.cpp b/src/grid_study.cpp
@@ -90,7 +90,7 @@ namespace PHiLiP
 
             // p0 tends to require a finer grid to reach asymptotic region
             unsigned int n_grids = n_grids_input;
-            if (poly_degree == 0) n_grids = n_grids_input + 2;
+            if (poly_degree <= 1) n_grids = n_grids_input + 1;
 
             std::vector<int> n_1d_cells(n_grids);
             n_1d_cells[0] = initial_grid_size;
@@ -168,6 +168,8 @@ namespace PHiLiP
                 //dg.set_physics(physics);
                 dg.allocate_system ();
 
+                dg.evaluate_inverse_mass_matrices();
+
                 //ODESolver<dim, double> *ode_solver = ODESolverFactory<dim, double>::create_ODESolver(parameters.solver_type);
                 ODESolver<dim, nstate, double> *ode_solver = ODESolverFactory<dim, nstate, double>::create_ODESolver(&dg);
 

diff --git a/src/linear_solver.cpp b/src/linear_solver.cpp
@@ -20,15 +20,18 @@ namespace PHiLiP
 
         //system_matrix.print(std::cout, true);
         //solution.print(std::cout);
-        //right_hand_side.print(std::cout);
-        //FullMatrix<double> fullA(system_matrix.m());
-        //fullA.copy_from(system_matrix);
-        //std::cout<<"Dense matrix:"<<std::endl;
-        //fullA.print_formatted(std::cout, 3, true, 10, "0", 1., 0.);
 
 
         using LinSolvParam = Parameters::LinearSolver;
         if (param->linear_solver_type == LinSolvParam::LinearSolverType::direct) {
+            if (param->linear_solver_output == Parameters::OutputType::verbose) {
+                right_hand_side.print(std::cout);
+                FullMatrix<double> fullA(system_matrix.m());
+                fullA.copy_from(system_matrix);
+                std::cout<<"Dense matrix:"<<std::endl;
+                fullA.print_formatted(std::cout, 3, true, 10, "0", 1., 0.);
+            }
+
             SolverControl solver_control(1, 0);
             TrilinosWrappers::SolverDirect::AdditionalData data(false);
             //TrilinosWrappers::SolverDirect::AdditionalData data(parameters.output == Parameters::Solver::verbose);

diff --git a/src/numerical_flux/CMakeLists.txt b/src/numerical_flux/CMakeLists.txt
@@ -1,4 +1,5 @@
 set(SOURCE
+    viscous_numerical_flux.cpp
     numerical_flux.cpp
     )
 

diff --git a/src/numerical_flux/numerical_flux.cpp b/src/numerical_flux/numerical_flux.cpp
@@ -1,6 +1,7 @@
 #include <Sacado.hpp>
 #include <deal.II/differentiation/ad/sacado_product_types.h>
 #include "numerical_flux.h"
+#include "viscous_numerical_flux.h"
 
 namespace PHiLiP
 {
@@ -87,42 +88,13 @@ namespace PHiLiP
         return numerical_flux_dot_n;
     }
 
-    template <int dim, int nstate, typename real>
-    NumericalFluxDissipative<dim,nstate,real>::~NumericalFluxDissipative() {}
-
-    template<int dim, int nstate, typename real>
-    void SymmetricInternalPenalty<dim,nstate,real>
-    ::evaluate_flux (
-        const std::array<real, nstate> &soln_int,
-        const std::array<real, nstate> &soln_ext,
-        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_int,
-        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_ext,
-        const Tensor<1,dim,real> &normal1,
-        std::array<real, nstate> &soln_flux,
-        std::array<Tensor<1,dim,real>, nstate> &grad_flux) const
-    {
-        std::array<real,nstate> soln_avg = array_average<nstate,real>(soln_int, soln_ext);
-
-        //std::array<Tensor<2,dim,real>,nstate> diffusion_matrix_int = 
-        //    pde_physics->diffusion_matrix(soln_int);
-
-        //std::array<Tensor<2,dim,real>,nstate> diffusion_matrix_ext = 
-        //    pde_physics->diffusion_matrix(soln_ext);
-    }
-
-
     // Instantiation
     template class NumericalFluxConvective<PHILIP_DIM, 1, double>;
     template class NumericalFluxConvective<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
 
     template class LaxFriedrichs<PHILIP_DIM, 1, double>;
     template class LaxFriedrichs<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
 
-    template class NumericalFluxDissipative<PHILIP_DIM, 1, double>;
-    template class NumericalFluxDissipative<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
-
-    template class SymmetricInternalPenalty<PHILIP_DIM, 1, double>;
-    template class SymmetricInternalPenalty<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
 
     template class NumericalFluxFactory<PHILIP_DIM, 1, double>;
     template class NumericalFluxFactory<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;

diff --git a/src/numerical_flux/numerical_flux.h b/src/numerical_flux/numerical_flux.h
@@ -3,6 +3,7 @@
 
 #include <deal.II/base/tensor.h>
 #include "physics/physics.h"
+#include "numerical_flux/viscous_numerical_flux.h"
 namespace PHiLiP
 {
     using namespace dealii;
@@ -46,56 +47,6 @@ namespace PHiLiP
 
     };
 
-    /// Numerical flux associated with dissipation
-    template<int dim, int nstate, typename real>
-    class NumericalFluxDissipative
-    {
-    public:
-    /// Base class destructor.
-    /// Abstract classes required virtual destructors.
-    /// Even if it is a pure function, its definition is still required.
-    virtual ~NumericalFluxDissipative() = 0;
-
-    virtual void evaluate_flux (
-        const std::array<real, nstate> &soln_int,
-        const std::array<real, nstate> &soln_ext,
-        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_int,
-        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_ext,
-        const Tensor<1,dim,real> &normal1,
-        std::array<real, nstate> &soln_flux,
-        std::array<Tensor<1,dim,real>, nstate> &grad_flux) const = 0;
-
-    };
-
-    template<int dim, int nstate, typename real>
-    class SymmetricInternalPenalty: public NumericalFluxDissipative<dim, nstate, real>
-    {
-    public:
-    /// Constructor
-    SymmetricInternalPenalty(Physics<dim, nstate, real> *physics_input)
-    :
-    pde_physics(physics_input)
-    {};
-    /// Destructor
-    ~SymmetricInternalPenalty() {};
-
-    /// Evaluate solution and gradient flux
-    /*  $\hat{u} = {u_h}$, 
-     *  $ \hat{A} = {{ A \nabla u_h }} - \mu {{ A }} [[ u_h ]] $
-     */
-    void evaluate_flux (
-        const std::array<real, nstate> &soln_int,
-        const std::array<real, nstate> &soln_ext,
-        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_int,
-        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_ext,
-        const Tensor<1,dim,real> &normal1,
-        std::array<real, nstate> &soln_flux,
-        std::array<Tensor<1,dim,real>, nstate> &grad_flux) const;
-    protected:
-    const Physics<dim, nstate, real> *pde_physics;
-
-    };
-
 
     template <int dim, int nstate, typename real>
     class NumericalFluxFactory