Using the adaptmesh Function for a Time-Dependent Source Function

6 次查看（过去 30 天）

Safiya 2024-11-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2168323-using-the-adaptmesh-function-for-a-time-dependent-source-function

编辑： Torsten 2024-11-25

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 !

12 个评论
显示 10更早的评论隐藏 10更早的评论

Torsten 2024-11-25

编辑：Torsten 2024-11-25

在 MATLAB Online 中打开

@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 2024-11-25

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 2024-11-25

编辑：Torsten 2024-11-25

在 MATLAB Online 中打开

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)).

请先登录，再进行评论。

请先登录，再回答此问题。