Re: [Getfem-users] Pure Advection Stationary Problem

[Yves Renard]

Dear Simone Di Gregorio,

You are of course right saying that the approximation of transport equation with a finite element method may induce some instabilities. This is true except in certain cases where the coercivity is preserved (for instance with a Dirichlet condition at the right end point and a with monotone velocity field). Adding the viscous term  solves the problem of course for beta large enough as you explain.

Your formulation seems to be perfectly correct to me. The only thing is that it seems to me that it is transposed. I would instead write
asm_real_or_complex_1_param(M,mim,mf,mf_data,V,rg, "a=data$1(#2);" "M$1(#1,#1)+=comp(Grad(#1).Base(#1).Base(#2))(:,i,:,i,p).a(p);")
because there is a test function for each line of the linear system (not for each column).

Of course, this is the _expression_ for a 1D (straight) domain.


A suggestion : you should pass to the high level generic assembly which is more easy to work with (and is now faster).

Best regards,

Yves


Le 16/03/2017 à 10:11, Simone Di Gregorio a écrit :

Dear GetFem users,

I'm dealing with a pure stationary advection problem in 1D straight single branch domain. So the equation is:

 d( v(s)c(s))/ds = 0 

with v(s) velocity field (NOT constant, so it varies with s) and c(s) solute concentration, while s is the arc length of the vessel. The boundary conditions are Dirichlet B.C. in inlet of the vessel

The weak form is:

(\partial_s  \phi, v c)_{\Lambda} + v(\Gamma _{out})c(\Gamma _{out})\phi(\Gamma _{out}) = 0

where \phi is the piecewise linear continuous test function, \Lambda si the domain, \Gamma _{out} is the outlet boundary.

So I build the advection matrix with:

asm_real_or_complex_1_param(M,mim,mf,mf_data,V,rg, "a=data$1(#2);" "M$1(#1,#1)+=comp(Base(#1).vGrad(#1).Base(#2))(:,:,i,j,p).a(p);")

Then I subtract to the last element of the matrix (in position N,N) the value of v(\Gamma _{out}).

As expected the solution is unstable, so i added an artificial diffusion in the following way:

+D(\partial_s  \phi, \partial_s  c)_{\Lambda} 

with D= beta*max(v(s))*h, beta large enough and h dimension of the element, and imposing dc/ds=0 over the boundary. I assembled this with 

asm_real_or_complex_1_param(A,mim,mf,mf_data,D,rg, "a=data$1(#2);" "M$1(#1,#1)+=comp(vGrad(#1).vGrad(#1).Base(#2))(:,i,j,:,i,j,p).a(p);")

The solution becomes stable.

The problem is that:

-the solution is always the same for very different function of v(s):

-if D is not very high in the first elements of the vessel, the solution is constant, whereas it varies in the last element(this happens because add v(\Gamma _{out})

otherwise the solution would be costant everywhere)

-if D is quite high the solution is a simple exponential, but is not sensitive to variation of v

From an analytical solution I expect that v(s)c(s) remain constant so c(s) should vary along the vessel dependently on value of v(s).

So it's like the problem doesn't feel the variation of velocity and I can't understand what i'm doing wrong.

Thanking you in advance, i hope i was clear enough in my explanation

Best Regards

Simone Di Gregorio 

_______________________________________________
Getfem-users mailing list
address@hidden
https://mail.gna.org/listinfo/getfem-users

-- 

  Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
  Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE
  http://math.univ-lyon1.fr/~renard

---------