navierstokes demo

demo/CMakeLists.txt

@@ -105,6 +105,16 @@ if (TOOLS_DIR)
     add_executable("drivenCavity" ${drivenCavity})
     target_link_libraries("drivenCavity" ${BASIS_LIBS})
 
+    include_directories(/home/spraetor/projects/src/common)
+    set(navierStokes ${TOOLS_DIR}/diffuseDomain/POperators.cc
+# 		    ${TOOLS_DIR}/baseProblems/NavierStokes_TaylorHood.cc
+		    ${TOOLS_DIR}/baseProblems/NavierStokes_TaylorHood_RB.cc
+		    ${TOOLS_DIR}/baseProblems/time/ExtendedRosenbrockStationary.cc
+		    /home/spraetor/projects/src/common/GeometryTools.cc
+		    src/navierStokes.cc)
+    add_executable("navierStokes" ${navierStokes})
+    target_link_libraries("navierStokes" ${BASIS_LIBS})
+
   endif()
 else()
     message(WARNING "No tools directory specified! Some demos will not be build.")

demo/init/elliptImplicit.dat.2d

@@ -12,8 +12,8 @@ mesh->refinement->interface width: 0.01
 ellipt->mesh:                   elliptMesh
 ellipt->dim:                    2
 ellipt->components:             1
-ellipt->polynomial degree[0]:   1
- 
+ellipt->polynomial degree[0]:   2
+
 ellipt->solver:                 cg
 ellipt->solver->ell:        1
 ellipt->solver->max iteration:  1000
@@ -33,7 +33,7 @@ ellipt->marker[0]->MSGamma:        0.5
 ellipt->adapt[0]->tolerance:          1e-6
 ellipt->adapt[0]->refine bisections:  2
 
-ellipt->adapt->max iteration:         5
+ellipt->adapt->max iteration:         10
 
 ellipt->output->filename:       output/elliptImplicit.2d
 ellipt->output->ParaView format: 1

demo/init/navierStokes.dat.2d

+dimension of world: 2
+
+% ====================== VARIABLES ========================
+output_folder: output
+output_postfix: _ns
+mesh_name: mesh
+polynomial-degree: 1
+
+% ====================== MESH =============================
+mesh->H: 4.1
+
+% ====================== MESH =============================
+${mesh_name}->macro file name: macro/kanal_msh.2d
+${mesh_name}->global refinements: 0
+${mesh_name}->check: 0
+
+% ====================== INCLUDES =========================
+#include "init/navierStokes_TaylorHood.inc.2d"
+
+% ====================== USER_PARAMETER - NS ==============
+
+ns->viscosity: 1/100
+ns->theta: 0.5
+%ns->force: [0.0, 9.81]		% gravitational force [m/s^2]
+
+ns->force dirichlet bc: 0
+
+ns->Um: 1.5
+
+ns->initial velocity: 0
+ns->initial velocity value: 0.0
+ns->laplace operator: 0  	% 0... div(nu*grad(u)), 1... div(0.5*nu*(grad(u)+grad(u)^T)) [sehr langsam]
+ns->non-linear term: 2 		% 0... u^old*grad(u_i^old), 1... u'*grad(u_i^old), 2... u^old*grad(u'_i)
+
+% ====================== TIMESTEPS ========================
+adapt->max iteration: 1
+adapt->max timestep iteration: 1
+adapt->max time iteration: 1
+
+adapt->timestep: 5.e-4
+adapt->max timestep: 1e+10
+adapt->min timestep: 1e-6
+adapt->start time: 0.0
+adapt->end time: 10
+
+% ====================== ESTIMATORS =======================
+adapt->strategy: 0   % 0=explicit, 1=implicit
+
+WAIT: 1
+

demo/init/navierStokes_Chorin.inc.2d

