copss.h

Go to the documentation of this file.

// Copyright (C) 2015-2016 Xujun Zhao, Jiyuan Li, Xikai Jiang

// This code 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 2.1 of the License, or (at your option) any later version.


// This code 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
// Lesser General Public License for more details.


// You should have received a copy of the GNU General Public
// License along with this code; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
#pragma once

// C++ includes
#include <iostream>
#include <algorithm>
#include <cstring>
#include <math.h>
#include <time.h>
#include <fstream> 
#include <tuple>
#include <stdlib.h>

// Libmesh includes 
#include "libmesh/libmesh.h"
#include "libmesh/libmesh_config.h"
#include "libmesh/getpot.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/enum_solver_package.h"

#include "libmesh/serial_mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/mesh_refinement.h"

#include "libmesh/dof_map.h"
#include "libmesh/linear_solver.h"
#include "libmesh/equation_systems.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/periodic_boundary.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/slepc_macro.h"
#include "libmesh/perf_log.h"

// User defined includes
#include "point_particle.h"
#include "particle_mesh.h"
#include "point_mesh.h"
#include "pm_system_stokes.h"
#include "brownian_system.h"
#include "pm_periodic_boundary.h"
#include "chebyshev.h"
#include "pm_toolbox.h"
#include "polymer_chain.h"
#include "random_generator.h"
#include "stokes_solver.h"
#include "ggem_system.h"
#include "copss_init.h"
#include "fix/fix_factory.h"
#include "fix/fix.h"



// Bring in everything from the libMesh namespace


namespace libMesh{

class Copss
{
public: 
  // PETSC MPI communicator
  Parallel::Communicator *comm_in;
  // error message string
  std::string error_msg;
  // output message string
  std::string output_msg;
  // control file object
  GetPot input_file;
  // control file name
  std::string control_fileName;
  // test name
  std::string test_name;
  bool print_info;
  // physical parameters
  const Real kB = 1.380662E-17;//1.380662E-23(J/K) = 1.3806623E-23(N*m/K) = 1.380662E-17(N*um/K)
  const Real PI = libMesh::pi;
  Real T; // simulation temperature (K)
  Real kBT;// (N*um)
  Real viscosity; // viscosity of the fluid (cp = N*s/um^2)
  Real Rb; // radius of the bead
  Real drag_c; // Drag coefficient (N*s/um)
  Real Db;
  std::string particle_type;

  // characteristic variables
  Real tc; // characteristic time (diffusion time) (s)
  Real uc; // characteristic velocity (um/s)
  Real fc; // characteristic force (N)
  Real muc; // non-dimension viscosity

  // Geometry information
  unsigned int dim; // dimension of the box
  std::string wall_type; // wall_type (slit or sphere)
  std::vector<Real> wall_params; // (wall parameters; e.g. slit = '-50,50,-50,50,-50,50' or sphere = '10')
  std::vector<bool> periodicity; // periodicity of the box
  std::vector<bool> inlet; // inlet direction of the box
  std::vector<Real> inlet_pressure; // inlet pressure
  std::vector<bool> shear; // apply shear or not on boundary pairs
  std::vector<Real> shear_rate; // shear rate on boundary pairs
  std::vector<unsigned int> shear_direction; // shear direction on boundary pairs

  // Mesh information
  bool generate_mesh; // flag to generate mesh or load mesh
  std::string domain_mesh_file; // domain mesh filename
  std::vector<unsigned int> n_mesh; // mesh size in all directions

  // Fix
  std::vector<Fix*> fixes;
  FixFactory* fix_factory;
  unsigned int numForceTypes;
  std::vector<std::string> forceTypes;
  std::vector<Fix::type_force> forces;
  // map force_type to Fix Pointers


  // GGEM information
  Real alpha;

  // Solver information
  int max_linear_iterations;
  Real linear_solver_rtol;
  Real linear_solver_atol;
  bool user_defined_pc;
  bool schur_user_ksp;
  Real schur_user_ksp_rtol;
  Real schur_user_ksp_atol;
  std::string schur_pc_type;
  std::string stokes_solver_type;
  StokesSolverType solver_type; 

