FreeFEM++ converting equation into code

Question

I am trying to solve the following problem nuemrically:

$$u_t = \Delta u + \sin t$$

To that effect I scanned the documentation of FreeFEM, the closest example I can find to my problem is the Thermal diffusion problem.

I am trying to modify the existing tutorial code to solve my problem instead of the one in the example.

My current biggest problem is that I don't know how to introduce the $\sin t$ part into the equation.

I tried this:

// Parameters
func u0 = 10. + 90.*x/6.;
real ue = 25.;
real alpha=0.25;
real T=5.;
real dt=0.1 ;

// Define mesh boundary
border C(t=-pi, pi){x=cos(t) -1/(1 + 4 * t * t); y=sin(t);}

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

// Fespace
fespace Vh(Th, P1);
Vh u=u0, v, uold, t;

// Problem
problem thermic(u, v)
    = int2d(Th)(
          u*v/dt
        + 1.0 *(
              dx(u) * dx(v)
            + dy(u) * dy(v)
        )
    )
    - int2d(Th)(
          uold*v/dt
    )
//     + int1d(Th)(
//           t
//     )
    + int1d(Th, 1, 3)(
          alpha*u*v
    )
    - int1d(Th, 1, 3)(
          alpha*ue*v
    )
    + on(3, 4, u=u0)
    ;

// Time iterations
ofstream ff("thermic.dat");
for(real t = 0; t < T; t += dt){
    uold = u; //equivalent to u^{n-1} = u^n
    thermic; 
    ff << u(3., 0.5) << endl;
    plot(u);
}

But if I comment out:

//     + int1d(Th)(
//           t
//     )

I cannot compile. I also tried adding a t to the problem declaration and that also leads to compilation errors. How are you supposed to introduced parts of the problem that depend on time?

For example if I try adding $t$ to the problem declaration I get this error:

-- FreeFem++ v4.12 (Tue Dec  6 19:14:11 CET 2022 - git v4.12)
   file : mycode.edp
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : // Parameters
    2 : func u0 = 10. + 90.*x/6.;
    3 : real ue = 25.;
    4 : real alpha=0.25;
    5 : real T=5.;
    6 : real dt=0.1 ;
    7 : 
    8 : // Define mesh boundary
    9 : border C(t=-pi, pi){x=cos(t) -1/(1 + 4 * t * t); y=sin(t);}
   10 : 
   11 : // The triangulated domain Th is on the left side of its boundary
   12 : mesh Th = buildmesh(C(100));
   13 : 
   14 : // Fespace
   15 : fespace Vh(Th, P1);
   16 : Vh u=u0, v, uold, t;
   17 : 
   18 : // Problem
   19 : problem pde(u, v, t)
 Error line number 19, in file mycode.edp, before  token )
 Error in test or unkwon function (odd number of function) 
  current line = 19
Compile error :  Error in test or unkwon function (odd number of function) 
        line number :19, )
error Compile error :  Error in test or unkwon function (odd number of function) 
        line number :19, )
 code = 1 mpirank: 0

All I get from this is that my syntax is wrong but it doesn't help me much to understand how to define the rhs/sin conditions.

Answer 1

In the piece of code that you mentioned from the FreeFEM documentation, the Galerkin Finite Element Method (FEM) is used for the spatial discretization and the Finite Difference Method (FDM) is used to discretize the time domain. Reffering to the Thermal diffusion problem in the FreeFEM documentation, the variational formulation of your problem is ($v$ is a trial function): $$ \int_\Omega d\Omega(\frac{u^n-u^{n-1}}{\delta t}v+\nabla u^n\nabla v)-\sin{t^n}\int_\Omega d\Omega v=0,\quad n=0,1,\dots\text{,} $$ where $t^n=n\delta t$ and $\delta t$ is the time step.

The following is the corrected version of your FreeFEM script:

// Parameters
func u0 = 0; // initial condition: u is zero everywhere
real T = 1.0; // final time
real dt = 0.1; // time step

// Define mesh boundary
border C(t=-pi, pi){x=cos(t) -1/(1 + 4 * t * t); y=sin(t);}

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

// Fespace
fespace Vh(Th, P1);
Vh u, v, uold;
real S;

// Problem
problem thermic(u,v) = int2d(Th)( u*v/dt + ( dx(u)*dx(v) + dy(u)*dy(v) ) )
    - int2d(Th)( uold*v/dt )
    - int2d(Th)( S*v ) // the source term: sin(t)*v
    + on(C, u=0); // Dirichlet B.C.: on Curve C the field variable u is zero.

// apply the initial condition:
u = u0;

// Time iterations
ofstream ff("thermic.dat");
for(real t = 0; t < T; t += dt){
    uold = u; //equivalent to u^{n-1} = u^n
    S = sin(t);
    thermic; 
    ff << u(3., 0.5) << endl; // write the value of u at x=3.0 and y=0.5 on the domain
    plot(u);
}

It is to be noted that this isn't the most efficient implementation of time-dependent problems (especially large problems) in the FreeFEM, according to its documentation.