+${ns}->predictor->mesh: ${mesh_name}
+${ns}->pressure->mesh:  ${mesh_name}
+${ns}->corrector->mesh: ${mesh_name}
+
+% ============ USER_PARAMETERS - NS =================================
+${ns}->viscosity: 1
+${ns}->beta: 1
+${ns}->sigma: 0.072
+${ns}->exponent: 2
+%${ns}->force: [0.0, -9.81] % gravitational force
+
+${ns}->viscosity1: 1
+${ns}->viscosity2: ${${ns}->viscosity1}/1000
+
+${ns}->force dirichlet bc: 0
+${ns}->poisson problem pertubation: 0 % 1..applySingularPertubation, 2..applySingularDBC, 0..nothing
+${ns}->simple algorithm: 1
+${ns}->calculate pressure: 0
+${ns}->initial velocity: 3
+${ns}->initial velocity value: 0.1
+${ns}->laplace operator: 1  % 0... div(nu*grad(u)), 1... div(0.5*nu*(grad(u)+grad(u)^T)) [sehr langsam]
+${ns}->non-linear term: 2 % 0... u^old*grad(u_i^old), 1... u'*grad(u_i^old), 2... u^old*grad(u'_i)
+
+% =========== OUTPUT ==============================================
+
+${ns}->predictor->output->filename: ${output_folder}/velocity/predictor${output_postfix}_
+${ns}->pressure->output->filename: ${output_folder}/velocity/pressure${output_postfix}_
+${ns}->corrector->output->filename: ${output_folder}/velocity/velocity${output_postfix}_
+
+% ============= PROBLEM-SPACES ==================================
+${ns}->predictor->components: 2
+${ns}->pressure->components: 1
+${ns}->corrector->components: 2
+
+${ns}->predictor->polynomial degree[0]: ${polynomial-degree}
+${ns}->predictor->polynomial degree[1]: ${polynomial-degree}
+${ns}->pressure->polynomial degree[0]: ${polynomial-degree}
+${ns}->corrector->polynomial degree[0]: ${polynomial-degree}
+${ns}->corrector->polynomial degree[1]: ${polynomial-degree}
+
+${ns}->predictor->dim: 2
+${ns}->pressure->dim: 2
+${ns}->corrector->dim: 2
+
+${ns}->pressure->name: pressure
+${ns}->corrector->name: velocity
+
+% ================== SOLVER ======================================
+${ns}->predictor->solver: umfpack
+${ns}->predictor->solver->symmetric strategy: 0
+${ns}->predictor->solver->ell:               3
+${ns}->predictor->solver->max iteration:     500
+${ns}->predictor->solver->restart:           10  %  only used for GMRES
+${ns}->predictor->solver->tolerance:         1.e-8
+${ns}->predictor->solver->info:              1
+${ns}->predictor->solver->left precon:       ilu
+
+${ns}->predictor->petsc->ksp_type: preonly
+${ns}->predictor->petsc->pc_type: lu
+${ns}->predictor->petsc->pc_factor_mat_solver_package: mumps
+
+${ns}->pressure->solver: umfpack
+${ns}->pressure->solver->symmetric strategy: 0
+${ns}->pressure->solver->ell:               3
+${ns}->pressure->solver->max iteration:     500
+${ns}->pressure->solver->tolerance:         1.e-8
+${ns}->pressure->solver->info:              1
+${ns}->pressure->solver->left precon:       ilu
+
+${ns}->pressure->petsc->ksp_type: preonly
+${ns}->pressure->petsc->pc_type: lu
+${ns}->pressure->petsc->pc_factor_mat_solver_package: mumps
+
+${ns}->corrector->solver: umfpack
+${ns}->corrector->solver->symmetric strategy: 0
+${ns}->corrector->solver->ell:               3
+${ns}->corrector->solver->max iteration:     500
+${ns}->corrector->solver->tolerance:         1.e-8
+${ns}->corrector->solver->info:              1
+${ns}->corrector->solver->left precon:       ilu
+
+${ns}->corrector->petsc->ksp_type: preonly
+${ns}->corrector->petsc->pc_type: lu
+${ns}->corrector->petsc->pc_factor_mat_solver_package: mumps
+
+% =================== OUTPUT =========================================
+${ns}->predictor->output->ParaView animation: 0
+${ns}->predictor->output->ParaView format: 0
+${ns}->predictor->output->write every i-th timestep: 1
+%${ns}->predictor->output->compression:  gzip
+${ns}->predictor->output->append index: 1
+${ns}->predictor->output->index length: 7
+${ns}->predictor->output->index decimals: 5
+
+${ns}->pressure->output->ParaView animation: 0
+${ns}->pressure->output->ParaView format: 0
+${ns}->pressure->output->write every i-th timestep: 1
+%${ns}->pressure->output->compression:  gzip
+${ns}->pressure->output->append index: 1
+${ns}->pressure->output->index length: 7
+${ns}->pressure->output->index decimals: 5
+
+${ns}->corrector->output->ParaView animation: 1
+${ns}->corrector->output->ParaView vector format: 1
+${ns}->corrector->output->write vector as 3d vector: 1
+${ns}->corrector->output->write every i-th timestep: 1
+%${ns}->corrector->output->compression:  gzip
+${ns}->corrector->output->append index: 1
+${ns}->corrector->output->index length: 7
+${ns}->corrector->output->index decimals: 5
+

