Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
Crystal_Growth_Phase_Field_Model.h
Go to the documentation of this file.
1
118 *  
119 *   #include "PhaseFieldSolver.h"
120 *  
121 *   void InitialValues::vector_value(const Point<2> &p,
122 *   Vector<double> & values) const
123 *   {
124 *   values(0)= 0.0; //Initial p value of domain
125 *   values(1)= 0.2; //Initial temperature of domain
126 *   }
127 * @endcode
128
129
130<a name="ann-PhaseFieldSolver.cpp"></a>
131<h1>Annotated version of PhaseFieldSolver.cpp</h1>
132 *
133 *
134 *
135 *
136 * @code
137 *   /* -----------------------------------------------------------------------------
138 *   *
139 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
140 *   * Copyright (C) 2024 by Umair Hussain
141 *   *
142 *   * This file is part of the deal.II code gallery.
143 *   *
144 *   * -----------------------------------------------------------------------------
145 *   */
146 *  
147 *   #include "PhaseFieldSolver.h"
148 *  
149 *   PhaseFieldSolver::PhaseFieldSolver()
150 *   : mpi_communicator(MPI_COMM_WORLD)
151 *   , n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator))
152 *   , this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator))
153 *   , pcout(std::cout, (this_mpi_process == 0))
154 *   , fe(FE_Q<2>(1), 2)
155 *   , dof_handler(triangulation)
156 *   , time(0.0)
157 *   , final_time(1.)
158 *   , time_step(.0002)
159 *   , theta(0.5)
160 *   , epsilon(0.01)
161 *   , tau(0.0003)
162 *   , gamma(10.)
163 *   , latent_heat(1.4)
164 *   , alpha(0.9)
165 *   , t_eq(1.)
166 *   , a(0.01)
167 *   {}
168 * @endcode
169
170
171<a name="ann-PhaseFieldSolver.h"></a>
172<h1>Annotated version of PhaseFieldSolver.h</h1>
173 *
174 *
175 *
176 *
177
178 *
179 *
180 * @code
181 *   /* -----------------------------------------------------------------------------
182 *   *
183 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
184 *   * Copyright (C) 2024 by Umair Hussain
185 *   *
186 *   * This file is part of the deal.II code gallery.
187 *   *
188 *   * -----------------------------------------------------------------------------
189 *   */
190 *  
191 *   #ifndef KOBAYASHI_PARALLEL_PHASEFIELDSOLVER_H
192 *   #define KOBAYASHI_PARALLEL_PHASEFIELDSOLVER_H
193 *  
194 *   #include <deal.II/base/quadrature_lib.h>
195 *   #include <deal.II/base/function.h>
196 *   #include <deal.II/base/utilities.h>
197 *   #include <deal.II/lac/vector.h>
198 *   #include <deal.II/lac/full_matrix.h>
199 *   #include <deal.II/lac/sparse_matrix.h>
200 *   #include <deal.II/lac/sparse_direct.h>
201 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
202 *   #include <deal.II/lac/solver_cg.h>
203 *   #include <deal.II/lac/precondition.h>
204 *   #include <deal.II/lac/affine_constraints.h>
205 *   #include <deal.II/grid/tria.h>
206 *   #include <deal.II/grid/grid_generator.h>
207 *   #include <deal.II/dofs/dof_handler.h>
208 *   #include <deal.II/dofs/dof_tools.h>
209 *   #include <deal.II/fe/fe_q.h>
210 *   #include <deal.II/fe/fe_values.h>
211 *   #include <deal.II/fe/fe_system.h>
212 *   #include <deal.II/numerics/vector_tools.h>
213 *   #include <deal.II/numerics/matrix_tools.h>
214 *   #include <deal.II/numerics/data_out.h>
215 *   #include <deal.II/grid/grid_in.h>
216 *  
217 * @endcode
218 *
219 * For Parallel Computation
220 *
221 * @code
222 *   #include <deal.II/base/conditional_ostream.h>
223 *   #include <deal.II/base/mpi.h>
224 *   #include <deal.II/lac/petsc_vector.h>
225 *   #include <deal.II/lac/petsc_sparse_matrix.h>
226 *   #include <deal.II/lac/petsc_solver.h>
227 *   #include <deal.II/lac/petsc_precondition.h>
228 *   #include <deal.II/grid/grid_tools.h>
229 *   #include <deal.II/dofs/dof_renumbering.h>
230 *  
231 *   #include <deal.II/distributed/solution_transfer.h>
232 *  
233 *   #include <fstream>
234 *   #include <iostream>
235 *  
236 *   using namespace dealii;
237 *  
238 *   class PhaseFieldSolver {
239 *   public:
240 *   PhaseFieldSolver();
241 *   void run();
242 *  
243 *   private:
244 *   void make_grid_and_dofs();
245 *   void assemble_system();
246 *   void solve();
247 *   void output_results(const unsigned int timestep_number) const;
248 *   double compute_residual();
249 *   void applying_bc();
250 *   float get_random_number();
251 *  
252 *   MPI_Comm mpi_communicator;
253 *   const unsigned int n_mpi_processes;
254 *   const unsigned int this_mpi_process;
255 *   ConditionalOStream pcout;
256 *  
258 *   FESystem<2> fe;
259 *   DoFHandler<2> dof_handler;
260 *   GridIn<2> gridin;
261 *  
262 *   PETScWrappers::MPI::SparseMatrix jacobian_matrix;
263 *  
264 *   double time;
265 *   const double final_time, time_step;
266 *   const double theta;
267 *   const double epsilon, tau, gamma, latent_heat, alpha, t_eq, a; //as given in Ref. [1]
268 *  
269 *   PETScWrappers::MPI::Vector conv_solution; //solution vector at last newton-raphson iteration
270 *   PETScWrappers::MPI::Vector old_solution; //solution vector at last time step
271 *   PETScWrappers::MPI::Vector solution_update; //increment in solution or delta solution
272 *   PETScWrappers::MPI::Vector system_rhs; //to store residual
273 *   Vector<double> conv_solution_np, old_solution_np; //creating non parallel vectors to store data for easy access of old solution values by all processes
274 *  
275 *   };
276 *  
277 * @endcode
278 *
279 * Initial values class
280 *
281 * @code
282 *   class InitialValues : public Function<2>
283 *   {
284 *   public:
285 *   InitialValues(): Function<2>(2)
286 *   {}
287 *   virtual void vector_value(const Point<2> &p,
288 *   Vector<double> & value) const override;
289 *   };
290 *  
291 *  
292 *   #endif //KOBAYASHI_PARALLEL_PHASEFIELDSOLVER_H
293 * @endcode
294
295
296<a name="ann-applying_bc.cpp"></a>
297<h1>Annotated version of applying_bc.cpp</h1>
298 *
299 *
300 *
301 *
302
303 *
304 *
305 * @code
306 *   /* -----------------------------------------------------------------------------
307 *   *
308 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
309 *   * Copyright (C) 2024 by Umair Hussain
310 *   *
311 *   * This file is part of the deal.II code gallery.
312 *   *
313 *   * -----------------------------------------------------------------------------
314 *   */
315 *  
316 *   #include "PhaseFieldSolver.h"
317 *  
318 *   void PhaseFieldSolver::applying_bc(){
319 *   FEValuesExtractors::Scalar phase_parameter(0);
320 *   FEValuesExtractors::Scalar temperature(1);
321 *  
322 *   QGauss<2> quadrature_formula(fe.degree + 1);
323 *   FEValues<2> fe_values(fe,
324 *   quadrature_formula,
326 *  
327 *   ComponentMask p_mask = fe.component_mask(phase_parameter);
328 *   ComponentMask t_mask = fe.component_mask(temperature);
329 *  
330 *   std::map<types::global_dof_index,double> boundary_values;
331 *  
332 * @endcode
333 *
334 * Prescribing p=1 at the left face (this will be maintained in the subsequent iterations when zero BC is applied in the Newton-Raphson iterations)
335 *
336 * @code
338 *   1,
340 *   boundary_values,p_mask);
341 *  
342 * @endcode
343 *
344 * To apply the boundary values only to the solution vector without the Jacobian Matrix and RHS Vector
345 *
346 * @code
347 *   for (auto &boundary_value : boundary_values)
348 *   old_solution(boundary_value.first) = boundary_value.second;
349 *  
350 *   old_solution.compress(VectorOperation::insert);
351 *  
352 *   }
353 * @endcode
354
355
356<a name="ann-assemble_system.cpp"></a>
357<h1>Annotated version of assemble_system.cpp</h1>
358 *
359 *
360 *
361 *
362 * @code
363 *   /* -----------------------------------------------------------------------------
364 *   *
365 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
366 *   * Copyright (C) 2024 by Umair Hussain
367 *   *
368 *   * This file is part of the deal.II code gallery.
369 *   *
370 *   * -----------------------------------------------------------------------------
371 *   */
372 *  
373 *   #include "PhaseFieldSolver.h"
374 *   #include <cmath>
375 *  
376 *   void PhaseFieldSolver::assemble_system() {
377 * @endcode
378 *
379 * Separating each variable as a scalar to easily call the respective shape functions
380 *
381 * @code
382 *   FEValuesExtractors::Scalar phase_parameter(0);
383 *   FEValuesExtractors::Scalar temperature(1);
384 *  
385 *   QGauss<2> quadrature_formula(fe.degree + 1);
386 *   FEValues<2> fe_values(fe,
387 *   quadrature_formula,
389 *  
390 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
391 *   const unsigned int n_q_points = quadrature_formula.size();
392 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
393 *  
394 *   Vector<double> cell_rhs(dofs_per_cell);
395 *  
396 * @endcode
397 *
398 * To copy values and gradients of solution from previous iteration
399 * Old Newton iteration
400 *
401 * @code
402 *   std::vector<Tensor<1, 2>> old_newton_solution_gradients_p(n_q_points);
403 *   std::vector<double> old_newton_solution_values_p(n_q_points);
404 *   std::vector<Tensor<1, 2>> old_newton_solution_gradients_t(n_q_points);
405 *   std::vector<double> old_newton_solution_values_t(n_q_points);
406 * @endcode
407 *
408 * Old time step iteration
409 *
410 * @code
411 *   std::vector<Tensor<1, 2>> old_time_solution_gradients_p(n_q_points);
412 *   std::vector<double> old_time_solution_values_p(n_q_points);
413 *   std::vector<Tensor<1, 2>> old_time_solution_gradients_t(n_q_points);
414 *   std::vector<double> old_time_solution_values_t(n_q_points);
415 *  
416 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
417 *   jacobian_matrix.operator=(0.0);
418 *   system_rhs.operator=(0.0);
419 *  
420 *   for (const auto &cell : dof_handler.active_cell_iterators()){
421 *   if (cell->subdomain_id() == this_mpi_process) {
422 *   cell_matrix = 0;
423 *   cell_rhs = 0;
424 *  
425 *   fe_values.reinit(cell);
426 *  
427 * @endcode
428 *
429 * Copying old solution values
430 *
431 * @code
432 *   fe_values[phase_parameter].get_function_values(conv_solution_np,old_newton_solution_values_p);
433 *   fe_values[phase_parameter].get_function_gradients(conv_solution_np,old_newton_solution_gradients_p);
434 *   fe_values[temperature].get_function_values(conv_solution_np,old_newton_solution_values_t);
435 *   fe_values[temperature].get_function_gradients(conv_solution_np,old_newton_solution_gradients_t);
436 *   fe_values[phase_parameter].get_function_values(old_solution_np,old_time_solution_values_p);
437 *   fe_values[phase_parameter].get_function_gradients(old_solution_np,old_time_solution_gradients_p);
438 *   fe_values[temperature].get_function_values(old_solution_np,old_time_solution_values_t);
439 *   fe_values[temperature].get_function_gradients(old_solution_np,old_time_solution_gradients_t);
440 *  
441 *   for (unsigned int q = 0; q < n_q_points; ++q){
442 *   double khi = get_random_number();
443 * @endcode
444 *
445 * Old solution values
446 *
447 * @code
448 *   double p_on = old_newton_solution_values_p[q]; //old newton solution
449 *   auto grad_p_on = old_newton_solution_gradients_p[q];
450 *   double p_ot = old_time_solution_values_p[q]; //old time step solution
451 *   auto grad_p_ot = old_time_solution_gradients_p[q];
452 *   double t_on = old_newton_solution_values_t[q];
453 *   auto grad_t_on = old_newton_solution_gradients_t[q];
454 *   double t_ot = old_time_solution_values_t[q];
455 *   auto grad_t_ot = old_time_solution_gradients_t[q];
456 *   for (unsigned int i = 0; i < dofs_per_cell; ++i){
457 * @endcode
458 *
459 * Shape Functions
460 *
461 * @code
462 *   double psi_i = fe_values[phase_parameter].value(i,q);
463 *   auto grad_psi_i = fe_values[phase_parameter].gradient(i,q);
464 *   double phi_i = fe_values[temperature].value(i,q);
465 *   auto grad_phi_i = fe_values[temperature].gradient(i,q);
466 *   for (unsigned int j = 0; j < dofs_per_cell; ++j){
467 * @endcode
468 *
469 * Shape Functions
470 *
471 * @code
472 *   double psi_j = fe_values[phase_parameter].value(j,q);
473 *   auto grad_psi_j = fe_values[phase_parameter].gradient(j,q);
474 *   double phi_j = fe_values[temperature].value(j,q);
475 *   auto grad_phi_j = fe_values[temperature].gradient(j,q);
476 *  
477 *   double mp = psi_i*(tau*psi_j);
478 *   double kp = grad_psi_i*(std::pow(epsilon,2)*grad_psi_j);
479 *   double m = (alpha/M_PI)*std::atan(gamma*(t_eq - t_on));
480 *   double t1 = (1-p_on)*(p_on-0.5+m);
481 *   double t2 = -(p_on)*(p_on-0.5+m);
482 *   double t3 = (p_on)*(1-p_on);
483 *   double nl_p = psi_i*((t1+t2+t3)*psi_j);
484 * @endcode
485 *
486 * Adding random noise at the interface
487 *
488 * @code
489 *   nl_p -= a*khi*psi_i*((1.0 - 2*(p_on))*psi_j);
490 *   double f1_p= mp + time_step*theta*kp - time_step*theta*nl_p; // doh f1 by doh p (first Jacobian terms)
491 *  
492 *   double t4 = (p_on)*(1-p_on)*(-(alpha*gamma/(M_PI*(1+std::pow((gamma*(t_eq-t_on)),2)))));
493 *   double nl_t = psi_i*(t4*phi_j);
494 *   double f1_t = -time_step*theta*nl_t; // doh f1 by doh t (second Jacobian terms)
495 *  
496 *   double mpt = phi_i*(latent_heat*psi_j);
497 *   double f2_p = -mpt; // doh f2 by doh p (third Jacobian terms)
498 *  
499 *   double mt = phi_i*(phi_j);
500 *   double kt = grad_phi_i*(grad_phi_j);
501 *   double f2_t = mt + time_step*theta*kt; // doh f2 by doh t (fourth Jacobian terms)
502 *  
503 * @endcode
504 *
505 * Assembling Jacobian matrix
506 *
507 * @code
508 *   cell_matrix(i,j) += (f1_p + f1_t + f2_p + f2_t)*fe_values.JxW(q);
509 *  
510 *   }
511 * @endcode
512 *
513 * Finding f1 and f2 at previous iteration for rhs vector
514 *
515 * @code
516 *   double mp_n = psi_i*(tau*p_on);
517 *   double kp_n = grad_psi_i*(std::pow(epsilon,2)*grad_p_on);
518 *   double m_n = (alpha/M_PI)*std::atan(gamma*(t_eq-t_on));
519 *   double nl_n = psi_i*((p_on)*(1-p_on)*(p_on-0.5+m_n));
520 *   double mp_t = psi_i*(tau*p_ot);
521 *   double kp_t = grad_psi_i*(tau*grad_p_ot);
522 *   double m_t = (alpha/M_PI)*std::atan(gamma*(t_eq-t_ot));
523 *   double nl_t = psi_i*(p_ot)*(1-p_ot)*(p_ot-0.5+m_t);
524 * @endcode
525 *
526 * Adding random noise at the interface
527 *
528 * @code
529 *   nl_n -= psi_i*(a*khi*(p_on)*(1-p_on));
530 *   nl_t -= psi_i*(a*khi*(p_ot)*(1-p_ot));
531 *  
532 *   double f1n = mp_n + time_step*theta*kp_n - time_step*theta*nl_n - mp_t + time_step*(1-theta)*kp_t - time_step*(1-theta)*nl_t; //f1 at last newton iteration
533 *  
534 *   double mt_n = phi_i*(t_on);
535 *   double kt_n = grad_phi_i*(grad_t_on);
536 *   double mpt_n = phi_i*(latent_heat*p_on);
537 *   double mt_t = phi_i*(t_ot);
538 *   double kt_t = grad_phi_i*(grad_t_ot);
539 *   double mpt_t = phi_i*(latent_heat*p_ot);
540 *  
541 *   double f2n = mt_n + time_step*theta*kt_n - mpt_n - mt_t + time_step*(1-theta)*kt_t + mpt_t; //f2 at last newton iteration
542 *  
543 * @endcode
544 *
545 * Assembling RHS vector
546 *
547 * @code
548 *   cell_rhs(i) -= (f1n + f2n)*fe_values.JxW(q);
549 *   }
550 *   }
551 *  
552 *   cell->get_dof_indices(local_dof_indices);
553 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
554 *   {
555 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
556 *   jacobian_matrix.add(local_dof_indices[i],
557 *   local_dof_indices[j],
558 *   cell_matrix(i, j));
559 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
560 *   }
561 *   }
562 *   }
563 *  
564 *   jacobian_matrix.compress(VectorOperation::add);
565 *   system_rhs.compress(VectorOperation::add);
566 *  
567 * @endcode
568 *
569 * Applying zero BC
570 *
571 * @code
572 *   std::map<types::global_dof_index, double> boundary_values;
574 *   1,
576 *   boundary_values);
577 *   MatrixTools::apply_boundary_values(boundary_values,
578 *   jacobian_matrix,
579 *   solution_update,
580 *   system_rhs, false);
581 *  
582 *   jacobian_matrix.compress(VectorOperation::insert);
583 *  
584 *   system_rhs.compress(VectorOperation::insert);
585 *   }
586 * @endcode
587
588
589<a name="ann-get_random_number.cpp"></a>
590<h1>Annotated version of get_random_number.cpp</h1>
591 *
592 *
593 *
594 *
595 * @code
596 *   /* -----------------------------------------------------------------------------
597 *   *
598 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
599 *   * Copyright (C) 2024 by Umair Hussain
600 *   *
601 *   * This file is part of the deal.II code gallery.
602 *   *
603 *   * -----------------------------------------------------------------------------
604 *   */
605 *  
606 *   #include "PhaseFieldSolver.h"
607 *   #include <random>
608 *  
609 *   float PhaseFieldSolver::get_random_number()
610 *   {
611 *   static std::default_random_engine e;
612 *   static std::uniform_real_distribution<> dis(-0.5, 0.5); // returns a random number in the range of -0.5 to 0.5
613 *   return dis(e);
614 *   }
615 * @endcode
616
617
618<a name="ann-grid_dof.cpp"></a>
619<h1>Annotated version of grid_dof.cpp</h1>
620 *
621 *
622 *
623 *
624 * @code
625 *   /* -----------------------------------------------------------------------------
626 *   *
627 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
628 *   * Copyright (C) 2024 by Umair Hussain
629 *   *
630 *   * This file is part of the deal.II code gallery.
631 *   *
632 *   * -----------------------------------------------------------------------------
633 *   */
634 *  
635 *   #include "PhaseFieldSolver.h"
636 *  
637 *   void PhaseFieldSolver::make_grid_and_dofs() {
638 * @endcode
639 *
640 * Reading mesh
641 *
642 * @code
643 *   gridin.attach_triangulation(triangulation);
644 *   std::ifstream f("mesh/Kobayashi_mesh100x400.msh");
645 *   gridin.read_msh(f);
646 *  
647 *  
649 *   dof_handler.distribute_dofs(fe);
650 *   DoFRenumbering::subdomain_wise(dof_handler);
651 *  
652 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
653 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
654 *  
655 *   const std::vector<IndexSet> locally_owned_dofs_per_proc =
657 *   const IndexSet locally_owned_dofs =
658 *   locally_owned_dofs_per_proc[this_mpi_process];
659 *   jacobian_matrix.reinit(locally_owned_dofs,
660 *   locally_owned_dofs,
661 *   dsp,
662 *   mpi_communicator);
663 *   old_solution.reinit(locally_owned_dofs, mpi_communicator);
664 *   system_rhs.reinit(locally_owned_dofs, mpi_communicator);
665 *   conv_solution.reinit(locally_owned_dofs, mpi_communicator);
666 *   solution_update.reinit(locally_owned_dofs, mpi_communicator);
667 *  
668 *   conv_solution_np.reinit(dof_handler.n_dofs());
669 *   old_solution_np.reinit(dof_handler.n_dofs());
670 *   }
671 * @endcode
672
673
674<a name="ann-main.cpp"></a>
675<h1>Annotated version of main.cpp</h1>
676 *
677 *
678 *
679 *
680 * @code
681 *   /* -----------------------------------------------------------------------------
682 *   *
683 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
684 *   * Copyright (C) 2024 by Umair Hussain
685 *   *
686 *   * This file is part of the deal.II code gallery.
687 *   *
688 *   * -----------------------------------------------------------------------------
689 *   */
690 *  
691 *   #include <iostream>
692 *   #include "PhaseFieldSolver.h"
693 *  
694 *   int main(int argc, char **argv) {
695 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
696 *   PhaseFieldSolver phasefieldsolver;
697 *   phasefieldsolver.run();
698 *   return 0;
699 *   }
700 * @endcode
701
702
703<a name="ann-output_results.cpp"></a>
704<h1>Annotated version of output_results.cpp</h1>
705 *
706 *
707 *
708 *
709 * @code
710 *   /* -----------------------------------------------------------------------------
711 *   *
712 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
713 *   * Copyright (C) 2024 by Umair Hussain
714 *   *
715 *   * This file is part of the deal.II code gallery.
716 *   *
717 *   * -----------------------------------------------------------------------------
718 *   */
719 *  
720 *   #include "PhaseFieldSolver.h"
721 *  
722 *   void PhaseFieldSolver::output_results(const unsigned int timestep_number) const {
723 *   const Vector<double> localized_solution(old_solution);
724 *  
725 * @endcode
726 *
727 * using only one process to output the result
728 *
729 * @code
730 *   if (this_mpi_process == 0)
731 *   {
732 *   DataOut<2> data_out;
733 *   data_out.attach_dof_handler(dof_handler);
734 *  
735 *   std::vector<std::string> solution_names;
736 *   solution_names.emplace_back ("p");
737 *   solution_names.emplace_back ("T");
738 *  
739 *   data_out.add_data_vector(localized_solution, solution_names);
740 *   const std::string filename =
741 *   "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtk";
742 *   DataOutBase::VtkFlags vtk_flags;
743 *   vtk_flags.compression_level =
744 *   DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed;
745 *   data_out.set_flags(vtk_flags);
746 *   std::ofstream output(filename);
747 *  
748 *   data_out.build_patches();
749 *   data_out.write_vtk(output);
750 *   }
751 *   }
752 * @endcode
753
754
755<a name="ann-run.cpp"></a>
756<h1>Annotated version of run.cpp</h1>
757 *
758 *
759 *
760 *
761 * @code
762 *   /* -----------------------------------------------------------------------------
763 *   *
764 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
765 *   * Copyright (C) 2024 by Umair Hussain
766 *   *
767 *   * This file is part of the deal.II code gallery.
768 *   *
769 *   * -----------------------------------------------------------------------------
770 *   */
771 *  
772 *   #include "PhaseFieldSolver.h"
773 *   #include <time.h>
774 *  
775 *   void PhaseFieldSolver::run() {
776 *   make_grid_and_dofs();
777 *   pcout << "Processors used: " << n_mpi_processes << std::endl;
778 *   pcout << " Number of degrees of freedom: " << dof_handler.n_dofs()
779 *   << " (by partition:";
780 *   for (unsigned int p = 0; p < n_mpi_processes; ++p)
781 *   pcout << (p == 0 ? ' ' : '+')
783 *   p));
784 *   pcout << ")" << std::endl;
785 * @endcode
786 *
787 * Initialise the solution
788 *
789 * @code
790 *   InitialValues initial_value;
791 *   VectorTools::interpolate(dof_handler,
792 *   initial_value,
793 *   old_solution);
794 *   VectorTools::interpolate(dof_handler,
795 *   initial_value,
796 *   conv_solution);
797 * @endcode
798 *
799 * Applying Boundary Conditions at t=0
800 *
801 * @code
802 *   applying_bc();
803 * @endcode
804 *
805 * Plotting initial solution
806 *
807 * @code
808 *   output_results(0);
809 *  
810 * @endcode
811 *
812 * Time steps begin here:
813 *
814 * @code
815 *   unsigned int timestep_number = 1;
816 *   for (; time <= final_time; time += time_step, ++timestep_number) {
817 *  
818 *   pcout << "Time step " << timestep_number << " at t=" << time+time_step
819 *   << std::endl;
820 *  
821 *   conv_solution.operator=(old_solution); // initialising the newton solution
822 *  
823 * @endcode
824 *
825 * Newton-Raphson iterations begin here:
826 *
827 * @code
828 *   for (unsigned int it = 1; it <= 100; ++it) {
829 *   pcout << "Newton iteration number:" << it << std::endl;
830 *  
831 *   if (it == 100) {
832 *   pcout << "Convergence Failure!!!!!!!!!!!!!!!" << std::endl;
833 *   std::exit(0);
834 *   }
835 * @endcode
836 *
837 * Saving parallel vectors as non-parallel ones
838 *
839 * @code
840 *   conv_solution_np = conv_solution;
841 *   old_solution_np = old_solution;
842 * @endcode
843 *
844 * Initialise the delta solution as zero
845 *
846 * @code
847 *   VectorTools::interpolate(dof_handler,
849 *   solution_update);
850 *   solution_update.compress(VectorOperation::insert);
851 * @endcode
852 *
853 * Assemble Jacobian and Residual
854 *
855 * @code
856 *   assemble_system();
857 * @endcode
858 *
859 * Solving to get delta solution
860 *
861 * @code
862 *   solve();
863 * @endcode
864 *
865 * Checking for convergence
866 *
867 * @code
868 *   double residual_norm = system_rhs.l2_norm(); //the norm of residual should converge to zero as the solution converges
869 * @endcode
870 *
871 * pcout << "Nothing wrong till here!!!!!!" << std::endl;
872 *
873 * @code
874 *   pcout << "the residual is:" << residual_norm << std::endl;
875 *   if (residual_norm <= (1e-4)) {
876 *   pcout << "Solution Converged!" << std::endl;
877 *   break; //Break to next time step if the N-R iterations converge
878 *   }
879 *   }
880 * @endcode
881 *
882 * Transfer the converged solution to the old_solution vector to plot output
883 *
884 * @code
885 *   old_solution.operator=(conv_solution);
886 *   old_solution.compress(VectorOperation::insert);
887 * @endcode
888 *
889 * output the solution at only specific number of time steps
890 *
891 * @code
892 *   if (timestep_number%10 == 0)
893 *   output_results(timestep_number);
894 *   }
895 *   }
896 * @endcode
897
898
899<a name="ann-solve.cpp"></a>
900<h1>Annotated version of solve.cpp</h1>
901 *
902 *
903 *
904 *
905 * @code
906 *   /* -----------------------------------------------------------------------------
907 *   *
908 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
909 *   * Copyright (C) 2024 by Umair Hussain
910 *   *
911 *   * This file is part of the deal.II code gallery.
912 *   *
913 *   * -----------------------------------------------------------------------------
914 *   */
915 *  
916 *   #include "PhaseFieldSolver.h"
917 *  
918 *   void PhaseFieldSolver::solve(){
919 * @endcode
920 *
921 * Using a direct parallel solver
922 *
923 * @code
924 *   SolverControl cn;
925 *   PETScWrappers::SparseDirectMUMPS A_direct(cn, mpi_communicator);
926 *   A_direct.solve(jacobian_matrix, solution_update, system_rhs);
927 * @endcode
928 *
929 * Updating the solution by adding the delta solution
930 *
931 * @code
932 *   conv_solution.add(1, solution_update);
933 *   conv_solution.compress(VectorOperation::add);
934 *   }
935 * @endcode
936
937
938*/
std::vector< bool > component_mask
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
Definition fe_q.h:554
Definition point.h:111
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
const Event initial
Definition event.cc:64
void subdomain_wise(DoFHandler< dim, spacedim > &dof_handler)
void random(DoFHandler< dim, spacedim > &dof_handler)
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
unsigned int count_dofs_with_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, const types::subdomain_id subdomain)
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
@ matrix
Contents is actually a matrix.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * begin(VectorType &V)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:92
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:107
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={})
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
long double gamma(const unsigned int n)
void copy(const T *begin, const T *end, U *dest)
int(& functions)(const void *v1, const void *v2)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
inline ::VectorizedArray< Number, width > atan(const ::VectorizedArray< Number, width > &x)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DataOutBase::CompressionLevel compression_level