Implementing Robin (Dirichlet+Neumann) boundary condition
 

Hello
I have attached a description of the project I am working on.

I don't get my expected results, I doubt that my way of implementation of the bottom boundary condition might not be correct.

we already discussed your problem in other posts and I still do not understand how the number of BC.s you used can close your fourth-order PDE...
Then, what about a stability analysis of your discretization?

Yes, I mean we have two initial conditions for time and two boundary conditions at top and bottom, and periodic boundary condition in the xi direction, that's how it closes the equation.

I did von neuman analysis and it showed that for the time steps and grid size that I am using it is stable, no worries for that, even I tried crank nicolson method it gave the same results.

Accepting "per fides" that your problem is mathematically well posed, the fact that by decreasing the steps dx and dz the solution becomes wrong is possible if:

1) the dt is too large and produces a numerical instabilty. I have not checked carefully your discrete equation, since your multi-step method produces a 2x2 amplification matrix are you sure that both eigenvalues ensure a conditional stability? However, that seems not your case as you wrote that the solution is stable and dissipate the wave.

2) the discretized problem (equation + bc.s) is not consistent, that means your local truncation error does not vanish for vanishing mesh sizes.

Could you plot your numerical solutions and the analytical one at a certain time-step?

Here I have attached the time evolution of the solution and the exact analytical expression, as you see they perfectly mach, and the other plot shows the convergence as you can see the difference goes to almost zero 1e-6, although I should have plotted it logarithmically to show it better, but it well converges and follows that analytical solution. This is for the case without bottom corrugation.

1) from the plot I cannot see the shift of the wave from the position at t=0, could you superimpose the initial condition?

2) the plot is just for a value of w in one point, could you plot the global solution w(x,z)?

However, I searched in the literature but I did not find any PDE resembling that in your problem...:confused:

I agree with Prof Denaro. The splitting that we discussed before assumes that you can solve a forced wave equation for phi that is defined as Laplacian(w), extract w from a Poisson equation, and then compute the next timestep driving function for the wave equation from that solution. And, this may be a favorable splitting for the boundary conditions in the initial test problem, but not for more general BCs.

There are two related and unresolved issues...characteristic/information flow directions/domains of dependence...elliptic, parabolic, etc...and well-posedness. I have no experience with equations of this sort...and your reticence to tell us about the physical system it purports to represent doesn't offer me or the others here any impetus to develop it. Quite frankly, we can write down PDEs at random that probably have very interesting and convoluted behaviors, but we have better things to think about.

Having said all of that, let me just say that you can likely get to the bottom of at least the information flow issue by building a fully implicit time discretization for w. That would entail using an implicit time integration scheme (backward Euler, Adams Moulton, or implicit RK) and apply it to the first-order split system that we outlined in previous threads. The system may be block implicit in w and dw/dt. But, this will give you a finite difference formulation that does NOT assume any directional information flow. The next timestep information will depend completely on old time and new time information. This does assume that the system is an evolution-type equation...the domain of dependence for new time lies in the old time domain and includes only THE NEXT time information...not all future time as well. That is probably a safe assumption. And this should give you results that don't rely on SPATIAL direction assumptions. If the linear system that you build is not singular or indeterminate, that (at least) gives you an indication that the system is possibly well-posed.

I will tell you that using implicit timestepping will likely introduce diffusion into the system and your waves will be damped over time. There may be energy conserving time integration schemes (especially implicit RK) that may help with this.

https://drive.google.com/file/d/0B4b83ymv5MKZc3ZZWm5Qb01HY28/view?usp=sharing

@FMDenaro I have simulated the one with normal boundary conditions, and this is the video of it

