Re: [deal.II] Assemble Righthand Side for vector-valued problem

[Wolfgang Bangerth]

On 02/09/2017 04:45 AM, Jaekwang Kim wrote:


The manufactured solution of velocities has not satisfies the continuity
equation.
(i.e. manufactured solution does not satisfy, div.u=0)...


Yes, that would do it :-)
Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Assemble Righthand Side for vector-valued problem

[Jaekwang Kim]

Dr. Bangerth 


Thank you, 
I just fixed up what was wrong... 
in most cases, it is usually easy problems. 

The manufactured solution of velocities has not satisfies the continuity 
equation. 
(i.e. manufactured solution does not satisfy, div.u=0)...

Jaekwang Kim  

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Assemble Righthand Side for vector-valued problem

[Jaekwang Kim]

Thank you for your advice. 

I would make the problem simpler. Do you get the right solution if you 
> have a constant viscosity? If you iterate, do you get the solution after 
> one iteration? 
>

Yes, I checked this. 
After one iteration, I get the solution after one iteration when constant 
viscosity case.
 

>
> I will also note that your manufactured solution is symmetric, but your 
> computational solution is not. This provides you with a powerful way 
> because you don't actually need to look at the convergence, it's enough 
> to check that on a *coarse* mesh the solution is *symmetric*. It may be 
> that solution is already non-symmetric on 1 cell, or on a 2x2 mesh -- in 
> which case you already know that something is wrong. The question then 
> is why that is so -- is your rhs correct, for example? 
>

After modified some of mistakes, now I have following result for pressure 
field.  
With same manufactured pressure solution...


  




<2x2 mesh> 





<4x4 mesh> 







For velocity field, I still have at least qualitatively good result as I 
posted previous comment.
Tho it does not converge with mesh refinement when it is compared to 
manufactured solution in L2_norm base. 


 
 

> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Assemble Righthand Side for vector-valued problem

[Wolfgang Bangerth]

On 02/08/2017 01:15 PM, Jaekwang Kim wrote:


I am using deal-ii to analyze Non-newtonian fluids, which varies viscosity.
For example, my viscosity eta is function of (u) field. I am solving
this with iteration method.. (Picard)
After coding, I am testing my code with method of manufacture.

For squared domain, I implemented all boundary condition with my assumed
solution and impose source term also to numerically generate my assumed
solution
however, the problem I got is also follow


While my velocity field is qualitatively similar to manufactured solution,
the pressure field is not following the exact form as seen below figure.
It is very strange that..I can approach u-soltution form while my
pressure does not


In addition , any mesh refinement would not make more accurate solution
now.

If you were me, what would you check more to find out what is problem?


I would make the problem simpler. Do you get the right solution if you 
have a constant viscosity? If you iterate, do you get the solution after 
one iteration?


I will also note that your manufactured solution is symmetric, but your 
computational solution is not. This provides you with a powerful way 
because you don't actually need to look at the convergence, it's enough 
to check that on a *coarse* mesh the solution is *symmetric*. It may be 
that solution is already non-symmetric on 1 cell, or on a 2x2 mesh -- in 
which case you already know that something is wrong. The question then 
is why that is so -- is your rhs correct, for example?


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Assemble Righthand Side for vector-valued problem

[Jaekwang Kim]

Thanks, I understand what you meant!!
I fixed it. 

Add up this.. can I ask your intuition for my another problem?

I am using deal-ii to analyze Non-newtonian fluids, which varies viscosity. 
For example, my viscosity eta is function of (u) field. I am solving this 
with iteration method.. (Picard) 
After coding, I am testing my code with method of manufacture. 

For squared domain, I implemented all boundary condition with my assumed 
solution and impose source term also to numerically generate my assumed 
solution 
however, the problem I got is also follow 


While my velocity field is qualitatively similar to manufactured solution, 
the pressure field is not following the exact form as seen below figure. 
It is very strange that..I can approach u-soltution form while my 
pressure does not


In addition , any mesh refinement would not make more accurate solution 
now. 

If you were me, what would you check more to find out what is problem?







