The initial condition was inconsistent and needed to exclude the boundary:

iC = Rationalize[
   P[x1, x2,0] == 
    Piecewise[{{1/
        Area[\[CapitalOmega]], (RegionMember[
           RegionBoundary[\[CapitalOmega]], {x1, x2}] == False \[And] 
         RegionMember[\[CapitalOmega], {x1, x2}])}}, 0], 0];

I know nothing about Kolmogorov equation. Your initial condition for P is a uniform value of 1/((Sqrt[3] \[Eta]^2)/4) over the triangle. Hence at t==0 all the spacial derivatives vanish. Now the initial time derivative of P is positive everywhere on the triangle:

fwrdKol /. t -> 0 /.
 {P[x1, x2, 0] -> 1/((Sqrt[3]   \[Eta]^2)/4),
  Derivative[_, _, 0][P][__] :> 0}

No surprise that the integral of P increases at the start.

An unrelated question: what is the purpose of all those Rationalize? Why don't you start with exact values for the parameters?

It seems reasonable that the mesh needs to be very fine only near the boundary, where the initial datum is discontinuous. You may try the option MeshRefinementFunction.

The purpose to make the input exact instead of floating-point is good. In your case it could be done more easily by giving exact parameters at the start. Instead of floating point values, as in

\[Eta] = 5.; 
\[Kappa] = .75; 
Tmax = 5.;

I would suggest

\[Eta] = 5; 
\[Kappa] = 3/4; 
Tmax = 5;

This way there is no need for Rationalize. Even more: with your original floating-point \[Eta] = 5. the triangle definition

\[CapitalOmega] = 
  Polygon[Rationalize[{{0, 0}, {\[Eta], 
      0}, {\[Eta]/2, (\[Eta]   Sqrt[3])/2}}, 
    0]]

actually decreases the precision, because \[Eta] Sqrt[3])/2 becomes first floating-point, and then Rationalize freezes the approximation instead of the exact value.

I would suggest to discretize the triangle so that the elements are smaller near the boundary:

\[CapitalOmega]Discretized = 
 DiscretizeRegion[\[CapitalOmega], 
  MeshRefinementFunction -> 
   Function[{vertices, area}, 
    area > 0.1 (RegionDistance[RegionBoundary[\[CapitalOmega]], 
        Mean[vertices]])]]

Instead of 0.1 you can use a smaller parameter, with longer computing time and hopefully better precision. Try using \[CapitalOmega]Discretized instead of \[CapitalOmega].

The time derivative of P is strongly negative near the boundary, enough to counteract the positive value inside:

Plot3D[Derivative[0, 0, 1][Psol][x1, x2, .01],
 Element[{x1, x2}, \[CapitalOmega]Discretized]]

Update 2: My formulation of the differential equation matches that produced automatically by Mathematica

eqn = ItoProcess[{\[DifferentialD]x1[
      t] == \[Kappa] (x1opt - 
        x1[t])  \[DifferentialD]t + \[DifferentialD]w1[
       t], \[DifferentialD]x2[
      t] == \[Kappa] (x2opt - 
        x2[t])   \[DifferentialD]t + \[DifferentialD]w2[t]}, {x1[t], 
   x2[t]}, {{x1, x2}, {x10, x20}}, 
  t, {w1 \[Distributed] WienerProcess[], 
   w2 \[Distributed] WienerProcess[]}]

eqn["KolmogorovForwardEquation"] // TraditionalForm 

Hi Cameron, I agree with you that your formulation of the problem and your PDEs look correct. If you plot your solution over time at one of the boundary points such as P(0,0,t), it looks like Mathematica is setting the boundary conditions to be P(0,0,t)*e^(-t) with the coefficient (-1) being exact. This is quite annoying, but you can remedy it by explicitly setting the boundary condition to P(0,0,t)*e^(-10000*t).

bC = Rationalize[
   DirichletCondition[
    P[x1, x2, t] == Exp[-10000  t] 1/((Sqrt[3]   \[Eta]^2)/4), True], 
   0]; 

In the study of linear PDEs this discontinuity (setting P(x,y,0)=const. on the domain and P(x,y,t)=0 on the boundary) isn't much of an issue and it's standard to do, so it's frustrating that Mathematica makes such a big deal about it! At least with this fix applied you can see that the derivative of the total probability mass approaches -infinity near t=0, so that hints at some potential numerical difficulties.