Parabolic problem in 3D

Hi,
i write the Following code to calculate the solution of two problems
du/dt-Delta u+ qdelta(x/eps,y/eps,z/eps) u/(1+u)= 0.25sin(2pi*x), (x,y,z) \in (0,1)^3, t \in [0,1)

avec la condition initiale u(x,0)=0 et les conditions au ord périodique.

qdelta(t,s,c)= 1 si (x,y,z) \in B_{delta}+ \Z^3 et 0 sinon.

le rayon de la boule B_{delta}= (3./4*(1./(pi*w)))^1/3.

et l’autre problème est

dw/dt-Delta w+w/(1+w)=0.25sin(2pi*x)

avec condition au bord periodique et condition initiale w(x,0)=0.

So i try with Following code:

Is my code correct? Please. Beacause when i visualize uh and wh there are static! is not normal. And if i try eps=0.5 and w=3, i see the ball in qdelta but the boundary of ball isn’t regular
Kind regards

 load "UMFPACK64"
defaulttoUMFPACK64();

load "iovtk"

load  "msh3"


verbosity=1;

func real qdelta (real t, real s, real c, real w) {
  real t1 = t - rint(t);
  real s1 = s - rint(s);
  real c1 = c -rint(c);
  real delta= ((3./4.)*(1./pi*w))^1/3;
  if ( square(s1)+square(t1)+square(c1) <= square(delta)) 
    return w;
  else  
    return 0.0;
  };


func real q (real t, real s, real c,real w) {
  real t1 = t - rint(t);
  real s1 = s - rint(s);
  real c1 = c -rint(c);
   real delta= ((3./4.)*(1./pi*w))^1/3;
  if ( square(s1)+square(t1)+square(c1) <= square(delta)) 
    return 1;
  else  
    return 0.0;
  };



func real f (real x)
  {
  return 0.25*sin(2*pi*x);
  };

// initial values for the spatial domain
int N=20;
real eps=0.1, w=2000;

real t=0.00;
real dt=0.1;


real[int] h;     // value of h 
real[int] error; // error in norm H1 
ofstream fout("error.txt");
fout << "#  h error" << endl;


 ofstream fouti("errH1.txt");
fouti << "# t errH1" << endl;

  mesh3 Th=cube(20,20,20);
  fespace Vh(Th,P1,periodic=[[1,x,z],[3,x,z],[2,y,z],[4,y,z],[5,x,y],[6,x,y]]);
  fespace Wh(Th,P2,periodic=[[1,x,z],[3,x,z],[2,y,z],[4,y,z],[5,x,y],[6,x,y]]);
  Vh uh=0.0;
  Vh wh=0.0;
  Vh vh;
  Vh oldU=0.0;
  Vh oldW=0.0;
  real T=1.0;// arbitraire
  int[int] Order= [1];


  Vh fh=f;

  Vh nith = 0.0; // for nonlinear iterations
  Vh g0,g1,g2,g3;
  Wh uplot;

  // Nonlinear part is calculated in previous iteration (nith) and becomes a part of the right hand side
  problem het(uh,vh,init=t) =
    int3d(Th)((dx(uh)*dx(vh)+dy(uh)*dy(vh)+dz(uh)*dz(vh))*dt)
    +int3d(Th)(uh*vh)
    -int3d(Th)(oldU*vh)
    +int3d(Th)(qdelta(x/eps,y/eps,z/eps,w)*g0*vh*dt) +int3d(Th)(qdelta(x/eps,y/eps,z/eps,w)*g1*uh*vh*dt)
    -int3d(Th)(fh*vh*dt)
    ;

  problem hom(wh,vh,init=t) =
    int3d(Th)((dx(wh)*dx(vh)+dy(wh)*dy(vh)+dz(wh)*dz(vh))*dt)
    +int3d(Th)(wh*vh)
    -int3d(Th)(oldW*vh)
    +int3d(Th)(g2*vh*dt) +int3d(Th)(g3*wh*vh*dt)
    -int3d(Th)(fh*vh*dt)
    ;
    
  uplot=q(x/eps,y/eps,z/eps,w); 