2017년 2월 8일 수요일 오후 1시 40분 46초 UTC-6, Wolfgang Bangerth 님의 말:
>
> On 02/08/2017 12:35 PM, Jaekwang Kim wrote: 
> >*/_local_rhs(i) += fe_values.shape_value(i,q) *_/* 
> > 
> > */_rhs_values[q] *_/* 
> > 
> > */_fe_values.JxW(q);_/* 
> > 
>
> rhs_values[q] is a Vector. You need to say which component of it 
> you want to multiply with. (Hint: The correct component is 
> fe.system_to_component_index(i).first). 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Assemble Righthand Side for vector-valued problem

[Wolfgang Bangerth]

On 02/08/2017 12:35 PM, Jaekwang Kim wrote:

   */_local_rhs(i) += fe_values.shape_value(i,q) *_/*

*/_rhs_values[q] *_/*

*/_fe_values.JxW(q);_/*



rhs_values[q] is a Vector. You need to say which component of it 
you want to multiply with. (Hint: The correct component is 
fe.system_to_component_index(i).first).


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Assemble Righthand Side for vector-valued problem

[Jaekwang Kim]

Hi all, 

I was testing my code, with method of manufactured solution. 

I have a source term, which is vector value, defined as .. 


  template 

  class RightHandSide : public Function

  {

  public:

RightHandSide () : Function(dim+1) {}


virtual void vector_value (const Point ,

   Vector   ) const;


  };


  template 

  void

  RightHandSide::vector_value (const Point ,

Vector   ) const

  {

  

  double n = 0.5 ;

  

  double x =p[0]; double y=p[1];

  

  double f_x=2*x + 7*y - 4*Power(E,3*(Power(-0.5 + x,2) + Power(-0.5 + 
y,2)))*(-0.5 + n)*

  (1 + 2*Power(-0.5 + x,2))*(-0.5 + x)*

  (240. - 320.*x + 448.*Power(x,2) - 256.*Power(x,3) + 128.*Power(x,4) 
- 192.*y +

   320.*Power(y,2) - 256.*Power(y,3) + 128.*Power(y,4))*

  Power(1 + (Power(E,2*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*

 (64.*Power(0.75 - 1.*x + Power(x,2),2) + 64.*Power(0.75 - 
1.*y + Power(y,2),2)))

/2.,-1.5 + n) - 16.*Power(E,Power(-0.5 + x,2) + Power(-0.5 + y,2
))*(-0.5 + x)*

  Power(1 + (Power(E,2*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*

 (64.*Power(0.75 - 1.*x + Power(x,2),2) + 64.*Power(0.75 - 
1.*y + Power(y,2),2)))

/2.,-0.5 + n) - 8*Power(E,Power(-0.5 + x,2) + Power(-0.5 + y,2
))*

  (1 + 2*Power(-0.5 + x,2))*(-0.5 + x)*

  Power(1 + (Power(E,2*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*

 (64.*Power(0.75 - 1.*x + Power(x,2),2) + 64.*Power(0.75 - 
1.*y + Power(y,2),2)))

/2.,-0.5 + n);

  

  double f_y=7*x + 8.*Power(E,3*(Power(-0.5 + x,2) + Power(-0.5 + y,2
)))*(-0.5 + n)*(-0.5 + y)*

  (0.75 - 1.*y + Power(y,2))*(240. - 192.*x + 320.*Power(x,2) - 256.
*Power(x,3) +

  128.*Power(x,4) - 320.*y + 448.*Power(y,2) 
- 256.*Power(y,3) + 128.*Power(y,4))*

  Power(1 + (Power(E,2*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*

 (64.*Power(0.75 - 1.*x + Power(x,2),2) + 64.*Power(0.75 - 
1.*y + Power(y,2),2)))

/2.,-1.5 + n) + 16.*Power(E,Power(-0.5 + x,2) + Power(-0.5 + y,2
))*(-0.5 + y)*

  Power(1 + (Power(E,2*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*

 (64.*Power(0.75 - 1.*x + Power(x,2),2) + 64.*Power(0.75 - 
1.*y + Power(y,2),2)))

/2.,-0.5 + n) + 16.*Power(E,Power(-0.5 + x,2) + Power(-0.5 + y,2
))*(-0.5 + y)*

  (0.75 - 1.*y + Power(y,2))*Power(1 +

   (Power(E,2*(Power(-0.5 + x,2) + 
Power(-0.5 + y,2)))*

(64.*Power(0.75 - 1.*x + Power(x,2),
2) + 64.*Power(0.75 - 1.*y + Power(y,2),2)))

   /2.,-0.5 + n);

  

  

  values(0) =  f_x;

  values(1) =  f_y;

  values(2) =  0;

  }
 


and I wanted to assemble RHS in assemble_system function


 template   // Described in Cartesian Coordinate

  void ShearProblem::assemble_system ()

  {

system_matrix=0;

system_rhs=0;



  

  // Variable Coefficient

double shear_rate;

double viscosity;



const MappingQ mapping (degree);

QGauss   quadrature_formula(4*degree+1);


FEValues 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   local_matrix (dofs_per_cell, dofs_per_cell);

Vector   local_rhs (dofs_per_cell);


std::vector local_dof_indices (dofs_per_cell);


const RightHandSide  right_hand_side;

  

std::vector  rhs_values (n_q_points, Vector(dim+1));


const FEValuesExtractors::Vector velocities (0);

const FEValuesExtractors::Scalar xvel (0);

const FEValuesExtractors::Scalar pressure (dim);


   

std::vector > symgrad_phi_u (dofs_per_cell);

std::vector  div_phi_u   (dofs_per_cell);

std::vector  phi_p   (dofs_per_cell);


  

std::vector > local_previous_symgrad_phi_u 
(n_q_points);

std::vector  local_previous_solution (n_q_points);



typename DoFHandler::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;



right_hand_side.vector_value_list(fe_values.get_quadrature_points(),rhs_values);

fe_values[velocities].get_function_symmetric_gradients 
(previous_solution, local_previous_symgrad_phi_u);

  

  

for (unsigned int q=0; q