@mprinkey
Thanks for the reply, if you have a look at my uploaded file on this post you see that I did not do any splitting and I solved the whole equation. Worth mentioning even with splitting, no matter in what way, the results are the same, upon similar resolutions.
I am not reticent in giving information I have mentioned all I know, all I was assigned by my course supervisor was to solve the PDE, I do not know more that what I have written.
Besides, I as well tried crank nicolson scheme, but still I had the limitation of spatial resolution, and it is weird, when I increase the spatial resolution, say more than 11 points per wave length the wave does not move at all, and gets damped away.

Yes, I see now, it appears nothing else that like rigid moving wave along the z direction ... it is like solving the equation dw/dt + c* dw/dz=0 at several positions along x with sinusoidal initial condition.

That increases my doubts ... can a second order in time PDE produce a single wave solution? should not be unstable (or presenting wiggles) the solution with central discretization for the case of a wave propagation?
As I wrote previously, despite my searching, I have not found nothing in literature like your PDE.

you mean a central difference might not be a good way of discretizing a wave like equation like this? since it cannot capture instability?

if you work with hyperbolic equations, the solution has a domain of dependence defined by the characteristic curves, the central discretization for a first order PDE is clearly a bad choice as it violates the domain of dependence. Conversely, if you have high order PDE you can have solution with two waves coming from different regions and the central discretization could get reasonable result.

To tell the true, I suppose that could be the particular choice of your initial condition to have the form of an eigenfunction of the operators in your PDE...

Then I will as well try with a forward and backward discretization (depending on closeness to each side) as well to see whether it fixes the problem that happens with increase in resolution.

no, this is not the issue...in case of first order hyperbolic equations, forward/backward stencils are automaticaly determined by the sign of the eigenvalues of the PDE, they are not determined by the limit of the domain.

Then it comes again the issue, that my equation is not any of the conventional PDE types. Then I should probably switch to finite element method for instance

no, finite element is just a different way to discretize the PDE, does not change the original nature of your PDE... what is more, the FE require to produce a suitable functional form of your PDE with proper shape function (in space). And the time integration is often nothing else that a FD discretization.

But a colleague has successfully with high resolution with finite element method has solved this PDE with no such issues that I have, and he has not examined the stability or has not found out the PDE type and etc, just has discretized and solved it.

The formulation in the pdf at the top of this thread is NOT implicit in time even though you are solving a linear system for the new time values of w. You used central differencing of the d^2/dt^2 operator, so the new time value is known explicitly by solving for phi_n+1 = laplacian(w_n+1). The fact that you have implicitly discretized the spatial operator and built a linear system does not make this an implicit time stepping scheme. So, this assumes some information flow in the system, just like the earlier splitting approach did.

My recommendations from before stand. Adams Moulton, backward Euler, or Implicit RK for the time advancement.

Yes in the pdf there is an explicit scheme, but I have also tried crank nicolson, but it had a similar behavior, limitations and problems described above.

Here is the source of these equations. See Appendix A.

http://arxiv.org/pdf/1604.07308v1.pdf

I had to find you at Berkeley and then find your advisor. You could have made this a lot easier for us if you told us this was some shallow wave model and gave us something to go on. I will look at this again, as you asked.

You said that the waves are damped, even in the flat bottom case as you resolve the spatial grid. Are you reducing the timestep proportionally? Also, are you solving the linear system to high tolerance?

The goal of the second problem seems to be the construction of the wave solution with a corrugated seabed. So rather than w = 0 at z = 0, it needs to be w = 0 at z = A sin(k x) where A is small compared to the water height. I don't see the value in trying to impose some complicated boundary condition while keeping z = 0. It is certainly possible and would be the direction one might take if pursuing an analytical or perturbation solution, but it discards one of the most important features of computational methods--the ability to have variable spaced grids/control volumes/elements. So, I will offer three approaches that should all work for the problem as defined. These are enumerated from EASY to HARD. Here is a visual aid for the first two approaches. The third, I mention only for the sake of completeness:

http://aeolustec.com/download/4thOrderWave.pdf

