Is my weak form right?

Asked by cutejeff

To solve nonlinear PDE

du/dt+3*lmbda*u^2*du/dx*(sina-cosa*du/dx)-lmbda*u^3*cosa*d2u/dx2=0

is the following weak form right?

L = v*u*dx-v*u0*dx+dt*lmbda*v*3*u*u*u.dx(0)*(sin(alpha)-cos(alpha)*u.dx(0))*dx+dt*lmbda*u*u*u*cos(alpha)*u.dx(0)*v.dx(0)*dx
a = derivative(L, u, du)

I am sure the rest code for newton's method and time steps should be right. But I got the wrong result. I notice that there is a
(du/dx)^2 in the 2nd term. I guess it may cause some problems. Because if I get rid of the (du/dx)^2, the result is exactly same with my finite difference code (which should be right), If I keep it there, it's different.

Anyone know why? Thank you.

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
cutejeff
Solved:
Last query:
Last reply:
Revision history for this message
Garth Wells (garth-wells) said :
#1

What do you mean by 'wrong'? Does it fail to converge or do you get an unexpected result?

Revision history for this message
cutejeff (illw84u) said :
#2

Thank you for your help.

The PDE was derived using conservation of mass, but the results from Fenics showed the mass actually kept growing with time. Because my finite difference code doesn't grow, I consider my finite element code got some problems somewhere.

Interestingly, if I get rid of the (du/dx)^2 term in both Fenics and my finite difference code, the results match each other. Because I haven't seen any examples on manipulating the power of derivatives, I think it might be the cause. I don't think I can take integration by parts for v*(du/dx)^2*dx, I don't know what to do with it, except keep it in this way.

I also tried the convergence study, it's getting a little better if I improve the mash refinement and the polynomial degree. But it doesn't change much, the mass still grows.

Revision history for this message
Garth Wells (garth-wells) said :
#3

Pay special attention to the boundary conditions, especially the expression that you're implicitly setting equal to zero when deriving your weak from. Make sure that that term is the mass flux, otherwise you may have mass coming in across the boundary.

Revision history for this message
Johan Hake (johan-hake) said :
#4

On Thursday April 29 2010 08:57:52 cutejeff wrote:
> Question #107819 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/107819
>
> Status: Answered => Open
>
> cutejeff is still having a problem:
> Thank you for your help.
>
> The PDE was derived using conservation of mass, but the results from
> Fenics showed the mass actually kept growing with time. Because my
> finite difference code doesn't grow, I consider my finite element code
> got some problems somewhere.
>
> Interestingly, if I get rid of the (du/dx)^2 term in both Fenics and my
> finite difference code, the results match each other. Because I haven't
> seen any examples on manipulating the power of derivatives, I think it
> might be the cause. I don't think I can take integration by parts for
> v*(du/dx)^2*dx, I don't know what to do with it, except keep it in this
> way.

Are you having (du/dx)^2 or d2u/dx2 as you state in your original question?

If it is the latter we do integrate that by part in several demos, Poisson for
example. But you have u^3*d2u/dx2 which makes the integration by parts non
trivial.

You need to get that correct in addition to Garth's suggestions, I think.

Johan

Revision history for this message
cutejeff (illw84u) said :
#5

Thank you for both of your suggestions.
Dr. Wells: you meant the concomitant term from integration by parts might not be zero, but I implicitly set to zero. Do I understand it correctly?

 What I did for the B.C. was that I made the domain big enough and set 0 for B.C. at both ends. I understand u^3*d2u/dx2 need to be taken integration by parts, but according to my BC setting, the concomitant term is indeed zero.

Thank you again for your time to read my problem. Let me describe it in details.

I want to use this evolution equation to simulate the thin film flow on an incline. The PDE is derived using 2D Navier-Stokes equation, coupled with continuity equation and kinematic condition at free surface. See Acheson's book for detailed derivations.

Acheson, D.J., Elementary Fluid Dynamics. Oxford University Press, 1990.

In this 1D PDE,

du/dt+3*lmbda*u^2*du/dx*(sina-cosa*du/dx)-lmbda*u^3*cosa*d2u/dx2=0

u is the height of the free surface. I use a parabola “cap” as the initial condition
u=0.3[1-〖(x-2)〗^2] for1<x<3
u=0 for the other domain

