How to solve 2 coupled PDEs with inner boundary condition?

I want to solve 2 coupled PDEs with inner region boundary conditions.

Diffuse equation

There are 2 kinds of particles in my system, $$p$$ and $$a$$. Their diffusion equations are as follows: $$\hat{D}=\begin{pmatrix} D_T& 0& 0 \\ 0 & D_T &0\\0 &0 &D_r \end{pmatrix}\\ \quad$$

$$D_T$$: transportation coefficient

$$D_r$$: rotation coefficient

A:

$$a = a (x,y,\theta)$$

$$\frac{\partial a }{\partial t}=\nabla\cdot \vec{j_A}$$

$$\vec{j_A}=\hat{D}\nabla a -\vec{v} a$$

$$\vec{v}=(v_0 \cos\theta,v_0 \sin\theta,0)$$

P can be seen as A whose $$v_0=0$$:

$$p = p (x,y,\theta)$$

$$\frac{\partial p }{\partial t}=\nabla\cdot \vec{j_P}$$

$$\vec{j_P}=\hat{D}\nabla p$$

Boundary condition

P will change to A at circle 1, and A will transfer to P at circle 2.

$$\hat{n}=\vec{r}/|r|$$

On the wall $$r=R$$

$$\hat{n}\cdot \nabla p (x,y,\theta)=0$$

On the circle2 ($$r=\sqrt{x^2+y^2}=r_2$$):

$$\hat{n}\cdot (\nabla p (x,y,\theta)-(\vec{j_A}+\vec{j_P}))=0$$

$$a =0$$

On the circle1 ($$r=\sqrt{x^2+y^2}=r_1$$):

$$\hat{n}\cdot (\nabla a (x,y,\theta)-(\vec{j_A}+\vec{j_P}))=0$$

$$p =0$$

On the center($$r=\sqrt{x^2+y^2}=0$$):

$$\rho(x,y,\theta)=1$$

code

I want to get a stable solution.

I know how to solve single PDE using Mathematica, but i don't know how to combine $$a$$ and $$p$$.

This program assumes that $$a(x,y,\theta)==1$$ on circle 1, in fact it should be solved by coupling $$a(x,y,\theta)$$ and $$p(x,y,\theta)$$.

{R, r1, r2, v0, Difft, Diffr} = {29.5, 15, 25, 0.1, 1/29.5^2, 0.01};

Needs["NDSolveFEM"]
\[CapitalOmega] =
RegionDifference[Cylinder[{{0, 0, 0}, {0, 0, 2 \[Pi]}}, r2/R],
Cylinder[{{0, 0, 0}, {0, 0, 2 \[Pi]}}, r1/R]];
cylinderMesh = ToElementMesh[\[CapitalOmega], MaxCellMeasure -> 0.001];
c = {{Difft, 0, 0}, {0, Difft, 0}, {0, 0, Diffr}};
v = {v0/R  Cos[\[Theta]], v0/R  Sin[\[Theta]], 0};

With[{a = a[x, y, \[Theta]]},
op = 0 ==
Div[c . Grad[a, {x, y, \[Theta]}] - v a, {x, y, \[Theta]}]];

res = NDSolveValue[{op,
DirichletCondition[a[x, y, \[Theta]] == 1, x^2 + y^2 == (r1/R)^2],
DirichletCondition[a[x, y, \[Theta]] == 0, x^2 + y^2 == (r2/R)^2]},
a, {x, y, \[Theta]} \[Element] cylinderMesh]


I guess that I should decompose the disk into 3 regions($$rr_2$$) and respectively solved them.

How shall I solve it? thank you!