demo/init/navierStokes_rb.dat.2d

+dimension of world: 2
+
+% ====================== VARIABLES ========================
+output_folder: output
+output_postfix: _refinement3
+mesh_name: mesh
+polynomial-degree: 1
+
+% ====================== MESH =============================
+mesh->H: 4.1
+
+obstacle->num vertices: 7
+obstacle->vertex[0]: [1.5, 2]
+obstacle->vertex[1]: [2, 2.5]
+% obstacle->vertex[2]: [4.5, 2]
+% obstacle->vertex[3]: [2, 1.5]
+obstacle->vertex[2]: [2.5, 2.05]
+obstacle->vertex[3]: [6.5, 2.05]
+obstacle->vertex[4]: [6.5, 1.95]
+obstacle->vertex[5]: [2.5, 1.95]
+obstacle->vertex[6]: [2, 1.5]
+
+mesh->refinement->initial level: 1
+mesh->refinement->level on interface: 6
+mesh->refinement->level in inner domain: 1
+mesh->refinement->level in outer domain: 1
+mesh->refinement->interface width: 0.5
+mesh->refinement->fade out width: 1.0
+
+% ====================== MESH =============================
+${mesh_name}->macro file name: macro/kanal_square_fin.2d
+${mesh_name}->global refinements: 0
+${mesh_name}->check: 0
+
+% ====================== INCLUDES =========================
+#include "init/navierStokes_TaylorHood.inc.2d"
+
+% ====================== USER_PARAMETER - NS ==============
+
+ns->viscosity: 1/500
+ns->theta: 0.5
+%ns->force: [0.0, 9.81]		% gravitational force [m/s^2]
+
+ns->force dirichlet bc: 0
+
+ns->Um: 1.5
+
+ns->initial velocity: 0
+ns->initial velocity value: 0.0
+ns->laplace operator: 0  	% 0... div(nu*grad(u)), 1... div(0.5*nu*(grad(u)+grad(u)^T)) [sehr langsam]
+ns->non-linear term: 2 		% 0... u^old*grad(u_i^old), 1... u'*grad(u_i^old), 2... u^old*grad(u'_i)
+
+% ====================== TIMESTEPS ========================
+
+adapt->rosenbrock method: rodasp
+adapt->fix first timesteps: 0
+adapt->rosenbrock->timestep study: 0
+adapt->rosenbrock->timestep study steps: 0
+adapt->rosenbrock->error weights: [1,1,1]
+
+adapt[0]->time tolerance: 1.e-3
+
+adapt->timestep:      5.e-4
+adapt->max timestep:  1e+1
+adapt->min timestep:  1e-4
+adapt->start time:    0.0
+adapt->end time:      25.0
+
+% ====================== ESTIMATORS =======================
+adapt->strategy: 0   % 0=explicit, 1=implicit
+
+%ns->space->estimator[0]:              residual
+%ns->space->estimator[0]->error norm:  1   % 1: H1_NORM, 2: L2_NORM
+%ns->space->estimator[0]->C0:          0.1 % constant of element residual
+%ns->space->estimator[0]->C1:          0.1 % constant of jump residual
+
+ns->space->marker[0]->strategy:       0   % 0: no adaption 1: GR 2: MS 3: ES 4:GERS
+
+adapt[0]->tolerance:          1e-3
+adapt->max iteration:      0
+
+WAIT: 1
+