(1) Do a Cartesian-type grid that intersects the corrugated seabed by making the first line of points in the z direction have a z-coordinate = A sin(k x). The w value at that boundary is still set at 0. The 2nd derivative in the z-direction will need to be adjusted for the first computed row of values where z = dz because there will be two different delta values. The distance from that row to the seabed will NOT necessarily be dz. The formula for d^2/dz^2 will need to be replaced in all locations with this:

https://mathformeremortals.wordpress...-three-points/

This is fairly EASY to implement starting with what you have already done. It is really just modifying the z-derivatives formulae. This implementation will only be correct if the amplitude of the seabed corrugation is less than the scale of dz. Maybe A < 0.3 dz.

(2) Proper coordinate transformation of the corrugated seabed domain. The transformation would look like:

x' = x
z' = (z - A sin(k x)/(1 - A sin(k x))

So, at z' = 0, z = A sin(k x) and the other boundary conditions locations and values remain the same and you can discretize the x',z' system using constant values of dx' and dz'.. But the coordinate transformation does require substantial modification to the governing equations because the coordinate lines are no longer straight lines and orthogonal. See Equation 39 in this paper:

http://www.csun.edu/~lcaretto/me692/...formations.pdf

to see how to write the Laplacian in the transformed domain. This is an difficult process because it requires analytical differentiation and much more code to get right, but this is a typical structured-grid approach from finite difference codes of days gone by. This does have the merits of being valid for any reasonable choice of corrugation height. This is MODERATELY difficult to implement.

(3) Unstructured finite volume formulation...this is likely the enabling fact that allowed your colleague to succeed with finite elements. The seabed it meshed with unstructured finite volumes in much the same way one might for finite elements. The transport equations then are built by integrating over each control volume...it will lead to a lot of complicated and potentially critical non-orthogonal gradient computations at the faces to do those surface integrals. If you want to pursue this approach, you are likely better off starting with OpenFOAM and build on their meshing and unstructured formulation. Starting from scratch here will be HARD.

Thank you for the reply.

I will proceed with your first solution probably, which seems more straightforward than the others. Yes the other guys have done this with finite element method using unstructured grid at the bottom, seems that the finite difference scheme cannot handle this type of PDE quite well, and has limitations.

But, I got to know that the one even without corrugation, only runs for specific time step and resolution.
I decreased time step for instance, and decreased the grid spacing proportionally, but it did not work and soon got damped away. (proportional in the CFL sense)

The optimum value for grid resolution came out 11 points per wavelength (a peak and trough) of the initial wave, that I use from the analytical expression as initial condition.
Suppose I initialize with a wave with 2 peaks and troughs in both x and z direction with 22 points only in x and z direction. It works fine and waves move for long run without any damping, and solution converges fast and follows theoretical expression.
But if with the same resolution and time step and domain, I just make the initial wave, one with 1 peak and trough, wave never moves, and gets damped away soon, with this case even if I make the time step half it does not move again and damps away. Only if I change number of points to 11 it works fine again.

It is so weird how it behaves, and I have no idea for that. So that's why I never moved forward for corrugated bottom.

Thank you for your time.

Just some ideas according to the Michael's comments

1) Owing to the linearized model, the BC.s can be set at z=0 providing a perturbation velocity along the plane. That is quite common for linearized problem base on small pertubation. Thus, FD are quite simple to be used.

2) The 4th order PDE comes from a standard hyperbolic system. I suppose it should somehow respect the mathematical constraint for hyperbolic problems in term of BC.s

3) the solution cannot be optimal only on a grid... It looks like when the 1D linear wave equation is solved using the Courant number equal to 1 and produces the exact solution for those time and grid steps. But this is just a particular case...

P.S.: sorry to say that from the reading of Appendix A I still have trouble understanding what is done... The procedure starts from the Euler equation (2.1a) that is a balance of accelerations ([v]/[t]), then (after linearization) that produces an equation for the balance of terms dimensionally equivalent to [v]/([t][L])^2.
That is not clear to me at all...