Matlab function 'parabolic' for 3d problem does not accept my source-term input as a vector

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!

Answer 1

[MYBLOG]

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

[superju]

Hm, but that's basically what I did as well? My Fmat = zeros(nt,L); was just an example, but has the same size than what you provided. Running this yields the same error for me: "Error using pde.EquationModel/assema (line 52) The number of rows in the f-coefficient matrix, 7836, is not equal the number of equations in the system, 1."

Also, I was not able to find information about what exactly are the requirements of the 'parabolic' function for 3D calculations and why or how they differ from the 2D case.

[Torsten]

What is your Fmat supposed the represent ? I only know source terms as a function of position and time.

[superju]

You can provide the input for the source term f either as a function (by a function handle f = @fun) or as a vector, where each entry of the vector defines the evaluation of the source term on an element of the finite element mesh. On this page it is said that f can be "specified as a scalar, matrix, character vector, character array, string scalar, string vector, or coefficient function. f represents the f coefficient in the scalar PDE", also that "(...) parabolic (...) produces the solution to the FEM formulation of the scalar PDE problem on a 2-D or 3-D region".

[Torsten]

And how do you address the element number of the grid where you want to specify the source term over time ? And how do you get its volume to divide the source by it ? I think by "matrix" it is meant that you define constant source terms for more than one PDE to be solved. In your case: How should MATLAB know that your vector Fmat(1,:) is the evolution of the source term over time in a specific element of your grid ?

[superju]

I have a defined sample position and starting from my triangulization, I can find out to which element this position belongs by e.g.

tn = pointLocation(TR, samplePosition'); 

and its volume by

[ar,~,~,~]=pdetrg(p,t);

The evolution of the source term is done by iterating through the Fmat matrix since every column there stands for the source term at a certain time step. As I said, in 2D this works well (not only the minimal working example above, but also my full implementation) and I get the correct results by doing so. I just don't understand why it won't work with this kind of input for my 3D problem, and the Matlab page on the 'parabolic' does not provide me enough information about the correct form of the input for f.

[Torsten]

Where so you supply the information about the location of the source ? Shouldn't this happen somewhere in the call to "parabolic" ? Is it hidden in "model" ?

[superju]

The location of the source in my full implementation is only encoded in the Fmat, where everything but the entry at the element (which contains the sample position) is 0. Does this answer your question?

Anyways, I still have the problem that I don't know which kind of input is accepted for f in the 3D case for the parabolic function.

[Torsten]

What if you provide c,a and d of the same size as f ?

And if you want to supply the source term at a certain time and not over the full integration period, shouldn't F be of size nT x nE (T: time, E: elements) or something similar ? Is your 2D model also time-dependent ?

[superju]

I did this, also for g, and I also provided u0 as a vector of length #points, but I still only get this error about wrong input for f.

Yes, my 2D model is also time-dependent and the whole matrix Fmat is of size #elements nt x #timesteps L. The evolution in my full implementation is performed in a loop over the number L of timesteps, but unfortunately, for the 3D case, I get this error already in the first iteration step (that's why I've shortened this to get the minimal working example).

[Torsten]

Could you supply a working code in 2d ?

[superju]

The full code is quite complex, but can be broken down to the 2D minimal working example I provided in my original post, just without this loop for the timesteps:

for j=1:(L-1)
    ttlist = [tlist(j) tlist(j+1)];
    adjResults = parabolic(u0,ttlist,model,c,a,Fmat(:,j)',d);
    adjResultsPDE = createPDEResults(model,adjResults,ttlist,'time-dependent');
    % etc. etc.
    adjSol(:,j+1) = adjResults(:,end);
end

[Torsten]

Well, if this really worked for the 2d-case, I'm out of ideas.

Maybe MATLAB support can tell about the differences between inputs for the 2d- and the 3d- case.

[superju]

Yeah, would be good to know about these differences, I couldn't find any other hint than the one on the Matlab page for the 'parabolic' function.

Maybe someone else here has a clue?