demo/init/pfc.dat.2d

@@ -39,7 +39,11 @@ adapt->timestep:      1.e-1
 adapt->max timestep:  1e+10
 adapt->min timestep:  1e-10
 adapt->start time:    0.0
-adapt->end time:      50.0
+adapt->end time:      10000.0
+
+adapt->sequence->calc end time: 0
+adapt->sequence->timesteps:  [1.e-2,2.e-2,5.e-2,1.e-1]
+adapt->sequence->number of timesteps:  [500,1500,1500,1500]
 
 
 % =================== OUTPUT =========================================

demo/src/navierStokes.cc

+#include "AMDiS.h"
+#include "NavierStokes_TaylorHood_RB.h"
+#include "navierStokes.h"
+#include "time/ExtendedRosenbrockAdaptInstationary.h"
+#include "Refinement.h"
+#include "MeshFunction_Level.h"
+
+#include "boost/date_time/posix_time/posix_time.hpp"
+
+using namespace AMDiS;
+using namespace boost::posix_time;
+
+
+
+class NS_Channel : public NavierStokes_TaylorHood_RB
+{
+public:
+  NS_Channel(std::string name_) : NavierStokes_TaylorHood_RB(name_) {}
+
+  void solveInitialProblem(AdaptInfo *adaptInfo)
+  {
+    super::solveInitialProblem(adaptInfo);
+
+    size_t nVertices = 0;
+    WorldVector<double> x;
+    Parameters::get("obstacle->num vertices",nVertices);
+    std::vector<WorldVector<double> > v(nVertices,x);
+    for (size_t i = 0; i < nVertices; i++)
+      Parameters::get("obstacle->vertex["+boost::lexical_cast<std::string>(i)+"]",v[i]);
+    v.push_back(v[0]);
+
+    SignedDistRefinement refFunction(getMesh());
+    RefinementLevelCoords2 refinement(
+      getFeSpace(),
+      &refFunction,
+      new Polygon(v));
+
+    // initial refinement
+    refinement.refine(10);
+  }
+protected:
+  void fillBoundaryConditions()
+  { FUNCNAME("NS_DrivenCavity::fillBoundaryConditions()");
+    
+    AbstractFunction<double, WorldVector<double> > *zero = new AMDiS::Const<double, WorldVector<double> >(0.0);
+    size_t dow = Global::getGeo(WORLD);
+
+    // +------ 5 ------+
+    // |               |
+    // 2   # <--1      3
+    // |               |
+    // +------ 4 ------+
+
+
+//     x[0] = 2.0;
+//     x[1] = 2.0;
+#if DEBUG != 0
+  DOFVector<double>* polygon = new DOFVector<double>(getFeSpace(0), "polygon");
+  polygon->interpol(new Polygon(v));
+  VtkWriter::writeFile(polygon, "polygon.vtu");
+#endif
+    
+    /// at rigid wall: no-slip boundary condition
+    for (size_t i = 0; i < dow; i++) {
+//       getProblem(0)->addImplicitDirichletBC(*(new Polygon(v)), i, i, *zero);
+//       getProblem(0)->addSingularDirichletBC(x, i, i, *zero);
+      getProblem(0)->addDirichletBC(1, i, i, zero);
+      getProblem(0)->addDirichletBC(4, i, i, zero);
+      getProblem(0)->addDirichletBC(5, i, i, zero);
+    }
+
+    double H = 4.1;
+    double Um = 1.5;
+    Parameters::get("mesh->H",H);
+    Parameters::get("ns->Um",Um);
+
+    /// at left wall: prescribed velocity
+    getProblem(0)->addDirichletBC(2, 0, 0, new InflowBC(H,Um));
+    getProblem(0)->addDirichletBC(2, 1, 1, zero);
+  }
+};
+
+
+int main(int argc, char** argv)
+{ FUNCNAME("main");
+
+  AMDiS::init(argc, argv);
+
+  NS_Channel nsProb("ns");
+  nsProb.initialize(INIT_ALL | INIT_EXACT_SOLUTION);
+
+  // Adapt-Infos
+  AdaptInfo adaptInfo("adapt", nsProb.getNumComponents());
+
+  // adaption loop - solve ch-prob and ns-prob
+//   AdaptInstationary adaptInstat("adapt", nsProb, adaptInfo, nsProb, adaptInfo);
+  ExtendedRosenbrockAdaptInstationary<NS_Channel> adaptInstat("adapt", nsProb, adaptInfo, nsProb, adaptInfo);
+
+  ptime start_time = microsec_clock::local_time();
+  nsProb.initTimeInterface();
+  int error_code = adaptInstat.adapt(); 
+  time_duration td = microsec_clock::local_time()-start_time;
+
+  MSG("elapsed time= %d sec\n", td.total_seconds());
+
+  AMDiS::finalize();
+
+  return error_code;
+};

