Thm seg parallel

CMakeLists.txt

@@ -92,7 +92,6 @@ DEAL_II_INITIALIZE_CACHED_VARIABLES()
 PROJECT(${TARGET})
 DEAL_II_INVOKE_AUTOPILOT()
 
-
 #add_library(lthm SHARED ${thm_src})
 
 # add_executable(geothm ${thm_SOURCE_DIR}/src/main.cc)

README.md

@@ -40,6 +40,8 @@ Please refer to [CB-Geo thm Documentation](https://cb-geo.github.io/thm-doc) for
 #login TACC
 ssh [email protected]
 
+nano ~/.bashrc
+# Write the following commands into the .bashrc file
 module load impi-largemem
 export CC=icc
 export CXX=icpc
@@ -66,10 +68,10 @@ wget http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/metis-5.1.0.tar.gz && \
 module load boost
 
 # Clone and compile dealii
-cd $HOME
+cds
 git clone https://github.com/dealii/dealii --depth=1 dealii-src
 cd dealii-src/ && mkdir build &&  cd build
-cmake -DCMAKE_INSTALL_PREFIX=$HOME/dealii -DPETSC_DIR=$PETSC_DIR -DPETSC_ARCH=clx -DDEAL_II_WITH_PETSC=On  -DDEAL_II_WITH_METIS=On -DMETIS_DIR=$HOME/metis/ -DDEAL_II_WITH_MPI=On ..
+cmake -DCMAKE_INSTALL_PREFIX=$SCRATCH/dealii -DPETSC_DIR=$PETSC_DIR -DPETSC_ARCH=clx -DDEAL_II_WITH_PETSC=On  -DDEAL_II_WITH_METIS=On -DMETIS_DIR=$HOME/metis/ -DDEAL_II_WITH_MPI=On ..
 make install -j4
 
 
@@ -79,7 +81,7 @@ git clone https://github.com/cb-geo/thm
 cd thm
 git checkout thm_seg_parallel
 mkdir build && cd build
-cmake -DCMAKE_BUILD_TYPE=Release -DDEAL_II_DIR=$HOME/dealii/ ..
+cmake -DCMAKE_BUILD_TYPE=Release -DDEAL_II_DIR=$SCRATCH/dealii/ ..
 make -j
 ```
 

include/boundaryvalues.h

@@ -17,13 +17,14 @@ class PressureDirichletBoundaryValues : public Function<dim> {
  private:
   const double period_;  // value
   int boundary_id_{-1};
+  int bd_i_;
 
  public:
   PressureDirichletBoundaryValues()
-      : Function<dim>(), period_(0.2) {}  // 之前的没有
+      : Function<dim>(), period_(0.2), bd_i_(0)  {}  // 之前的没有
   virtual double value(const Point<dim>& p,
                        const unsigned int component = 0) const;  // boundary
-
+  virtual void get_bd_i(int bd_i) {bd_i_ = bd_i;}
   // virtual void vector_value(const Point<dim>& p,  //放在这里没啥用
   //                           Vector<double>& value) const;
   virtual void set_boundary_id(int bnd_id) { boundary_id_ = bnd_id; }
@@ -56,10 +57,12 @@ class TemperatureDirichletBoundaryValues : public Function<dim> {
  private:
   const double period_;  // value
   int boundary_id_{-1};
+  int bd_i_;
 
  public:
   TemperatureDirichletBoundaryValues()
-      : Function<dim>(), period_(0.2) {}  // 之前的没有
+      : Function<dim>(), period_(0.2) ,bd_i_(0){}  // 之前的没有
+  virtual void get_bd_i(int bd_i) {bd_i_ = bd_i;}
   virtual double value(const Point<dim>& p,
                        const unsigned int component = 0) const;  // boundary
   // virtual void vector_value(const Point<dim>& p,  //放在这里没啥用
@@ -88,7 +91,8 @@ double TemperatureDirichletBoundaryValues<dim>::value(
   // }
   if (boundary_id_ == 4) {
     return g_Tb_well;
-  } else if (boundary_id_ == 3 || boundary_id_ == 8 || boundary_id_ == 12 ||
+  } 
+  if (boundary_id_ == 3 || boundary_id_ == 8 || boundary_id_ == 12 ||
              boundary_id_ == 13 || boundary_id_ == 14) {
     return g_Tb_top + g_T_grad * (0. - p[2]);
   }
@@ -106,10 +110,12 @@ class PressureNeumanBoundaryValues : public Function<dim> {
  private:
   const double period_;  // value
   int boundary_id_{-1};
+  int bd_i_;
 
  public:
   PressureNeumanBoundaryValues()
-      : Function<dim>(), period_(0.2) {}  // 之前的没有
+      : Function<dim>(), period_(0.2), bd_i_(0){}  // 之前的没有
+  virtual void get_bd_i(int bd_i) {bd_i_ = bd_i;}
   virtual double value(const Point<dim>& p,
                        const unsigned int component = 0) const;  // boundary
 
@@ -136,10 +142,12 @@ class TemperatureNeumanBoundaryValues : public Function<dim> {
  private:
   const double period;  // value
   int boundary_id_{-1};
+  int bd_i_;
 
  public:
   TemperatureNeumanBoundaryValues()
-      : Function<dim>(), period(0.2) {}  // 之前的没有
+      : Function<dim>(), period(0.2), bd_i_(0)  {}  // 之前的没有
+  virtual void get_bd_i(int bd_i) {bd_i_ = bd_i;}
   virtual double value(const Point<dim>& p,
                        const unsigned int component = 0) const;  // boundary
   // virtual void vector_value(const Point<dim>& p,  //放在这里没啥用
@@ -159,9 +167,11 @@ double TemperatureNeumanBoundaryValues<dim>::value(
   if (boundary_id_ == 4) {
     return g_QT_well;  // boundary value is set to zero in
                        // this case
-  } else if (boundary_id_ == 3) {
+  } 
+  if (boundary_id_ == 8) {
     return g_QT_top;
-  } else if (boundary_id_ == 8) {
+  } 
+  if (boundary_id_ == 3) {
     return g_QT_bottom;
   }
 }