Hi All, I figured out the problem with my diffusion equation. Now the next step for me is to test an advection-diffusion equation.
I am using a stabilized formulation i.e. the weak form looks like
(w, theta_dot) + (tau*{u}dot{grad w}, theta_dot) + (w, {u}dot{grad
theta}) + (tau*{u}dot{grad w}, {u}dot{grad theta}) + (k {grad w},
{grad theta}) = (w, f)
I am trying to solve rotating cone problem (a simple advection of a
cone in a rotating fluid)
{u} = (y, -x)
domain = [-1,1]x[-1,1]
k = 1e-6
cone {center: (0.5, 0), radius: 0.2, height: 1.0}
I see the following:
1. when no stabilization is added, i.e. tau = 0, I get a rotating
cone, however, with spurious oscillations. This is the correct behavior
2. when I use tau = h_K/ (2*||{u}||), the solution is completely wrong
I have attached the relevant parts of assemble matrix. Does my
implementation look correct.
Thanks,
Badri
NOTE:
h_K = element diameter, defined as sqrt(detJ*16/PI), where detJ is the
determinant of element Jacobi matrix
++++++++++++++++++++++++++++
double detJ = 0.0;
for(uint q=0; q<numItgPts; ++q)
{
detJ += fe_values.JxW(q);
}
double h = sqrt(16.0*detJ/pi);
for(uint q=0; q<numItgPts; ++q)
{
double tau = h/(2.0*velocity_values[q].norm());
//tau = 0.0;
for(uint a=0; a<nbf; ++a)
{
for(uint b=0; b<nbf; ++b)
{
Ae(a,b) +=
(
// M
(
fe_values.shape_value(a,q) *
fe_values.shape_value(b,q)
) +
// M_delta
(
tau *
(
velocity_values[q] *
fe_values.shape_grad(a,q)
) *
fe_values.shape_value(b,q)
) +
time_step *
// N
(
fe_values.shape_value(a,q) *
velocity_values[q] *
fe_values.shape_grad(b,q)
) +
// N_delta
time_step *
(
tau *
(
velocity_values[q] *
fe_values.shape_grad(a,q)
) *
(
velocity_values[q] *
fe_values.shape_grad(b,q)
)
) +
// K
time_step *
(
diffusion_coeff_values[q] *
fe_values.shape_grad(a,q) *
fe_values.shape_grad(b,q)
)
)* fe_values.JxW(q);
}
be(a) +=
(
fe_values.shape_value(a,q) *
(
time_step *
force_vector_values[q] +
solution_prev_values[q]
) +
(
tau * velocity_values[q] *
fe_values.shape_grad(a,q) *
time_step *
force_vector_values[q]
)
)* fe_values.JxW(q);
_______________________________________________