//plot(uplot,fill=true,value=true,ps="q".eps");
savevtk("q.vtu",Th,uplot,dataname="qdelta", order=Order);
 cout << "mean value of qdelta: " << int3d(Th)(uplot);
 

  for (t=0;t<=T;t+=dt) {
    cout << "Time = " << t << endl;
    fh=f(t);
    oldU=uh;
    oldW=wh;
    int k=0;        
    int noit = 10; // maximum no of nonlinear iterations
    real errnlhet = 0;  // error of nonlinear iteration of het
    real normuh = 0;
    real errnlhom =0; // error of nonlinear iteration of hom
    real normwh=0;
    
    // Nonlinear iterations for heterogeneous solution
    nith = uh; // start nonlinear iterations  with last solution for het
    g0=uh*uh/(1.0+uh)^2;
    g1=1.0/(1.0+uh)^2;
    cout << "    Solving nonlinear heterogeneous problem.\n";
    for(k=1;k<noit; ++k) {
      het;
      normuh = sqrt(int3d(Th)(square(uh))); // L^2 nor of the solution of het
      errnlhet = sqrt(int3d(Th)(square(uh - nith))); // L^2 norm of the error in nonlinear iteration of het
      cout << "        " << k << ": error = " << errnlhet << endl;
      if(errnlhet <= 1E-5*normuh) // stopping criteria. If relative error
      break;        // is less than 1E-5 stop nonlinear iterations of het
      nith = uh; // prepare new iteration
      g0=uh*uh/(1.0+uh)^2;
      g1=1.0/(1.0+uh)^2;
      }
    // Check convergence for het  
    if(k >= noit) {
      cout << "Nonlinear solver for het did not converged.\n";
      cout << "Error = " << errnlhet << endl;
      cout << "|| uh || = " << normuh << endl;
      exit(1); // exit the program with error code 1
      }

    nith= wh; // start nonlinear iteration with last solution for hom
    g2=wh*wh/(1.0+wh)^2;
    g3=1.0/(1.0+wh)^2;
    cout << "    Solving nonlinear homogeneous problem.\n";
    for(k=1;k<noit; ++k) {
      hom;
      normwh = sqrt(int3d(Th)(square(wh))); // L^2 nor of the solution of hom
      errnlhom = sqrt(int3d(Th)(square(wh - nith))); // L^2 norm of the error in nonlinear iteration of hom
      cout << "        " << k << ": error = " << errnlhom << endl;
      if (errnlhom <= 1E-5*normwh) // stopping criteria. If relative error
        break;        // is less than 1E-5 stop nonlinear iterations of het
      nith = wh; // prepare new iteration
      g2=wh*wh/(1.0+wh)^2;
      g3=1.0/(1.0+wh)^2;
      }
    // Check convergence for het  
    if(k >= noit) {
      cout << "Nonlinear solver for hom did not converged.\n";
      cout << "Error = " << errnlhom << endl;
      cout << "|| wh || = " << normwh << endl;
      exit(1); // exit the program with error code 1
      }
  real errH1 = sqrt(int3d(Th)(square(uh - wh)) + int3d(Th)(square(dx(uh)-dx(wh)) + square(dy(uh)-dy(wh))+ square(dz(uh)-dz(wh))));
  cout << " ==== errH1 ( t = " << t <<", errH1 = " << errH1 << ") = "<< errH1 <<endl;
  fout << t <<" "<< errH1 <<endl;
 }

savevtk("u.vtu",Th,uh,dataname="uh", order=Order);
savevtk("w.vtu",Th,wh,dataname="vh", order=Order);

   real errH1 = sqrt(int3d(Th)(square(uh - wh)) + int3d(Th)(square(dx(uh)-dx(wh)) + square(dy(uh)-dy(wh))+ square(dz(uh)-dz(wh))));
	
  cout << " ==== errH1 ( eps = " << eps <<", omega = " << w << ", eps^2*omega = " << eps^2*w << ") = "<< errH1 <<endl;
  

  h=eps^2*w;