pi-DoMUS: Parallel Deal.II MUltiphysics Solver
hydrogels_three_fields.h
Go to the documentation of this file.
1 
5 #ifndef _hydrogels_three_fields_h_
6 #define _hydrogels_three_fields_h_
7 
8 #include "pde_system_interface.h"
9 #include <deal2lkit/parsed_function.h>
10 
11 
12 #include <deal.II/fe/fe_values.h>
17 
21 
22 #include <deal.II/lac/solver_cg.h>
24 
25 #include "lac/lac_type.h"
26 
27 
28 
29 
30 template <int dim, int spacedim, typename LAC>
31 class HydroGelThreeFields : public PDESystemInterface<dim,spacedim, HydroGelThreeFields<dim,spacedim,LAC>, LAC>
32 {
33 public:
34 
36 
38 
41 
42  void set_matrix_couplings(std::vector<std::string> &couplings) const;
43 
44  template <typename EnergyType, typename ResidualType>
46  FEValuesCache<dim,spacedim> &scratch,
47  std::vector<EnergyType> &energies,
48  std::vector<std::vector<ResidualType> > &local_residuals,
49  bool compute_only_system_terms) const;
50 
51 
53  const std::vector<shared_ptr<LATrilinos::BlockMatrix> >,
56 
57 
58 
59 private:
60  double T;
61  double Omega;
62  double G;
63  double chi;
64  double l0;
65 
66  double mu;
67  double mu0;
68  double l02;
69  double l03;
70  double l0_3;
71  double l0_6;
72  const double R=8.314;
73 
74 
75  mutable shared_ptr<TrilinosWrappers::PreconditionAMG> U_prec;
76  mutable shared_ptr<TrilinosWrappers::PreconditionJacobi> p_prec;
77  mutable shared_ptr<TrilinosWrappers::PreconditionJacobi> c_prec;
78 
79 };
80 
81 template <int dim, int spacedim, typename LAC>
83  PDESystemInterface<dim,spacedim,HydroGelThreeFields<dim,spacedim,LAC>, LAC>("Free Swelling Three Fields",
84  dim+2,2,
85  "FESystem[FE_Q(2)^d-FE_Q(1)-FE_Q(2)]",
86  "u,u,u,p,c","1,0,0")
87 {}
88 
89 template <int dim, int spacedim, typename LAC>
90 void
92 set_matrix_couplings(std::vector<std::string> &couplings) const
93 {
94  couplings[0] = "1,1,0;1,0,1;0,1,1";
95  couplings[1] = "0,0,0;0,1,0;0,0,0";
96 }
97 
98 template <int dim, int spacedim, typename LAC>
99 template <typename EnergyType, typename ResidualType>
100 void
103  FEValuesCache<dim,spacedim> &fe_cache,
104  std::vector<EnergyType> &energies,
105  std::vector<std::vector<ResidualType> > &,
106  bool compute_only_system_terms) const
107 {
108  EnergyType alpha = this->alpha;
109  this->reinit(alpha, cell, fe_cache);
110 
111  const FEValuesExtractors::Vector displacement(0);
112  const FEValuesExtractors::Scalar pressure(dim);
113  const FEValuesExtractors::Scalar concentration(dim+1);
114 
115  auto &Fs = fe_cache.get_deformation_gradients("solution", "Fu", displacement, alpha);
116 
117  auto &ps = fe_cache.get_values("solution", "p", pressure, alpha);
118 
119  auto &cs = fe_cache.get_values("solution", "c", concentration, alpha);
120 
121  const unsigned int n_q_points = ps.size();
122 
123  auto &JxW = fe_cache.get_JxW_values();
124 
125  for (unsigned int q=0; q<n_q_points; ++q)
126  {
127  const Tensor<2, dim, EnergyType> &F = Fs[q];
128  const Tensor<2, dim, EnergyType> C = transpose(F)*F;
129  const EnergyType &c = cs[q];
130  const EnergyType &p = ps[q];
131 
132  const EnergyType I = trace(C);
133  const EnergyType J = determinant(F);
134 
135 
136  EnergyType psi = ( 0.5*G*l0_3*(l02*I - dim)
137 
138  + (l0_3*R*T/Omega)*((Omega*l03*c)*std::log((Omega*l03*c)/(1.+Omega*l03*c))
139  + chi*((Omega*l03*c)/(1.+Omega*l03*c)) )
140 
141  - (mu0)*c - p*(J-l0_3-Omega*c)) ;
142 
143  energies[0] += psi*JxW[q];
144 
145  if (!compute_only_system_terms)
146  energies[1] += 0.5*(p*p)*JxW[q];
147  }
148 
149 }
150 
151 template <int dim, int spacedim, typename LAC>
153 {
155  this->add_parameter(prm, &T, "T", "298.0", Patterns::Double(0.0));
156  this->add_parameter(prm, &Omega, "Omega", "1e-5", Patterns::Double(0.0));
157  this->add_parameter(prm, &chi, "chi", "0.1", Patterns::Double(0.0));
158  this->add_parameter(prm, &l0, "l0", "1.5", Patterns::Double(1.0));
159  this->add_parameter(prm, &G, "G", "10e3", Patterns::Double(0.0));
160  this->add_parameter(prm, &mu, "mu", "-1.8e-22", Patterns::Double(-1000.0));
161 }
162 
163 template <int dim, int spacedim, typename LAC>
165 {
166  l02 = l0*l0;
167  l03 = l02*l0;
168  l0_3 = 1./l03;
169  l0_6 = 1./(l03*l03);
170 
171  mu0 = R*T*(std::log((l03-1.)/l03) + l0_3 + chi*l0_6) + G*Omega/l0;
172 }
173 
174 
175 template <int dim, int spacedim, typename LAC>
176 void
179  const std::vector<shared_ptr<LATrilinos::BlockMatrix> > matrices,
182 {
183 
184  std::vector<std::vector<bool> > constant_modes;
185  FEValuesExtractors::Vector displacement(0);
186  DoFTools::extract_constant_modes (dh, dh.get_fe().component_mask(displacement),
187  constant_modes);
188 
189  p_prec.reset (new TrilinosWrappers::PreconditionJacobi());
190  c_prec.reset (new TrilinosWrappers::PreconditionJacobi());
191  U_prec.reset (new TrilinosWrappers::PreconditionAMG());
192 
194  Amg_data.constant_modes = constant_modes;
195  Amg_data.elliptic = true;
196  Amg_data.higher_order_elements = true;
197  Amg_data.smoother_sweeps = 2;
198  Amg_data.aggregation_threshold = 0.02;
199 
200 
201  U_prec->initialize (matrices[0]->block(0,0), Amg_data);
202  p_prec->initialize (matrices[1]->block(1,1));
203  c_prec->initialize (matrices[0]->block(2,2));
204 
205 
206  // SYSTEM MATRIX:
207  auto A = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(0,0) );
208  auto Bt = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(0,1) );
209  auto Z02 = 0*linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(0,2) );
210 
211  auto B = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(1,0) );
212  auto Z11 = 0*linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(1,1) );
213  auto C = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(1,2) );
214 
215  auto Z20 = 0*linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(2,0) );
216  auto D = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(2,1) );
217  auto E = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(2,2) );
218 
219  auto P0 = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(0,0));
220  auto P1 = linear_operator< LATrilinos::VectorType::BlockType >( matrices[1]->block(1,1));
221  auto P2 = linear_operator< LATrilinos::VectorType::BlockType >( matrices[0]->block(2,2));
222 
223 
224  static ReductionControl solver_control_pre(5000, 1e-3);
225  static SolverCG<LATrilinos::VectorType::BlockType> solver_CG(solver_control_pre);
226 
227  auto P0_inv = inverse_operator( P0, solver_CG, *U_prec);
228  auto P1_inv = inverse_operator( P1, solver_CG, *p_prec);
229  auto P2_inv = inverse_operator( P2, solver_CG, *c_prec);
230 
231  auto P0_i = P0_inv;
232  auto P1_i = P1_inv;
233  auto P2_i = P2_inv;
234 
235 
236  const std::array<std::array<LinearOperator<LATrilinos::VectorType::BlockType>, 3 >, 3 > matrix_array = {{
237  {{ A, Bt , Z02 }},
238  {{ B, Z11 , C }},
239  {{ Z20, D , E }}
240  }
241  };
242 
243  system_op = block_operator<3, 3, LATrilinos::VectorType >(matrix_array);
244 
245  const std::array<LinearOperator<LATrilinos::VectorType::BlockType>, 3 > diagonal_array = {{ P0_i, P1_i, P2_i }};
246 
247 
248  prec_op = block_diagonal_operator<3,LATrilinos::VectorType>(diagonal_array);
249 
250 }
251 
252 
253 
254 #endif
255 
void declare_parameters(ParameterHandler &prm)
Definition: hydrogels_three_fields.h:152
void compute_system_operators(const DoFHandler< dim, spacedim > &, const std::vector< shared_ptr< LATrilinos::BlockMatrix > >, LinearOperator< LATrilinos::VectorType > &, LinearOperator< LATrilinos::VectorType > &) const
Definition: hydrogels_three_fields.h:178
Definition: lac_type.h:33
std::vector< std::vector< bool > > constant_modes
ActiveSelector::active_cell_iterator active_cell_iterator
This is the class that users should derive from.
Definition: pde_system_interface.h:46
void set_matrix_couplings(std::vector< std::string > &couplings) const
Definition: hydrogels_three_fields.h:92
void parse_parameters_call_back()
Definition: hydrogels_three_fields.h:164
LinearOperator< typename Solver::vector_type, typename Solver::vector_type > inverse_operator(const LinearOperator< typename Solver::vector_type, typename Solver::vector_type > &op, Solver &solver, const Preconditioner &preconditioner)
~HydroGelThreeFields()
Definition: hydrogels_three_fields.h:35
Definition: hydrogels_three_fields.h:31
void energies_and_residuals(const typename DoFHandler< dim, spacedim >::active_cell_iterator &cell, FEValuesCache< dim, spacedim > &scratch, std::vector< EnergyType > &energies, std::vector< std::vector< ResidualType > > &local_residuals, bool compute_only_system_terms) const
Definition: hydrogels_three_fields.h:102
HydroGelThreeFields()
Definition: hydrogels_three_fields.h:82
void extract_constant_modes(const DoFHandlerType &dof_handler, const ComponentMask &component_mask, std::vector< std::vector< bool > > &constant_modes)
const FiniteElement< dim, spacedim > & get_fe() const