Setting different B.C. for different fields in NDSolve

Question

My question regards construction of a PDE solver for several fields, with different type of boundary conditions for several fields.

Say I have the following set of equations:

$$ \partial_{t}\rho+v\partial_{x}\rho=-\rho\partial_{x}v$$

$$\partial_{t}\chi+v\partial_{x}\chi=g_{1}\left(\frac{f}{\rho},\chi\right)$$

$$\partial_{t}f=\rho\left[1+\partial_{x}v-g_{2}\left(\frac{f}{\rho},\chi\right)\right]-f\partial_{x}v$$

$$\partial_{xx}v=\frac{1}{f-\rho}\left[\partial_{x}\rho\left(\partial_{x}v+1\right)-\partial_{x}\left(\rho\cdot g_{2}\left(\frac{f}{\rho},\chi\right)\right)\right]$$

where $\rho\equiv\rho\left(x,t\right)$, $\chi\equiv\chi\left(x,t\right)$,$v\equiv v\left(x,t\right)$ and $f\equiv f\left(t\right)$. I have here two PDE's (the first two equations), and the two other equations are consistent and ensure that $f$ indeed does not depend on space, and varies in time only.

Please note the following problems with this system of equations:

1. The functions $g_1$ and $g_2$ are exponential - this makes it hard to implement with the FEM package (as suggested here ) .

2. This is a set of a differential-algebraic equations (I am using the "IndexReduction" method to solve the system).

3. As explained, the last equation is derived from the third one (taking spatial derivative, which is very much like this discussion ).

In Mathematica I allow for $f$ to vary also in space, but the solution of the last equation ensures that $f$ is spatially independent. I impose periodic B.C. for all the fields (note that even though in principle I don't need to impose B.C. for $f$, because I let $f\equiv f\left(x,t\right)$ for simulation purposes, I have to put B.C. for this field as well).

I would like to impose periodic boundary conditions for $\rho,\chi$ and $f$, but demand that $v$ vanish at the boundaries. This makes much more sense for the system I am dealing with, and I also suspect that the $v$ dynamics with periodic B.C. generate some convergence problems.

Naive NDSolve with the mixed boundary conditions yields the following error message:

NDSolve`ProcessEquations::pdnbc: Mixed periodic and Dirichlet/Neumann boundary conditions are not allowed.

Is there a way to manipulate the boundary conditions (maybe using NDSolve's ProcessEquations function)?

My code (with periodic B.C. for all fields) looks as follows:

α = 1; β = 2; γ = 2;
g[ζ_, ξ_] = E^(-(1/ξ)) (E^ (-1 - α ζ) + E^ (-1 + α ζ)) (1 - β/ζ); 

eqns[S_, v_, f_, ξ_] := {
   D[S, t] + v D[S, x] == -D[v, x] S,
   D[ξ, t] + v D[ξ, x] == f/S  g[f/S, ξ] (γ - ξ), 
   D[v, x, x] == 1/(f - S) (D[v, x] D[S, x] + D[S, x] - D[S g[f/S, ξ], x]), 
   D[f, t] == S (1 + D[v, x] - g[f/S, ξ]) - f D[v, x]
 }; 

initc[ρinit_, finit_, χinit_, vinit_] := {
  ρ[0, x] == ρinit, f[0, x] == finit, χ[0, x] == χinit, v[0, x] == vinit
 };

boundc = {
   ρ[t, -1] == ρ[t, 1], f[t, -1] == f[t, 1], 
   χ[t, -1] == χ[t, 1], v[t, -1] == v[t, 1]
 };

sol = First@
   NDSolve[
    Join[
     eqns[ρ[t, x], v[t, x], f[t, x], χ[t, x]], 
     initc[1, β, 1 - 10^-4 Cos[x], 0], 
     boundc
    ], 
    {ρ, v, f, χ}, {x, -1, 1}, {t, 0, 0.05}, 
    Method -> {"IndexReduction" -> Automatic}
   ];

ρ = ρ /. sol[[1]];
v = v /. sol[[2]];
f = f /. sol[[3]];
χ = χ /. sol[[4]];