Define the nonlinearity f inside and outside the time loop gave different result

General Discussion
1

Hello everyone,
I am solving a simple heat equation of variable u with the source term f.
At first, function f only depends on u so I can define f outside the time loops and it gives a good result. But then I want f to depend on time t and I try to define it inside the time loops. I found that even though I haven’t changed the function f, the result I obtained was different from the previous case.
Does anyone have an idea why I got this difference and how I can avoid it?

Here comes my code:

real r=0.5, D = 1., nuE = 0.08, muE = 0.05 ,muF = 0.1,  b = 10 , K = 200, dt=0.01, Tf=5. ;
real r1 = 1, r2 = 5, c= 0.05, lbda = 500;
border C(z=0, 2*pi){x=10*cos(z); y=10*sin(z);}

// The triangulated domain Th is on the left side of its boundary
mesh Th = buildmesh(C(100));

fespace Vh(Th,P1);
Vh uh,vh,uh0,f;

//Function f outside the time loop
f = r*nuE*b*uh/(b*uh/K - muE - nuE)*(1-exp(-uh))- muF*uh;

func init = 200*(sqrt(x*x+y*y) >= 7);
macro Grad(u)[dx(u),dy(u)]//
problem chaleur(uh,vh) = int2d(Th)(uh*vh/dt + Grad(uh)'*Grad(vh)*D) - int2d(Th, optimize = 0)(uh0*vh/dt + f*vh);

// Time loop
real t = 0;
uh0 = init;
for (int m = 0; m <= Tf/dt; m++){
    // Update
    t = t+dt;
    // Function f inside the time loop
    //f = r*nuE*b*uh/(b*uh/K - muE - nuE)*(1-exp(-uh))- muF*uh;

    // Solve
    chaleur;
    uh0 = uh;

    // Plot
    if (!(m % 10)) {
        plot(uh, cmm="t="+(t-dt)+"[days]",  value = true, fill = true, wait=1);
    }
}

Thank you very much!

2

if you RHS depend on u the problem is no- linear
and f is not a function so the operator = do interpolation
if you put the f outside the loop f never change (it is f at time 0)
over wise in the loop it is f a time t-dt (previous time step).

it is better but you can also make a fixed point also to converge in f by adding a loop to solve the non linear problem

   for( int iter =0; iter < 10; ++i)
{
    // Function f inside the time loop
    f = r*nuE*b*uh/(b*uh/K - muE - nuE)*(1-exp(-uh))- muF*uh;
   Vh up;
    up[]=u[];: 
    // Solve
    chaleur;
up[] -= u[];
cout << " iter = " <<iter<< " err=" <<  up[].linfty << endl;
if( up[].linfty < 1e-5) break; 
}