demo/src/navierStokes.h

+/** \file navierStokes.h */
+
+#ifndef NAVIER_STOKES_H
+#define NAVIER_STOKES_H
+
+#include "AMDiS.h"
+#include "GeometryTools.h"
+
+struct InflowBC : AbstractFunction<double, WorldVector<double> >
+{
+  InflowBC(double H_=4.1, double Um_=1.5) : H(H_), Um(Um_) {}
+  double operator()(const WorldVector<double> &x) const {
+    return 4.0 * Um * x[1] * (H - x[1]) / sqr(H);
+  }
+protected:
+  double H;
+  double Um;
+};
+
+
+class Polygon : public AbstractFunction<double, WorldVector<double> >
+{
+public:
+
+  Polygon(WorldVector<double> x0_, WorldVector<double> x1_, WorldVector<double> x2_, WorldVector<double> x3_)
+  {
+    vertices.push_back(x0_);
+    vertices.push_back(x1_);
+    vertices.push_back(x2_);
+    vertices.push_back(x3_);
+    vertices.push_back(x0_);
+  }
+
+  Polygon(std::vector<WorldVector<double> > xi_) : vertices(xi_) { }
+
+  double operator()(const WorldVector<double>& x) const
+  {
+    double result = 1.e15;
+    for (size_t i = 0; i < vertices.size()-1; i++)
+      result = std::min(result, meshconv::distance_point_line_2d(x.begin(), vertices[i].begin(), vertices[i+1].begin()));
+    return result * (meshconv::point_in_polygon(x.begin(), vertices) ? -1.0 : 1.0);
+  };
+
+private:
+  std::vector<WorldVector<double> > vertices;
+};
+
+#endif

demo/src/pfc.cc

 #include "AMDiS.h"
 #include "PhaseFieldCrystal_Base.h"
 #include "Helpers.h"
+#include "getMaxima.h"
 
 #include "boost/date_time/posix_time/posix_time.hpp"
 
@@ -87,6 +88,9 @@ int main(int argc, char** argv)
 
   MSG("elapsed time= %d sec\n", td.total_seconds());
 
+  // calc maxima
+  getMaxima(&adaptInfo, pfcProb.getProblem(), "output");
+
   AMDiS::finalize();
 
   return error_code;