i am using the implimenation by Jaekwang Kim to model polymer flows. The major change that I want to do is to output convergence tables like in Step-7 and Step-22, which means I have to have an ExactSolution.
The code is encountering an error because the dimension of the solutions do not match at a some point when generating output: *An error occurred in line <820> of file </home/jurombo/Downloads/dealii-9.3.2/include/deal.II/numerics/vector_tools_integrate_difference.templates.h> in function void dealii::VectorTools::internal::do_integrate_difference(const dealii::hp::MappingCollection<dim, spacedim>&, const dealii::DoFHandler<dim, spacedim>&, const InVector&, const dealii::Function<spacedim, typename InVector::value_type>&, OutVector&, const dealii::hp::QCollection<dim>&, const dealii::VectorTools::NormType&, const dealii::Function<spacedim>*, double) [with int dim = 2; int spacedim = 2; InVector = dealii::BlockVector<double>; OutVector = dealii::Vector<double>; typename InVector::value_type = double] The violated condition was: exact_solution.n_components == n_components Additional information: Dimension 3 not equal to 7. Stacktrace: ----------- #0 /usr/local/share/dealii/lib/libdeal_II.g.so.9.3.2: #1 /usr/local/share/dealii/lib/libdeal_II.g.so.9.3.2: void dealii::VectorTools::integrate_difference<2, dealii::BlockVector<double>, dealii::Vector<double>, 2>(dealii::Mapping<2, 2> const&, dealii::DoFHandler<2, 2> const&, dealii::BlockVector<double> const&, dealii::Function<2, dealii::BlockVector<double>::value_type> const&, dealii::Vector<double>&, dealii::Quadrature<2> const&, dealii::VectorTools::NormType const&, dealii::Function<2, double> const*, double) #2 ./oldroyd_bn: Viscoelastic::StokesProblem<2>::output_results(unsigned int) #3 ./oldroyd_bn: Viscoelastic::StokesProblem<2>::run() #4 ./oldroyd_bn: main --------------------------------------------------------* I have checked the ExactSolution in the code to try and figure wher the error is occurring without success. Can anyone point me the direction to look. -- The information in this message is confidential and legally privileged. It is intended solely for the addressee(s). Access to this message by anyone else is unauthorized. If received in error, please accept our apologies and notify the sender immediately. You must also delete the original message from your machine. If you are not the intended recipient, any use, disclosure, copying, distribution or action taken in reliance of it, is prohibited and may be unlawful. The information, attachments, opinions or advice contained in this email are not the views or opinions of Harare Institute of Technology, its subsidiaries or affiliates. Although this email and any attachments are believed to be free of any virus or other defects which might affect any computer or IT system into which they are received, no responsibility is accepted by Harare Institute of Technology and/or its subsidiaries for any loss or damage arising in any way from the receipt or use thereof. -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/a7464d4b-87df-4e49-80cf-7a0900fe1a70n%40googlegroups.com.
// Developed by Jaekwang Kim // Welcome to the world of viscoelastic fluids // Coding from 18th May ~ // Viscoelastic Models // Oldroyd B model // Upper Convective Maxwel polymer(UCM) + a Newtonian solvent // Mixed Finite Element (also called viscous formulation) will be considered // (a) Starting from large viscous contribution from a Newtonian Solvent // (b) Low Wissenberg Number // (c) Stabilization Method // (d) Be careful on the LBB condition; Continuous Discretization of T_tensor? //what is the convention of grad_u_[0][1] ? #include <deal.II/base/quadrature_lib.h> #include <deal.II/base/logstream.h> #include <deal.II/base/function.h> #include <deal.II/base/utilities.h> #include <deal.II/base/convergence_table.h> #include <deal.II/lac/block_vector.h> #include <deal.II/lac/full_matrix.h> #include <deal.II/lac/block_sparse_matrix.h> #include <deal.II/lac/solver_cg.h> // CG used for viscous flow #include <deal.II/lac/solver_gmres.h> //Matrix solver for transport equation, unsymmetric matrix #include <deal.II/lac/precondition.h> #include <deal.II/lac/affine_constraints.h> #include <deal.II/lac/dynamic_sparsity_pattern.h> #include <deal.II/lac/sparse_direct.h> #include <deal.II/lac/sparse_ilu.h> #include <deal.II/grid/tria.h> #include <deal.II/grid/grid_generator.h> #include <deal.II/grid/tria_accessor.h> #include <deal.II/grid/tria_iterator.h> #include <deal.II/grid/grid_tools.h> #include <deal.II/grid/grid_refinement.h> #include <deal.II/grid/manifold_lib.h> #include <deal.II/grid/grid_in.h> #include <deal.II/grid/grid_out.h> #include <deal.II/dofs/dof_handler.h> #include <deal.II/dofs/dof_renumbering.h> #include <deal.II/dofs/dof_accessor.h> #include <deal.II/dofs/dof_tools.h> #include <deal.II/fe/fe_q.h> #include <deal.II/fe/fe_system.h> #include <deal.II/fe/fe_values.h> #include <deal.II/fe/mapping_q.h> #include <deal.II/numerics/vector_tools.h> #include <deal.II/numerics/matrix_tools.h> #include <deal.II/numerics/data_out.h> #include <deal.II/base/parameter_handler.h> #include <deal.II/numerics/error_estimator.h> #include <deal.II/numerics/solution_transfer.h> #include <iostream> #include <fstream> #include <sstream> #include "headers/post_processor.h" #include "headers/precondition.h" #include "headers/bcfunctions.h" #include "headers/fluidclass.h" #include "headers/flowclass.h" namespace Viscoelastic { using namespace dealii; template <int dim> class StokesProblem { public: StokesProblem (const unsigned int degree); void run (); private: void read_mesh (); void setup_dofs (); void setup_bc_constraints (); void assemble_stokes (); // Flow system void solve_stokes (); // Solve with CG void assemble_polymer_system (); // UCM advection void solve_polymer_system (); // Solve with GMRES void refine_mesh (const unsigned int refinement_cycle); void output_results (const unsigned int refinement_cycle); const unsigned int degree; Triangulation<dim> triangulation; FESystem<dim> fe; DoFHandler<dim> dof_handler; MappingQ<dim> mapping; AffineConstraints<double> constraints; BlockSparsityPattern sparsity_pattern; BlockSparseMatrix<double> system_matrix; BlockVector<double> solution; // BlockVector<double> exact_solution; BlockVector<double> interpolated_exact_solution; BlockVector<double> system_rhs; // Previous solution for Piccard iterations BlockVector<double> previous_solution; AdvectionBoundaryValues<dim> advection_boundary_values; ConstantU<dim> Uplate; // BC Function Oldroyd_FluidClass fluid; double lam1 = fluid.lam1; double lam2 = fluid.lam2; double etap = fluid.etap; double etas = fluid.etas; Oldroyd_FlowClass flow; double Ux = flow.Ux; double h = flow.h; double alpha = flow.alpha; double Wi = flow.Wi; std::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner; ConvergenceTable convergence_table; }; // @sect3{Boundary values and right hand side} template <int dim> class BoundaryValues : public Function<dim> { public: BoundaryValues () : Function<dim>(dim+1) {} virtual void vector_value (const Point<dim> &p, Vector<double> &value) const; const double U = 1.6; // const double Ux = flow.Ux; }; template <int dim> void BoundaryValues<dim>::vector_value (const Point<dim> & p, Vector<double> &values) const { const double rad = 1.0 - p.square(); const double Ux = .0125; values(0) = 0.0; values(1) = 0.0; values(2) = 0.0; values(3) = Ux*(1 - rad); values(4) = 0.0; values(5) = 0.0; values(6) = 0.0; } // the right hand side template <int dim> class RightHandSide : public Function<dim> { public: RightHandSide () : Function<dim>(dim+1) {} virtual void vector_value (const Point<dim> &p, Vector<double> &value) const; }; template <int dim> void RightHandSide<dim>::vector_value (const Point<dim> &p, Vector<double> &values) const { values(0) = -8.0*p[1] - 2.0*p[0]; values(1) = 8.0*p[0] - 2.0*p[1]; values(2) = 0.0; values(3) = 0.0; values(4) = 0.0; values(5) = 0.0; values(6) = 0.0; } // ExactSolution template <int dim> class ExactSolution : public Function<dim> { public: ExactSolution () : Function<dim>(dim+1) {} virtual void vector_value (const Point<dim> &p, Vector<double> &value) const; Tensor<1,dim> gradient (const Point<dim> &p, const unsigned int component) const; const double R = 1.0; const double alpha = alpha = etas/(etap+etas); const double Wi = lam1*(Ux/h); //const double U = 1.6; //const double Ux = .125; }; template <int dim> void ExactSolution<dim>::vector_value (const Point<dim> &p, Vector<double> &values) const { Assert (values.size() == dim+1, ExcDimensionMismatch (values.size(), dim+1)); const double rad = 1.0 - p.square(); values(0) = -p[1]*rad; values(1) = p[0]*rad; values(2) = rad; values(3) = -2*p[1]*(1 - alpha)*U/R*R; values(4) = 8*Wi*(1 - alpha)*p.square()*Ux*Ux/R*R; values(5) = 8*Wi*(1 - alpha)*p.square()*Ux*Ux/R*R; values(6) = 0.0; } template <int dim> Tensor<1,dim> ExactSolution<dim>::gradient (const Point<dim> &p, const unsigned int component) const { Assert (component < this->n_components, ExcIndexRange (component, 0, this->n_components)); double x = p[0]; double y = p[1]; Tensor<1,dim> gradient; switch (component) { //velocity case 0: gradient[0]=2*x*y; gradient[1]=-1+x*x+3*y*y; break; case 1: gradient[0]=1-3*x*x-y*y; gradient[1]=-2*x*y; break; //pressure case 2: gradient[0]=-2*x; gradient[1]=-2*y; break; // stress case 3: gradient[0]=-2*(1-alpha)*U/R*R; gradient[1]=16*Wi*(1-alpha)*p[1]*Ux*Ux/R*R; gradient[2]=16*Wi*(1-alpha)*p[1]*Ux*Ux/R*R; gradient[3]=0.0; break; } return gradient; } template <int dim> StokesProblem<dim>::StokesProblem (const unsigned int degree) : degree (degree), triangulation (Triangulation<dim>::allow_anisotropic_smoothing), fe (FE_Q<dim>(degree+1), dim, // u FE_Q<dim>(degree) , 1 , // p FE_Q<dim>(degree+1) , dim*dim ), //tau_p dof_handler (triangulation), mapping(2) {} template <int dim> void StokesProblem<dim>::setup_bc_constraints () { constraints.clear (); FEValuesExtractors::Vector velocities(0); FEValuesExtractors::Scalar xvel(0); FEValuesExtractors::Scalar yvel(1); DoFTools::make_hanging_node_constraints (dof_handler, constraints); // Inlet VectorTools::interpolate_boundary_values (dof_handler, 4, Uinlet<dim>(), constraints, fe.component_mask(velocities)); // Ends - parallel flow VectorTools::interpolate_boundary_values (dof_handler, 2, ZeroFunction<dim>(dim+1+dim*dim), constraints, fe.component_mask(yvel)); // Side - Uniform wall motion: VectorTools::interpolate_boundary_values (dof_handler, 5, Uplate, constraints, fe.component_mask(velocities)); // Fixed Wall --- 1 is flat sections, > 5 are rounded sections with manifolds VectorTools::interpolate_boundary_values (dof_handler, 1, ZeroFunction<dim>(dim+1+dim*dim), constraints, fe.component_mask(velocities)); VectorTools::interpolate_boundary_values (dof_handler, 6, ZeroFunction<dim>(dim+1+dim*dim), constraints, fe.component_mask(velocities)); VectorTools::interpolate_boundary_values (dof_handler, 7, ZeroFunction<dim>(dim+1+dim*dim), constraints, fe.component_mask(velocities)); VectorTools::interpolate_boundary_values (dof_handler, 8, ZeroFunction<dim>(dim+1+dim*dim), constraints, fe.component_mask(velocities)); constraints.close (); } template <int dim> void StokesProblem<dim>::setup_dofs () { std::cout << " Set_Up_DOF -- Start" <<std::endl; A_preconditioner.reset (); system_matrix.clear (); dof_handler.distribute_dofs (fe); DoFRenumbering::Cuthill_McKee (dof_handler); std::vector<unsigned int> block_component (dim+1+dim*dim,0); // initial 0 for u block_component[dim] = 1; // pressure for (unsigned int i=0;i<dim*dim;i++) block_component[dim+1+i] = 2; //block for polymer_T variable... DoFRenumbering::component_wise (dof_handler, block_component); std::vector<types::global_dof_index> dofs_per_block (3); DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component); const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1], n_s = dofs_per_block[2]; std::cout << " n_u: " << n_u << " n_p: " << n_p << " n_s: " << n_s << std::endl; // Apply BC Constraints setup_bc_constraints(); { BlockDynamicSparsityPattern dsp (3,3); dsp.block(0,0).reinit (n_u, n_u); dsp.block(1,0).reinit (n_p, n_u); dsp.block(0,1).reinit (n_u, n_p); dsp.block(1,1).reinit (n_p, n_p); dsp.block(1,2).reinit (n_p, n_s); dsp.block(2,1).reinit (n_s, n_p); dsp.block(2,0).reinit (n_s, n_u); dsp.block(0,2).reinit (n_u, n_s); dsp.block(2,2).reinit (n_s, n_s); dsp.collect_sizes(); DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, true); sparsity_pattern.copy_from (dsp); } system_matrix.reinit (sparsity_pattern); solution.reinit (3); solution.block(0).reinit (n_u); solution.block(1).reinit (n_p); solution.block(2).reinit (n_s); solution.collect_sizes (); interpolated_exact_solution.reinit (3); interpolated_exact_solution.block(0).reinit (n_u); interpolated_exact_solution.block(1).reinit (n_p); interpolated_exact_solution.block(2).reinit (n_s); interpolated_exact_solution.collect_sizes (); system_rhs.reinit (3); system_rhs.block(0).reinit (n_u); system_rhs.block(1).reinit (n_p); system_rhs.block(2).reinit (n_s); system_rhs.collect_sizes (); std::cout << " Active cells: "<< triangulation.n_active_cells() << std::endl << " DoFs: " << dof_handler.n_dofs() << " (" << n_u << '+' << n_p << '+'<< n_s <<')' << std::endl; } // Assemble Momentum and Continuity Equation template <int dim> void StokesProblem<dim>::assemble_stokes () { // Momentum and Continuity Equation system_matrix=0; system_rhs=0; const MappingQ<dim> mapping (degree); QGauss<dim> quadrature_formula(2*degree+1); FEValues<dim> fe_values (mapping, fe, quadrature_formula, update_values | update_quadrature_points | update_JxW_values | update_gradients); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell); Vector<double> local_rhs (dofs_per_cell); std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell); std::vector<Vector<double> > rhs_values (n_q_points, Vector<double>(dim+1+dim*dim)); const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); const FEValuesExtractors::Scalar Txx (dim+1); const FEValuesExtractors::Scalar Txy (dim+2); const FEValuesExtractors::Scalar Tyx (dim+3); const FEValuesExtractors::Scalar Tyy (dim+4); std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell); std::vector<Tensor<2,dim> > grad_phi_u (dofs_per_cell); std::vector<double> div_phi_u (dofs_per_cell); std::vector<double> phi_p (dofs_per_cell); std::vector<double> phi_ur (dofs_per_cell); std::vector<double> phi_i_s_xx(dofs_per_cell); std::vector<double> phi_i_s_xy(dofs_per_cell); std::vector<double> phi_i_s_yx(dofs_per_cell); std::vector<double> phi_i_s_yy(dofs_per_cell); // Variables related to previous solution std::vector<SymmetricTensor<2,dim> > local_symgrad_u (n_q_points); std::vector<SymmetricTensor<2,dim> > local_T (n_q_points); std::vector<double> local_Txx (n_q_points); std::vector<double> local_Txy (n_q_points); std::vector<double> local_Tyx (n_q_points); std::vector<double> local_Tyy (n_q_points); typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) { fe_values.reinit (cell); local_matrix = 0; local_rhs = 0; fe_values[velocities].get_function_symmetric_gradients (previous_solution, local_symgrad_u); fe_values[Txx].get_function_values(previous_solution,local_Txx); fe_values[Txy].get_function_values(previous_solution,local_Txy); fe_values[Tyx].get_function_values(previous_solution,local_Tyx); fe_values[Tyy].get_function_values(previous_solution,local_Tyy); for (unsigned int q=0; q<n_q_points; ++q) { for (unsigned int k=0; k<dofs_per_cell; ++k) { symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k, q); grad_phi_u[k] = fe_values[velocities].gradient (k, q); div_phi_u[k] = fe_values[velocities].divergence (k, q); phi_p[k] = fe_values[pressure].value (k, q); } for (unsigned int i=0; i<dofs_per_cell; ++i) { for (unsigned int j=0; j<=i; ++j) { local_matrix(i,j) += ( //Stokes Formulation //(a) Solvent Contribution 2.*etas*(symgrad_phi_u[i]*symgrad_phi_u[j]) //(b) presure term and continuity term - div_phi_u[i]*phi_p[j] - phi_p[i]*div_phi_u[j] + phi_p[i]*phi_p[j] ) * fe_values.JxW(q); } local_rhs(i) += ( grad_phi_u[i][0][0]*(local_Txx[q]) +grad_phi_u[i][1][0]*(local_Txy[q]) +grad_phi_u[i][0][1]*(local_Tyx[q]) +grad_phi_u[i][1][1]*(local_Tyy[q]) )*fe_values.JxW(q); } } for (unsigned int i=0; i<dofs_per_cell; ++i) for (unsigned int j=i+1; j<dofs_per_cell; ++j) local_matrix(i,j) = local_matrix(j,i); cell->get_dof_indices (local_dof_indices); constraints.distribute_local_to_global (local_matrix, local_rhs, local_dof_indices, system_matrix, system_rhs); } // apply Dirichlet boundary condition on velocity std::map<types::global_dof_index,double> boundary_values; A_preconditioner = std::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type()); A_preconditioner->initialize (system_matrix.block(0,0), typename InnerPreconditioner<dim>::type::AdditionalData()); } // Assemble Constitutive Equation template <int dim> void StokesProblem<dim>::assemble_polymer_system () { std::cout << " Assemble Polymer Constitutive " << std::endl; const MappingQ<dim> mapping (degree); QGauss<dim> quadrature_formula(2*degree+1); QGauss<dim-1> face_quadrature_formula(2*degree+1); FEValues<dim> fe_values (mapping, fe, quadrature_formula, update_values | update_quadrature_points | update_JxW_values | update_gradients); FEFaceValues<dim> fe_face_values (mapping, fe, face_quadrature_formula, update_values | update_normal_vectors | update_gradients | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); const unsigned int n_face_q_points = face_quadrature_formula.size(); std::vector<Vector<double> > solution_values_face(n_face_q_points, Vector<double>(dim+1+dim*dim)); std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell); FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell); Vector<double> local_rhs (dofs_per_cell); // Flow solution; std::vector<SymmetricTensor<2,dim> > local_symgrad_phi_u (n_q_points); std::vector<Tensor<2,dim> > local_grad_phi_u (n_q_points); std::vector<Tensor<1,dim> > local_phi_u (n_q_points); std::vector<Tensor<1,1> > local_pres (n_q_points); const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar Txx (dim+1); const FEValuesExtractors::Scalar Txy (dim+2); const FEValuesExtractors::Scalar Tyx (dim+3); const FEValuesExtractors::Scalar Tyy (dim+4); const FEValuesExtractors::Scalar pres (dim); std::vector<Tensor<1,dim>> face_advection_directions (n_face_q_points); std::vector<Vector<double>> face_boundary_values (n_face_q_points,Vector<double>(dim+1+dim*dim)); Tensor<1,dim> normal; double vel_magnitude; double phi_i_s_xx, phi_i_s_xy, phi_i_s_yy, phi_i_s_yx; double phi_j_s_xx, phi_j_s_xy, phi_j_s_yy, phi_j_s_yx; Tensor<1,dim> grad_phi_i_s_xx, grad_phi_i_s_xy, grad_phi_i_s_yy, grad_phi_i_s_yx; Tensor<1,dim> grad_phi_j_s_xx, grad_phi_j_s_xy, grad_phi_j_s_yy, grad_phi_j_s_yx; typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) { fe_values.reinit (cell); local_matrix = 0; local_rhs = 0; // Use of local solution fe_values[velocities].get_function_symmetric_gradients (solution, local_symgrad_phi_u); fe_values[velocities].get_function_gradients (solution, local_grad_phi_u); fe_values[velocities].get_function_values (solution, local_phi_u); // fe_values[pres].get_function_values (solution, local_pres); // Stabilization double delta = 0.1 * cell->diameter (); for (unsigned int q=0; q<n_q_points; ++q) { vel_magnitude = std::sqrt(local_phi_u[q]*local_phi_u[q]); for (unsigned int i=0; i<dofs_per_cell; ++i) { phi_i_s_xx = fe_values[Txx].value(i,q); phi_i_s_xy = fe_values[Txy].value(i,q); phi_i_s_yx = fe_values[Tyx].value(i,q); phi_i_s_yy = fe_values[Tyy].value(i,q); grad_phi_i_s_xx = fe_values[Txx].gradient (i, q); grad_phi_i_s_xy = fe_values[Txy].gradient (i, q); grad_phi_i_s_yx = fe_values[Tyx].gradient (i, q); grad_phi_i_s_yy = fe_values[Tyy].gradient (i, q); for (unsigned int j=0; j<dofs_per_cell; ++j) { phi_j_s_xx = fe_values[Txx].value(j,q); phi_j_s_xy = fe_values[Txy].value(j,q); phi_j_s_yx = fe_values[Tyx].value(j,q); phi_j_s_yy = fe_values[Tyy].value(j,q); grad_phi_j_s_xx = fe_values[Txx].gradient (j, q); grad_phi_j_s_xy = fe_values[Txy].gradient (j, q); grad_phi_j_s_yx = fe_values[Tyx].gradient (j, q); grad_phi_j_s_yy = fe_values[Tyy].gradient (j, q); local_matrix(i,j) += ( //self term + phi_i_s_xx * phi_j_s_xx + phi_i_s_xy * phi_j_s_xy + phi_i_s_yx * phi_j_s_yx + phi_i_s_yy * phi_j_s_yy + lam1*( //advection + phi_i_s_xx * (local_phi_u[q]*grad_phi_j_s_xx ) + phi_i_s_xy * (local_phi_u[q]*grad_phi_j_s_xy ) + phi_i_s_yx * (local_phi_u[q]*grad_phi_j_s_yx ) + phi_i_s_yy * (local_phi_u[q]*grad_phi_j_s_yy ) // rotation term // // Gradient Tensor convention (dux/dy)=local_grad_phi_u[q][0][1]) - phi_i_s_xx * ( 2.*local_grad_phi_u[q][0][0]*phi_j_s_xx + local_grad_phi_u[q][0][1]*phi_j_s_yx + local_grad_phi_u[q][0][1]*phi_j_s_xy ) - phi_i_s_xy * ( local_grad_phi_u[q][1][0]*phi_j_s_xx + local_grad_phi_u[q][0][0]*phi_j_s_xy + //(a) local_grad_phi_u[q][1][1]*phi_j_s_xy + //(b) local_grad_phi_u[q][0][1]*phi_j_s_yy ) - phi_i_s_yx * (local_grad_phi_u[q][1][0]*phi_j_s_xx + local_grad_phi_u[q][0][0]*phi_j_s_yx + //(c) local_grad_phi_u[q][1][1]*phi_j_s_yx + //(d) local_grad_phi_u[q][0][1]*phi_j_s_yy ) //actually (a+b+c+d) will be zero effectively - phi_i_s_yy * ( 2.*local_grad_phi_u[q][1][1]*phi_j_s_yy + local_grad_phi_u[q][1][0]*phi_j_s_yx + local_grad_phi_u[q][1][0]*phi_j_s_xy ) //std::cout << "local_" // stabilization //need to check + delta/vel_magnitude* ( local_phi_u[q]*grad_phi_i_s_xx //(u \cdot \nabla w) *local_phi_u[q]*grad_phi_j_s_xx + local_phi_u[q]*grad_phi_i_s_xy *local_phi_u[q]*grad_phi_j_s_xy + local_phi_u[q]*grad_phi_i_s_yy *local_phi_u[q]*grad_phi_j_s_yy + local_phi_u[q]*grad_phi_i_s_yx *local_phi_u[q]*grad_phi_j_s_yx ) ) )*fe_values.JxW(q); } //close 'j' cycle local_rhs(i) += etap*( + phi_i_s_xx * 2.*local_grad_phi_u[q][0][0] + phi_i_s_xy * (local_grad_phi_u[q][1][0] + local_grad_phi_u[q][0][1]) + phi_i_s_yy * 2.*local_grad_phi_u[q][1][1] + phi_i_s_yx * (local_grad_phi_u[q][1][0] + local_grad_phi_u[q][0][1]) )*fe_values.JxW(q); } // i } // q for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no) { if (cell->face(face_no)->boundary_id()==1) // strong inflow boundary condition { // inflow face fe_face_values.reinit (cell, face_no); fe_face_values.get_function_values (solution,solution_values_face); advection_boundary_values.vector_value_list (fe_face_values.get_quadrature_points(),face_boundary_values); for (unsigned int q=0; q<n_face_q_points; ++q) { Tensor<1,dim> present_u_face; for (unsigned int d=0; d<dim; ++d) present_u_face[d] = solution_values_face[q](d); //if-inflow condition if (fe_face_values.normal_vector(q) * present_u_face < 0) { for (unsigned int i=0; i<dofs_per_cell; ++i) { phi_i_s_xx = fe_face_values[Txx].value(i,q); phi_i_s_xy = fe_face_values[Txy].value(i,q); phi_i_s_yy = fe_face_values[Tyy].value(i,q); phi_i_s_yx = fe_face_values[Tyx].value(i,q); for (unsigned int j=0; j<dofs_per_cell; ++j) { phi_j_s_xx = fe_face_values[Txx].value(j,q); phi_j_s_xy = fe_face_values[Txy].value(j,q); phi_j_s_yy = fe_face_values[Tyy].value(j,q); phi_j_s_yx = fe_face_values[Tyx].value(j,q); local_matrix(i,j) -= (present_u_face *fe_face_values.normal_vector(q) *(+phi_i_s_xx*phi_j_s_xx +phi_i_s_xy*phi_j_s_xy +phi_i_s_yy*phi_j_s_yy +phi_i_s_yx*phi_j_s_yx ) )*fe_face_values.JxW(q); } //close cycle j local_rhs(i) -= present_u_face * fe_face_values.normal_vector(q)* ( + phi_i_s_xx * face_boundary_values[q][3] + phi_i_s_xy * face_boundary_values[q][4] + phi_i_s_yx * face_boundary_values[q][5] + phi_i_s_yy * face_boundary_values[q][6] )*fe_face_values.JxW(q); } // i } // if-inflow } // q } // cell } // face cell->get_dof_indices (local_dof_indices); constraints.distribute_local_to_global (local_matrix, local_rhs, local_dof_indices, system_matrix, system_rhs); } } template <int dim> void StokesProblem<dim>::solve_stokes () { const InverseMatrix<SparseMatrix<double>, typename InnerPreconditioner<dim>::type> A_inverse (system_matrix.block(0,0), *A_preconditioner); Vector<double> tmp (solution.block(0).size()); { Vector<double> schur_rhs (solution.block(1).size()); A_inverse.vmult (tmp, system_rhs.block(0)); system_matrix.block(1,0).vmult (schur_rhs, tmp); schur_rhs -= system_rhs.block(1); SchurComplement<typename InnerPreconditioner<dim>::type> schur_complement (system_matrix, A_inverse); SolverControl solver_control ( 10 * solution.block(1).size(), 1e-5*schur_rhs.l2_norm()); // solver_control here SolverCG<> cg (solver_control); SparseILU<double> preconditioner; preconditioner.initialize (system_matrix.block(1,1), SparseILU<double>::AdditionalData()); InverseMatrix<SparseMatrix<double>,SparseILU<double> > m_inverse (system_matrix.block(1,1), preconditioner); cg.solve (schur_complement, solution.block(1), schur_rhs, m_inverse); constraints.distribute (solution); std::cout << "Flow Solved:" << solver_control.last_step() << " outer CG Schur complement iterations for pressure" << std::endl; } { system_matrix.block(0,1).vmult (tmp, solution.block(1)); tmp *= -1; tmp += system_rhs.block(0); A_inverse.vmult (solution.block(0), tmp); constraints.distribute (solution); } } template <int dim> void StokesProblem<dim>::solve_polymer_system () { std::cout << " Solve Advection" << std::endl; SolverControl solver_control (std::pow(10,8), std::pow(10,-4)); unsigned int restart = 500; SolverGMRES< Vector<double> >::AdditionalData gmres_additional_data(restart+2); SolverGMRES< Vector<double> > solver(solver_control, gmres_additional_data); //'gmres_additional_data' means how much temporary vector you will going to use for grmres solver //the more the faster, but the more memory will be consummed. //make preconditioner SparseILU<double>::AdditionalData additional_data(0,1000); // (0 , additional diagonal terms) SparseILU<double> preconditioner; preconditioner.initialize (system_matrix.block(2,2), additional_data); solver.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2), preconditioner); constraints.distribute (solution); } template <int dim> void StokesProblem<dim>::output_results (const unsigned int refinement_cycle) { std::vector<std::string> solution_names(dim,"Velocity"); solution_names.push_back ("Pressure"); solution_names.push_back ("Txx"); solution_names.push_back ("Txy"); solution_names.push_back ("Tyx"); solution_names.push_back ("Tyy"); std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation (dim, DataComponentInterpretation::component_is_part_of_vector); data_component_interpretation .push_back (DataComponentInterpretation::component_is_scalar); for (unsigned int i=0; i<dim*dim; i++) data_component_interpretation .push_back (DataComponentInterpretation::component_is_scalar); DataOut<dim> data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, solution_names, DataOut<dim>::type_dof_data, data_component_interpretation); data_out.build_patches (); std::ostringstream filenameeps; filenameeps << "solution-"<< Utilities::int_to_string (refinement_cycle, 2)<< ".vtk"; std::ofstream output (filenameeps.str().c_str()); data_out.write_vtk (output); //data_out.write_tecplot (output); // Compute error const ComponentSelectFunction<dim> pressure_mask (dim, dim+1); const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim), dim+1); ExactSolution<dim> exact_solution; Vector<double> cellwise_errors (triangulation.n_active_cells()); // QTrapez<1> q_trapez; // QIterated<dim> quadrature (q_trapez, degree+2); QGauss<dim> quadrature (degree+3); const QTrapezoid<1> q_trapez; const QIterated<dim> q_iterated (q_trapez, 5); // L2 error (velocity) VectorTools::integrate_difference (mapping, dof_handler, solution, exact_solution, cellwise_errors, quadrature, VectorTools::L2_norm, &velocity_mask); const double u_L2_error = cellwise_errors.l2_norm(); // H1 semi norm (velocity) VectorTools::integrate_difference (mapping, dof_handler, solution, exact_solution, cellwise_errors, quadrature, VectorTools::H1_seminorm, &velocity_mask); const double u_H1_error = cellwise_errors.l2_norm(); // sup norm (velocity) VectorTools::integrate_difference (mapping, dof_handler, solution, exact_solution, cellwise_errors, q_iterated, VectorTools::Linfty_norm, &velocity_mask); const double u_Linfty_error = cellwise_errors.linfty_norm(); // L2 error (pressure) VectorTools::integrate_difference (mapping, dof_handler, solution, exact_solution, cellwise_errors, quadrature, VectorTools::L2_norm, &pressure_mask); const double p_L2_error = cellwise_errors.l2_norm(); // SUP error (pressure) VectorTools::integrate_difference (mapping, dof_handler, solution, exact_solution, cellwise_errors, q_iterated, VectorTools::Linfty_norm, &pressure_mask); const double p_Linfty_error = cellwise_errors.linfty_norm(); const unsigned int n_active_cells=triangulation.n_active_cells(); const unsigned int n_dofs=dof_handler.n_dofs(); convergence_table.add_value("cycle", refinement_cycle); convergence_table.add_value("cells", n_active_cells); convergence_table.add_value("dofs", dof_handler.n_dofs()); convergence_table.add_value("u_L2", u_L2_error); convergence_table.add_value("u_H1", u_H1_error); convergence_table.add_value("u_Linfty", u_Linfty_error); convergence_table.add_value("p_L2", p_L2_error); convergence_table.add_value("p_Linfty", p_Linfty_error); // output numerical solution DataOut<2> sol_out; sol_out.attach_dof_handler (dof_handler); sol_out.add_data_vector (solution, "numerical solution"); sol_out.build_patches (); std::ofstream sol_output ("num_sol.gpl"); sol_out.write_gnuplot (sol_output); // output nodal error if ( refinement_cycle == 0 ) { std::cout << "Output nodal error" << std::endl; VectorTools::interpolate (dof_handler, ExactSolution<dim>(), interpolated_exact_solution); interpolated_exact_solution -= solution; DataOut<2> err_out; err_out.attach_dof_handler (dof_handler); err_out.add_data_vector (interpolated_exact_solution, "pointwise_error"); err_out.build_patches (); std::ofstream err_output ("pointwise_err.gpl"); err_out.write_gnuplot (err_output); } Vector<double> interpolated_exact_solution (dof_handler.n_dofs()); VectorTools::interpolate (dof_handler, ExactSolution<dim>(), interpolated_exact_solution); for ( unsigned int i=0; i< solution.block(0).size(); i++ ) interpolated_exact_solution(i) -= solution.block(0)(i); DataOut<2> err_out; err_out.attach_dof_handler (dof_handler); err_out.add_data_vector (interpolated_exact_solution, "pointwise_error"); err_out.build_patches (); std::ofstream err_output ("pointwise_err.gpl"); err_out.write_gnuplot (err_output); // output result std::cout << " Number of active cells: " << n_active_cells << std::endl << " Number of degrees of freedom: " << n_dofs << std::endl; std::cout << "p L2 = " << p_L2_error << std::endl; std::cout << "p SUP= " << p_Linfty_error << std::endl; std::cout << "u L2 = " << u_L2_error << std::endl; std::cout << "u SUP= " << u_Linfty_error << std::endl; std::cout << "u H1 = " << u_H1_error << std::endl; } template <int dim> void StokesProblem<dim>::read_mesh () { { //Generate mesh and designate boundary /* Parallel plate const Point<2> end (5,1); const Point<2> start (0.0,0.0); GridGenerator::hyper_rectangle (triangulation, start, end, false); triangulation.refine_global (3); */ GridIn<dim> grid_in; grid_in.attach_triangulation (triangulation); //std::ifstream input_file("mesh/slant-rounded-mesh.inp"); std::ifstream input_file("mesh/new_mesh.inp"); std::cout <<"here" << std::endl; grid_in.read_ucd (input_file); std::cout <<"here" << std::endl; //Point<2> center1 (-0.2,-1.7); Point<2> center1 (-1.7,-0.2); static const SphericalManifold<dim> manifold_description1 (center1); triangulation.set_manifold (6, manifold_description1); //Point<2> center2 (-1,-1.3); Point<2> center2 (-1.3,-1.0); static const SphericalManifold<dim> manifold_description2 (center2); triangulation.set_manifold (7, manifold_description2); //Point<2> center3 (-0.6,1.5); Point<2> center3 (1.5,-0.6); static const SphericalManifold<dim> manifold_description3 (center3); triangulation.set_manifold (8, manifold_description3); triangulation.refine_global (3); } } template <int dim> void StokesProblem<dim>::run () { read_mesh(); //To read mesh file from outside for (unsigned int refinement_cycle = 0; refinement_cycle<3; ++refinement_cycle) { std::cout << "Refinement cycle " << refinement_cycle << std::endl; setup_dofs (); if (refinement_cycle < 1 ) { previous_solution = solution; previous_solution = 0.; } assemble_stokes (); solve_stokes (); assemble_polymer_system (); solve_polymer_system (); BlockVector<double> difference; int iteration_number=0 ; previous_solution =solution ; do{ iteration_number +=1; assemble_stokes (); solve_stokes (); assemble_polymer_system (); solve_polymer_system (); difference = solution; difference -= previous_solution; previous_solution=solution; std::cout << " Iteration Number: " << iteration_number << " Difference Norm: " << difference.l2_norm() << std::endl << std::flush; } while (difference.l2_norm()>pow(10,-9)* dof_handler.n_dofs()); output_results (refinement_cycle); // refine_mesh (refinement_cycle) ; } } } int main () { try { using namespace dealii; using namespace Viscoelastic; StokesProblem<2> flow_problem1(1); flow_problem1.run (); } catch (std::exception &exc) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Exception on processing: " << std::endl << exc.what() << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } return 0; }