Matlab function 'parabolic' for 3d problem does not accept my source-term input as a vector
4 views (last 30 days)
Show older comments
I am re-using a program in Matlab I once wrote for solving a PDE in 2D. I used the 'parabolic' function which accepts a vector of length #elements (number of elements in the mesh) as input for a source-term. Now, as I want to use this program for a 3D calculation, this vector input does not work. I am using Matlab version R2018b, here is a minimal example for the 2D version which is running just fine
clear; clc; close all;
model = createpde;
%% geometry 2d
R1 = [3;4;0;1;1;0;0;0;1;1];
gdm = R1;
sf = 'R1';
ns = char('R1')';
g = decsg(gdm,sf,ns);
geometryFromEdges(model,g);
hmax = 0.02;
mesh = generateMesh(model,'Hmax',hmax,'GeometricOrder','linear');
% pdemesh(model,'ElementLabels','off')
%% constants etc.
dt = 10;
tlist = 0:dt:14400;
L = length(tlist);
%% Fmat - source term
p = mesh.Nodes;
t = mesh.Elements;
np=size(p,2);
nt=size(t,2);
Fmat = zeros(nt,L);
%% equation
u0 = 0;
c = 10;
d = 1000*1000;
a = 0;
g = 100;
initialTemp = 0;
setInitialConditions(model,initialTemp);
applyBoundaryCondition(model,'neumann','Edge',1:4,'q',g,'g',0); % 2d
ttlist = [tlist(1) tlist(2)]
adjResults = parabolic(u0,ttlist,model,c,a,Fmat(:,1)',d);
and this is the code in 3D which gives the error "Error using pde.EquationModel/assema (line 52) The number of columns in the f-coefficient matrix, 7836, is not equal one."
clear; clc; close all;
model = createpde;
%% geometry 3d
[x,y,z] = meshgrid(0:1:1);
% Create the convex hull.
x = x(:);
y = y(:);
z = z(:);
K = convhull(x,y,z);
nodes = [x';y';z'];
elements = K';
hmax = 0.1;
geometryFromMesh(model,nodes,elements);
mesh = generateMesh(model,'hmax',hmax,'geometricorder','linear');
% pdeplot3D(model)
%% constants etc.
dt = 10;
tlist = 0:dt:14400;
L = length(tlist);
%% Fmat - source term
p = mesh.Nodes;
t = mesh.Elements;
np=size(p,2);
nt=size(t,2);
Fmat = zeros(nt,L);
%% equation
u0 = 0;
c = 10;
d = 1000*1000;
a = 0;
g = 100;
initialTemp = 0;
setInitialConditions(model,initialTemp);
applyBoundaryCondition(model,'neumann','Face',1:6,'q',g,'g',0); % 3d
ttlist = [tlist(1) tlist(2)]
adjResults = parabolic(u0,ttlist,model,c,a,Fmat(:,1)',d);
I didn't want to rewrite everything since I just wanted to briefly check something, but I cannot figure out what is different for the 3D implementation and why this input for the source term seems to be incorrect.
(The background of my implementation is to implement an adjoint equation with point source, i.e. with a dirac delta function for a certain sample point in space at a certain time as source term. Therefore I used the Fmat vector with length #elements, which is zero everywhere except at the element index where the sample point lies. Here, I scaled the source value by the volume of the element.
I also tried different implementations by using e.g. the @circlef function or by using other Matlab functions such as 'solvepde', but it was also quite complicated to come up with a working formulation for the source term.)
Can anyone help me with this?
Thanks a lot!
0 Comments
Answers (1)
MYBLOG
on 30 Aug 2023
Hello,
To resolve this issue, you might need to adjust how you create the source term vector for the 3D case. Consider whether the 'parabolic' function in 3D expects a different format for the source term compared to the 2D case. This could involve reshaping or reformatting the source term vector Fmat(:,1)' to meet the requirements of the 'parabolic' function for 3D calculations.
clear; clc; close all;
model = createpde;
%% geometry 3D
[x, y, z] = meshgrid(0:1:1);
x = x(:);
y = y(:);
z = z(:);
K = convhull(x, y, z);
nodes = [x'; y'; z'];
elements = K';
hmax = 0.1;
geometryFromMesh(model, nodes, elements);
mesh = generateMesh(model, 'Hmax', hmax, 'GeometricOrder', 'linear');
%% constants etc.
dt = 10;
tlist = 0:dt:14400;
L = length(tlist);
%% Fmat - source term
% In the 3D case, Fmat should be a 3D matrix, where each column corresponds to
% the source term vector for a specific time step.
Fmat = zeros(size(mesh.Elements, 2), L);
% Set the source term for the first time step (for example)
sample_element_index = 42; % Adjust this index to your desired sample point
Fmat(sample_element_index, 1) = 1; % Adjust the value as needed
%% equation
u0 = 0;
c = 10;
d = 1000 * 1000;
a = 0;
g = 100;
initialTemp = 0;
setInitialConditions(model, initialTemp);
applyBoundaryCondition(model, 'neumann', 'Face', 1:6, 'q', g, 'g', 0); % 3D
ttlist = [tlist(1) tlist(2)];
adjResults = parabolic(u0, ttlist, model, c, a, Fmat(:, 1), d); % Use Fmat(:, 1) as the source term
13 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!