I am trying to solve a differential equation by `NDSlove` for $h(x,t)$. It reads
$$h_t=h_{xx}-V_h-\lambda(t)$$
where $V_h$ is a given function of $h(x,t)$ denoted by `vdh[x_,t_]` in my code, $\lambda(t)$ is a time-dependent parameter which is determined by a definite integration.

$$\lambda(t)=-\frac{1}{L}\int_0^L V_h dx$$
In addition, $h(x,t)$ subjects to simple periodic boundary condition and initial condition, which is defined on the periodic interval $[0,L]$, see my code.


First, Defining constants, function `vdh[x_,t_]` and $\lambda(t)$

    ClearAll[x, t]
    amp = 5787/1000;
    tmax = 1000;
    L = 20;
    vdh[x_?NumericQ, t_?NumericQ] := 1/(1 + h[x, t]) - h[x, t] + Log[h[x, t]/(1 + h[x, t])];
    \[Lambda][t_Real] := -(1/L)*NIntegrate[vdh[x, t], {x, 0, L}];
   
Then, the main part 

     myfunel = First[h /. NDSolve[{
     D[h[x, t], t] == D[h[x, t], {x, 2}] - vdh[x, t] - \[Lambda][t],
     h[0, t] == h[L, t],
     Derivative[1, 0][h][0, t] == Derivative[1, 0][h][L, t],
     Derivative[2, 0][h][0, t] == Derivative[2, 0][h][L, t],
     h[x, 0] == 1 + 1/Sqrt[2*\[Pi]]*Exp[-((x - 10)^2/2)],
     WhenEvent[h[L/2, t] >= amp, "StopIntegration"]},
     h, {x, 0, L}, {t, 0, tmax},
     Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
     "MinPoints" -> 101, "MaxPoints" -> 101, 
     "DifferenceOrder" -> 4}}, AccuracyGoal -> 20, 
     WorkingPrecision -> MachinePrecision, StepMonitor :> Print[t]]]

Which dose not work. `NDSolve` prompts a mass of errors, such as:

>NIntegrate::inumr: "The integrand 1-h[x,y,t]-h[x,y,t]/(1+h[x,y,t])+Log[h[x,y,t]/(1+h[x,y,t])] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,20},{0,20}}."

**Note 1**: In every time step, I want to determine $\lambda(t) $ numerically using `\[Lambda][t_Real]`, where lies the problem. It is well known that when `NDSolve` uses `(N)Integrate` internally, we will suffer from such kind of difficult very often. However, the equation I want to solve contain parameter $\lambda(t)$ found by integrating a function of the solution of my `NDSlove`. Fortunately, the `NIntegrate` does not involve `t` variable, which can be regarded as the average value of ${dV \over dh}$. I believe this point causes these problems. I'm not sure whether I have some other mistakes.

**Note 2**: I also need to see the evolution of the following function:

    f[t_] := NIntegrate[h[x, t]^2 + 1/2*(D[h[x, t], x])^2, {x, 0, L}]
    Plot[Evaluate[f[t]/.myfunel], {t, 0, tstop}, PlotRange -> All]

where `tstop` is the stop moment of the main part. I have used `WhenEvent` to stop my integration, but how to get the stop time in my main part rather than check `StepMonitor :> Print[t]` after its stop?

**Note 3**: I also tried `LaplaceTransform`, but it is obvious that this equation does nto has a *closed form solution*, like [this problem][1]

Any ideas? Also any suggestions for code improvements are more than welcome!

  [1]: http://mathematica.stackexchange.com/questions/24626/how-to-plot-and-solve-the-numerical-solution-of-a-integro-differential-equation