To make the domain wide enough that in a short time period, the flow doesn't cross the boundary. I set the domain mesh = Interval(1000,0,16), at least in 12 sec, the flow won't touch the boundaries. Hence, u=0 at both ends all the time, which is also my BC.

According to the derivations, if I integrate u over the domain, I should get the mass of the fluids, it shouldn't change with time due to the conservation of mass. But the fact is that it grows for my Fenics code. I couldn't find any problem in my model so far.

I just start to learn the finite element processes, not very good at it. But I really have lots interests on this Fenics package. Any helps will be appreicated.

Revision history for this message
Garth Wells (garth-wells) said :
#6

On 30/04/10 17:58, cutejeff wrote:
> Question #107819 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/107819
>
> Status: Answered => Open
>
> cutejeff is still having a problem:
> Thank you for both of your suggestions.
> Dr. Wells: you meant the concomitant term from integration by parts might not be zero, but I implicitly set to zero. Do I understand it correctly?
>

Yes.

Garth

> What I did for the B.C. was that I made the domain big enough and set 0
> for B.C. at both ends. I understand u^3*d2u/dx2 need to be taken
> integration by parts, but according to my BC setting, the concomitant
> term is indeed zero.
>
> Thank you again for your time to read my problem. Let me describe it in
> details.
>
> I want to use this evolution equation to simulate the thin film flow on
> an incline. The PDE is derived using 2D Navier-Stokes equation, coupled
> with continuity equation and kinematic condition at free surface. See
> Acheson's book for detailed derivations.
>
> Acheson, D.J., Elementary Fluid Dynamics. Oxford University Press, 1990.
>
> In this 1D PDE,
>
> du/dt+3*lmbda*u^2*du/dx*(sina-cosa*du/dx)-lmbda*u^3*cosa*d2u/dx2=0
>
> u is the height of the free surface. I use a parabola “cap” as the initial condition
> u=0.3[1-〖(x-2)〗^2] for1<x<3
> u=0 for the other domain
>
> To make the domain wide enough that in a short time period, the flow
> doesn't cross the boundary. I set the domain mesh = Interval(1000,0,16),
> at least in 12 sec, the flow won't touch the boundaries. Hence, u=0 at
> both ends all the time, which is also my BC.
>
> According to the derivations, if I integrate u over the domain, I should
> get the mass of the fluids, it shouldn't change with time due to the
> conservation of mass. But the fact is that it grows for my Fenics code.
> I couldn't find any problem in my model so far.
>
> I just start to learn the finite element processes, not very good at it.
> But I really have lots interests on this Fenics package. Any helps will
> be appreicated.
>

Revision history for this message
cutejeff (illw84u) said :
#7

Thank for the helps.

As I said, the boundary conditions setting forced the the concomitant term from integration by parts be zero, so there shouldn't be any problem on B.Cs.

I still couldn't find where the problem is. I appreciate the suggestions. Any other thoughts about the mistake I could make?

Thank you.

Revision history for this message
Johan Hake (johan-hake) said :
#8

Have you figured out the integration by part?

Johan

On Monday May 3 2010 07:28:06 cutejeff wrote:
> Question #107819 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/107819
>
> cutejeff posted a new comment:
> Thank for the helps.
>
> As I said, the boundary conditions setting forced the the concomitant
> term from integration by parts be zero, so there shouldn't be any
> problem on B.Cs.
>
> I still couldn't find where the problem is. I appreciate the
> suggestions. Any other thoughts about the mistake I could make?
>
> Thank you.

Revision history for this message
cutejeff (illw84u) said :
#9

Johan:
Thank you for your reply.
For you last question, I have both (du/dx)^2 and d2u/dx2 terms in my PDE. I don't do integration by part for (du/dx)^2, I just make it multiply test function v. For d2u/dx2 term, I take integration by part as shown below,

d2u/dx2*v*dx=-du/dx*dv/dx*dx+du/dx*v@BC

Since I set u=0 at the boundaries, the concomitant term du/dx*v@BC is zero.

I think the FE processes on the power of the derivative may have some problem. I will open a question and make an easier example to describe it.

I close this one. Please refer to my new question. I appreciate the helps.