How to solve 2 coupled PDEs with inner boundary condition?

Question

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["NDSolve`FEM`"]
\[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($r<r_1,r_1<r<r_2,r>r_2$) and respectively solved them.

How shall I solve it? thank you!