  // Chebyshev information
  unsigned int max_n_cheb; // max order of Chebyshev polynomianl 
  Real tol_cheb; // tolerance of chebyshev convergence
  Real eig_factor; // factor to resize eigen value range
  Real tol_eigen; // tolerance of eigenvalue convergence
  bool compute_eigen; // if compute eigen value; false by default
  bool read_eigen; // if read eigen value from existed file; false by default
  // Run time info
  bool with_hi;
  bool with_brownian; // if consider brownian motion
  std::size_t random_seed;
  bool adaptive_dt; // if use adaptive time step (essential for brownian systems)
  std::vector<Real> max_dr_coeff; // max displacement per step
  Real max_dr;
  bool restart; // if restart
  std::size_t restart_step; // restart step
  unsigned int nstep; // totol number of steps to run
  unsigned int n_relax_step; // relaxation steps when Chebyshev cannot converge even after recompute eigenvalue
  bool debug_info;
  // ouput file 
  unsigned int write_interval; // output file write interval
  std::vector<std::string> output_file;
  // mesh
  SerialMesh* mesh;
  PointMesh<3>* point_mesh;
  //std::unique_ptr<SerialMesh> mesh;
  Real min_mesh_size, max_mesh_size;
  Real search_radius_p, search_radius_e;

  //periodic boundary
  PMPeriodicBoundary* pm_periodic_boundary;
  //std::unique_ptr<PMPeriodicBoundary> pm_periodic_boundary;

  // equation system
  unsigned int u_var, v_var, w_var, p_var;

  //integrate
  // paramters for dynamic process;
  bool reinit_stokes;
  unsigned int NP;
  unsigned int n_vec;
  Real hminf, hmaxf; // fluid mesh minimum
  Real hmin, hmax; // all mesh minimum
  bool cheb_converge = false;
  unsigned int n_chebyshev_failure = 0;
  Real eig_min = 0, eig_max = 0;
  const std::string out_system_filename = "output_pm_system.e";
  UniquePtr<NumericVector<Real>> v0_ptr;
  ExodusII_IO* exodus_ptr;

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   Create vectors and Shell Mat for use:
   U0:          particle velocity vector;
   R0/R_mid:    particle position vector;
   dw/dw_mid:   random vector;
   RIN/ROUT:    the initial and intermediate particle postion vector for msd output
   RIN will not change, and ROUT excludes pbc
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  Vec             U0, R0, R_mid, RIN,ROUT, dw, dw_mid;
  PetscViewer     viewer;
  Mat             M;
  PetscRandom     rand_ctx;
  PetscScalar     coef = 0.0;
  BrownianSystem* brownian_sys; 
  std::vector<Point> center0;

  //variables in integrator
  unsigned int istart;
  unsigned int o_step;
  Real real_time; // real time
  bool update_neighbor_list_everyStep;
  unsigned int neighbor_list_update_interval = 0; // how often to update neighbor list
  unsigned int timestep_duration; // how many steps elapsed since last neighbor_list update
  bool neighbor_list_update_flag;
  std::vector<Real> vel0;
  std::vector<Real> vel1;

  





  Copss(CopssInit& init);

  ~Copss();
  
  int check_libmesh();

  void start_time(struct tm * timeinfo);

  void end_time(struct tm * timeinfo);

  EquationSystems init_system(std::string input_file); // including the following 10 steps

  void read_input();

  
  virtual void create_object_mesh() = 0;

  EquationSystems create_equation_systems();

  virtual void run(EquationSystems& equation_systems) = 0;


  void destroy(); 


protected:
  void read_system_info(); 
  void read_physical_info(); 
  virtual void read_particle_info() = 0;
  void read_domain_info(); 
  void read_force_info();
  virtual void read_ggem_info() = 0;
  void read_stokes_solver_info();
  void read_chebyshev_info();
  void read_run_info();
  void read_restart_time();
  void read_restart_eigenvalue();

  void create_domain_mesh();
  void create_periodic_boundary();
  virtual void create_object() = 0;

  virtual void attach_object_mesh(PMLinearImplicitSystem& system) = 0;
  void attach_period_boundary(PMLinearImplicitSystem& system);
  virtual void set_parameters(EquationSystems& equation_systems) = 0;

  // force field
  void attach_fixes(PMLinearImplicitSystem& system);

  void solve_undisturbed_system(EquationSystems& equation_systems);

  void create_brownian_system(EquationSystems& equation_systems);

  void fixman_integrate(EquationSystems& equation_systems, unsigned int& i);

  void langevin_integrate(EquationSystems& equation_systems, unsigned int& i);

  // update object positions due to PBS
  virtual void update_object() {};
  virtual void write_object(unsigned int step_id) = 0;

  // output precision
  const int o_precision = 6;
};

} // end namespace libMesh