Time-dependent pde with periodic boundary conditions - NDSolve gives incorrect results

I'm trying to solve very simple differential equation with perdiodic boundary conditions and non-constant initial condition, however the results I'm getting don't seem to be correct. Here's my code:

size = 10 (* size of the domain *)
tMax = 20000000 (* time of the experiment *)
pde = D[h[x, y, t], t] == -0.7
BConds = {h[0, y, t] == h[size, y, t], h[x, 0, t] == h[x, size, t]} (* periodic boundary conditions *)
topography[x1_, x2_] = Piecewise[{{400./size*x1 - 100., size/4. < x1 < size/2.}, {300. - 400./size*x1, size/2. < x1 < 3.*size/4.}, {100., x1 == size/2.}, {0, x1 <= size/4.}, {0, x2 >= 3.*size/4.}}]
ICond = h[x, y, 0] == topography[x, y] (* initial condition as a function of x *)
sol = NDSolve[Join[{pde}, BConds, {ICond}], h, {x, 0, size}, {y, 0, size}, {t, 0, tMax}]
fun[x_, y_, t_] := h[x, y, t] /. sol[[1]]


I wound expect, because of how the initial condition is specified, the solution to be symmetrical with respect to size/2. That's true for e.g. x=size/4 and x=3*size/4, meaning .fun[size/4,0,0]=fun[3*size/4,0,0].

However that's not true for x=size/3 and x=2*size/3, for which the solution in my opinion also should have the same values, but it doesn't:

fun[size/3, 0, 0]
30.8642
fun[2*size/3, 0, 0]
31.6872


Why the values aren't the same?

The other problem is when I try to plot the solution:

Manipulate[ContourPlot[fun[x, y, t], {x, 0, 10}, {y, 0, 10}, ColorFunction -> (ColorData["Rainbow"][1 #] &), PlotLegends -> Automatic], {t, 0, tMax}]


The plot, for t=0, displays regions of fun[x,y,t]<0 even though the minimum value for the initial condition is 0. I would guess it's because of how the solution is interpolated, but could someone explain it to me in more details?

• I tried to use Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> "FiniteElement"}} and it gives much better results, but when I apply periodic boundary conditions the execution fails saying that FEM doesn't support that. Dec 2, 2015 at 21:15

sol = NDSolve[Join[{pde}, BConds, {ICond}],