Determine locations struct used in PDE solver

Hi,

I'm trying to solve a 2D PDE where the 'c' coefficient is a function of the (x,y) coordinate. In order to specify this coefficient, I'm trying to create a matrix that is of length N by M where N are the number of x points and M the number of y points used to solve the PDE. I can then pass this matrix to the PDE solver by calling a function that returns the value of the coefficient at (x,y).

Unfortunately I need to populate the matrix outside of the function call, since I am using the PDE solver to fit for some data so I need to be able to modify it in between function calls. Therefore, I need to know how many points I need to fit for.

I noticed I cannot access the locations struct before calling solvepde and that the number of (x,y) points is larger than the number of nodes in my mesh, so I can't form the matrix needed or find the value for the (x,y) pair for the C coefficient. Is it somehow possible to pre-compute the coordinates and number of points that will be used in the PDE solver?

Thanks in advance!

Answer 1

[Ravi Kumar]

You can write that code that you indented to use to "create a matrix that is of length N by M where ..." in the function itself. These points are Gauss point coordinates, you don't need to know them in advance, your funcion should use them to create the value of c at those points.

[Roberto Brenes]

Unfortunately I need to populate the matrix outside of the function call, since I am using the PDE solver to fit for some data so I need to be able to modify it in between function calls.

Therefore, I need to know how many points I need to fit for.

[Ravi Kumar]

Not sure if I understand your use case correctly, but if you want to modify some data between function calls and create c-cofficient based on the modified data, then it can be done by calling the modifying function from the c-coefficient function.

Can you give a simplified example, so that I can understand your difficulty.

[Roberto Brenes]

Sure, the following use case should hopefully explain what I'm trying to do. I omitted most of the PDE setup stuff apart from the part where I'm having an issue with. Basically, I'm just trying to fit for the 'c' coefficient through minimizing the error in the simulated PDE and the measured data through a least-squares minimization.

I can't give an initial guess if I don't know how big locations.x and locations.y are.

%Have tdata and ydata from experiment that correspond to the simulated
%points

%simulatePDE is detailed below
fun = @(cmat)simulatePDE(cmat,tdata,ydata);

%assume we have access to locations struct
c0 = repmat(5,numel(locations.x),numel(locations.y));

bestc = fminsearch(fun,c0);

function lsqsum = simulatePDE(cmat,tdata,ydata)
...
%Set-up PDE, PDE geometry, initial conditions and other required things for
%PDE
...
    

c = @(locations,~) cmat(x_idx,y_idx);
%x_idx is the index of the called location.x in the location.x
%matrix and similarly for y_idx
d = 1;
a = 1;
f = 0;

specifyCoefficients(pdem,'m',0,'d',d,'c',c,'a',a,'f',f);

result = solvepde(pdem,tdata);
u1 = result.NodalSolution;

lsqsum = sum((ydata - u1).^2);
end

[Ravi Kumar]

Thanks for the additional info.

The sample code indicates that you are solving a single PDE in 2-D domain. You are trying to find 'c' as a matrix in a pre-defined 2-D space (the geometric domain). That is, your final outcome would be a matrix of c at known x and y locations. If this is your goal, then the following approach might work for you.

For simplicity, lets say you are working with a unit square domain with original as one of its corner. Define in initial c value with, independent of PDE solvers call location. Say you have grid of 0.1, that is xg = 0:0.1:1, yg = 0:0.1:1, the define c0 as:

c0 = 5*ones(numel(xg));

The optimizer, fminsearch, calls simulatePDE with an updated cmat in each iteration, use this cmat and the KNOWN (xg,yg) to construct a gridded interpolant:

[xg,yg] = ndgrid(0:0.1:1);
cInterpolant = griddedInterpolant(xg,yg,cmat);

Then construct the c-coefficient function handle to interpolate and find the values of c at the points requested by solver (location.x,location.y):

c = @(locations,~) cInterpolant(location.x,location.y)

This should set up the problem the way you want.

Note 1: If your domain is not suitable to gridded data, then you can use the scatteredInterpolant.

Note 2: The way you have set the optimization problem might not be ideal, you have as many parameters as elements in cmat. I think a better option would be to define a function with undermined coefficient and optimize to determine the coefficient.

Regards,

Ravi

[Roberto Brenes]

Thanks so much for your reply! This makes sense and it should work for my problem. I wanted to determine whether 'c' varies spatially for my data, which is why I planned on having a coefficient for each location. I'll play around with the optimization to see whether this is feasible or to follow your tip.