Using the adaptmesh Function for a Time-Dependent Source Function

Hello !

I am trying to solve a time-dependent problem using the adaptmesh function in MATLAB, where the source term in Poisson's equation changes with time. Specifically, I want to update the source term at each time step and adapt the mesh dynamically to improve the solution’s accuracy at each stage.

The problem involves solving the Poisson equation with a time-dependent source function that has the following form:

where is a sinusoidal function that represents the moving source:

Here’s the basic structure of my code, but it’s not working as expected:

% Constants
R0 = 1e-15;  % Approximation of the distribution radius
L = 1e-13;   % Domain size (m)
v = 1e6;     % Proton speed (m/s)
tmax = 2*L/v; % Maximum time
dt = 1e-20;

A = 1e-14;         % Amplitude (m)
f = 1e9;           % Frequency (Hz)
phi = 0;           % Initial phase (radians)
x0 = @(t) A * sin(2 * pi * f * t + phi);

model = createpde;

% Geometry (rectangle)
R1 = [3,4,-L/2,L/2,L/2,-L/2,L/4,L/4,-L/3,-L/3]'; % Rectangle
g = decsg(R1);
geometryFromEdges(model,g);

figure; 
pdegplot(model,"EdgeLabels","on"); 
axis equal
title("Geometry With Edge Labels Displayed")

% Coefficients for Poisson's equation
c = 1;  % Diffusion coefficient
a = 0;  % Reaction term
sigma = 1e-15; % Standard deviation

% Time loop
t = 0;
while t <= tmax
    % Update the source term at each time step (f needs to be updated)
    f = @(location, state) exp(-(location.x - x0(t)).^2 / (2 * sigma^2)) / (sqrt(2 * pi) * sigma);

    % Specify coefficients for Poisson's equation
    specifyCoefficients(model,"m",0,"d",0,"c",c,"a",a,"f",f);

    % Boundary conditions
    applyBoundaryCondition(model,"dirichlet","Edge",(1:4),"u",0);

    % Adapt the mesh
    [u,p,e,t] = adaptmesh(g,model,c,a,f, 0, "tripick", ...
                          "circlepick", ...
                          "maxt",2000, ...
                          "par",1e-3);

    % Solve the PDE
    results = solvepde(model);

    % Extract the results
    u = results.NodalSolution;

    % Display the results
    figure;
    pdeplot(model,"XYData",u, "ZData",u);
    title(['Solution at time t = ', num2str(t)]);

    % Update time
    t = t + dt;
end

The issue is that the adaptmesh function does not accept my function f as it is written in the code. How can we dynamically adapt the mesh over time, specifically near a given location (such as the proton's position)?

Thank you in adavance for your help !

[Torsten]

Why do you set the first and second order time derivatives to 0

specifyCoefficients(model,"m",0,"d",0,"c",c,"a",a,"f",f);

, but specify a time-dependent source term

f = @(location, state) exp(-(location.x - x0(t)).^2 / (2 * sigma^2)) / (sqrt(2 * pi) * sigma);

?

What is the problem you are trying to solve ? Do you want to get u == 0 as a solution for all values for t ?

[Safiya]

Here is the initial equation I want to solve :

And I used a gaussian function to approximate the Dirac function.

[Torsten]

You can try to solve this equation for fixed values of t, but I'm quite sure that there should be a term in the equation that "connects" these stationary solutions over time.

Look up "Fundamental Solution" under

https://en.wikipedia.org/wiki/Laplace%27s_equation

Here, the stationary solution for a fixed value of t is explicitly given.

[Safiya]

The main issue that I have is with the adaptmesh function that doesn't take the source function f even for fixed value of t. I am not sure if I configured correclty the adaptmesh function or if there is a specific syntax for the f function.

[Torsten]

The problem as given is simple to solve and gives u == 0 as solution in the complete domain for all values of t.

If you want to solve a different problem (where you might need an adaptive grid or something similar), change your question accordingly.

[Paul]

@Torsten

What is u? I don't see that variable in the equation.

[Torsten]

@Paul

@Safiya tries to solve the problem

-(d^2u/dx^2 + d^2u/dy^2) = dirac([x y]-[x0 y0])

by approximating

dirac([x y]-[x0 y0])

by

% Update the source term at each time step (f needs to be updated)
f = @(location, state) exp(-(location.x - x0(t)).^2 / (2 * sigma^2)) / (sqrt(2 * pi) * sigma);

(for whatever reason, the y-dimension is ignored)

and setting boundary conditions as u == 0 on the edges of a rectangle.

I mistakenly wrote u == 0 is the solution, but it should be

u(x,y) = -log(sqrt((x-x0)^2+(y-y0)^2))/(2*pi)

according to

https://en.wikipedia.org/wiki/Laplace%27s_equation

(at least if the problem is defined without boundary condition).

[Paul]

Just to make sure I'm following, the phi in the equation shown in this comment is u in the solution. But now I'm confused about

d^2u/dx^2 + d^2u/dy^2 = dirac([x y]-[x0 y0])

The left hand side is a scalar, but the right hand side is vector?

[Torsten]

Mathematically, dirac(x) = 0 for x ~= 0 and dirac(x) = Inf for x = 0 where x can be any vector in IR^n.

Maybe "dirac" in MATLAB is differently defined.

[Torsten]

@Safiya

You used this example from the UserGuide to start with ?

c = -1;
a = 0;
f = @circlef;
numberOfPDE = 1;
model = createpde(numberOfPDE);
%Create a geometry and include it in the model.
g = @circleg;
%L = 1e-13;
%R1 = [3,4,-L/2,L/2,L/2,-L/2,L/4,L/4,-L/3,-L/3]'; % Rectangle
%g = decsg(R1);
geometryFromEdges(model,g);
%Plot the geometry and display the edge labels.
figure;
pdegplot(model,"EdgeLabels","on");
axis equal
title("Geometry With Edge Labels Displayed")

applyBoundaryCondition(model,"dirichlet","Edge",(1:4),"u",0);
[u,p,e,t] = adaptmesh(g,model,c,a,f,"tripick", ...
                      "circlepick", ...
                      "maxt",2000, ...
                      "par",1e-3);

Number of triangles: 258
Number of triangles: 515
Number of triangles: 747
Number of triangles: 1003
Number of triangles: 1243
Number of triangles: 1481
Number of triangles: 1705
Number of triangles: 1943
Number of triangles: 2155

Maximum number of triangles obtained.

figure;
pdemesh(p,e,t);
axis equal

%Plot the error values.
x = p(1,:)';
y = p(2,:)';
r = sqrt(x.^2+y.^2);
uu = log(r)/2/pi;
figure;
pdeplot(p,e,t,"XYData",u-uu,"ZData",u-uu,"Mesh","off");

figure;
pdeplot(p,e,t,"XYData",u,"ZData",u,"Mesh","off");

[Safiya]

Yes, that's right. I tried to adapt this code to my case where the refinement position of the mesh is updated with the position of the proton moving in time ( from ).

[Torsten]

For your code so far, you can just plot the result as

pdeplot([p(1,:)+x(time);p(2,:)+y(time)],e,t,"XYData",u-uu,"ZData",u-uu,"Mesh","off");

where [x(time),y(time)] is the position of the proton at t = time.

This works since the results for different values of t don't depend on one another and because the solution for (0,0) is just translated by the vector (x(time),y(time)).