Discretisation of a PDE in Lagrangian coordinates

[hunt_mat]

TL;DR

Sense check for a discretisation scheme for 2D Lagrangian scheme

I am writing a 2D hydrocode in Lagrangian co-ordinates. I have never done this before, so I am completely clueless as to what I'm doing. I have a route as to what I want to do, but I don't know if this makes sense or not. I've gone from Eulerian to Lagrangian co-ordinates using the Piola identities and then using the finite volume method to get a discrete model. One of the equations that I've been looking at is:
<br /> \frac{\partial D}{\partial X}-\frac{\partial C}{\partial Y}=0<br />
Choosing a rectangle with vertices (X,Y),(X+\delta X,Y),(X+\delta X,Y+\delta Y),(X+Y+\delta Y) to integrate over, and using Green's theorem in the plane, one is able to obtain:
<br /> \int_{X}^{X+\delta X}C(s,Y+\delta Y)+C(s,Y)ds+\int_{Y}^{Y+\delta Y}D(X,s)+D(X+\delta X,s)ds=0<br />
If I then want to discretise this, I would use the approximation f(x)+f(x+\delta x/2)\approx f(x+\delta x/2), and using the midpoint rule for integrals one then obtains:
<br /> C(X+\delta X/2,Y+\delta Y/2)\delta X+D(X+\delta X/2,Y+\delta Y/2)\delta Y\approx 0<br />
Does this sound reasonable to you?

 

[pasmith]

I tihnk you have some sign errors.

You are using <br /> \iint_R \frac{\partial D}{\partial X} - \frac{\partial C}{\partial Y}\,dX\,dY = \oint_{\partial R} (C\mathbf{e}_X + D\mathbf{e}_Y)\cdot d\mathbf{X} where the boundary is traversed counter-clockwise, so that <br /> \oint_{\partial R} = \int_{(X,Y)}^{(X + \delta X,Y)} + \int_{(X + \delta X, Y)}^{(X + \delta X, Y + \delta Y)} + \int_{(X + \delta X. Y + \delta Y)}^{(X,Y + \delta Y)} + \int_{(X,Y+ \delta Y)}^{(X,Y)} and hence <br /> \oint_{\partial R} (C\mathbf{e}_X + D\mathbf{e}_Y)\cdot d\mathbf{X} = \int_{X}^{X + \delta X} C(s,Y) - C(s,Y + \delta Y)\,ds + \int_Y^{Y + \delta Y} D(X + \delta X,s) - D(X, s)\,ds.

How does your discretisation work? Are you storing the values of C and D at the midpoint of the cell or at the vertices of the cell? If at the midpoint, then the corners of the cell are at (X_n \pm \frac12\delta X, Y_m \pm \frac12\delta Y) and you should approximate values on the cell boundary by interpolation: \begin{split}<br /> C(X_n \pm \tfrac12 \delta X, Y_m) &amp;= \frac{C_{n\pm 1,m} + C_{n,m}}{2} \\<br /> D(X_n, Y_m \pm \tfrac12 \delta Y) &amp;= \frac{D_{n,m\pm 1} + D_{n,m}}2.\end{split}

 

[hunt_mat]

The way you wrote Green's theorem is very odd. As you wrote it, it seems confused. It might be easier if you just wrote the normal \hat{\mathbf{n}}, and the flux, \mathbf{F}=(D,-C). It's I think that you've made an error with the signs of your normals.

I'm using the finite volume method, so I'm, storing things at the cell centres. The cell vertices are as I gave them.

 

[pasmith]

If you integrate <br /> \frac{\partial D}{\partial X} - \frac{\partial C}{\partial Y} over a rectangle, then applying the FTC to each term separately reduces it to \begin{split}<br /> \int_Y^{Y + \delta Y} \int_X^{X+\delta X} \frac{\partial D}{ \partial X}\,dx\,dy - \int_X^{X+\delta X} \int_Y^{Y + \delta Y} \frac{\partial C}{ \partial Y}\,dy\,dx \\<br /> = <br /> \int_Y^{Y + \delta Y} D(X + \delta X, s) - D(X,s)\,ds - \int_X^{X + \delta X} C(s, Y + \delta Y) - C(s,Y)\,ds\end{split} which is consistent with my expansion, but not with yours.

 

[hunt_mat]

Okay.