2
$\begingroup$

In the process of discretization of a 4D PDE, I need to construct a final sparse matrix $B$, which is very large (I denoted by $\text{size}$ here) and time-dependent, viz, some of its entries change with $t$. To move into the 4D discretization matrices, I use the command KroneckerProduct[].

A sample case is as follows:

ClearAll["Global`*"];
SeedRandom[1];

m1 = 8; m2 = 8; m3 = 8; m4 = 8; size = m1*m2*m3*m4; 
gamma = 0.3; etad = 0.7/100; etaf = 1.2/100;
rosv = -0.4; rosd = -0.15; rovd = 0.3; rosf = -0.15; rovf = 0.3; rodf = 0.25; 
e = 100.; κ = 0.5; η = 0.1; lambdad = 0.01; lambdaf = 0.05;

c1 = 0.074; c2 = 0.014; c3 = 2.10; thetad[t] := c1 - c2*Exp[-c3*t]
f1 = 1.0; f2 = 0.5; f3 = 0.5; thetaf[t] := f1 - f2*Exp[-f3*t]

nx = Sort@RandomReal[{0, 100}, m1];
ny = Sort@RandomReal[{0, 10}, m2];
nz = Sort@RandomReal[{-1, 1}, m3];
nw = Sort@RandomReal[{-1, 1}, m4];
origrid = Flatten[Outer[List, nx, ny, nz, nw], 3];

Idx = SparseArray[{{i_, i_} -> 1.}, {m1, m1}, 0]; Idy = 
 SparseArray[{{i_, i_} -> 1.}, {m2, m2}, 0]; 
Idz = SparseArray[{{i_, i_} -> 1.}, {m3, m3}, 0]; Idw = 
 SparseArray[{{i_, i_} -> 1.}, {m4, m4}, 0];
DS = KroneckerProduct[(SparseArray@DiagonalMatrix@nx), Idy, Idz, Idw]; 
DV = KroneckerProduct[Idx, (SparseArray@DiagonalMatrix@ny), Idz, Idw]; 
DRD = KroneckerProduct[Idx, Idy, (SparseArray@DiagonalMatrix@nz), Idw];
DRF = KroneckerProduct[Idx, Idy, Idz, (SparseArray@DiagonalMatrix@nw)];
Id = KroneckerProduct[Idx, Idy, Idz, Idw];

dudx = SparseArray[{Band[{1, 1}] -> 1., Band[{2, 1}] -> 1.5, 
    Band[{1, 2}] -> -0.3}, {m1, m1}];
d2udx2 = SparseArray[{Band[{1, 1}] -> 2., Band[{2, 1}] -> 1., 
    Band[{1, 2}] -> 0.6, Band[{3, 1}] -> 1.}, {m1, m1}];
dudy = SparseArray[{Band[{1, 1}] -> 0.75, Band[{2, 1}] -> 5., 
    Band[{1, 2}] -> 2.}, {m2, m2}];
d2udy2 = SparseArray[{Band[{1, 1}] -> 3.5, Band[{2, 1}] -> 1., 
    Band[{1, 2}] -> 2., Band[{3, 1}] -> 0.1}, {m2, m2}];
dudz = SparseArray[{Band[{1, 1}] -> 1.2, Band[{2, 1}] -> 1., 
    Band[{1, 2}] -> 2.}, {m3, m3}];
d2udz2 = SparseArray[{Band[{1, 1}] -> 1.5, Band[{2, 1}] -> 1., 
    Band[{1, 2}] -> 1., Band[{3, 1}] -> 0.3}, {m3, m3}];
dudw = SparseArray[{Band[{1, 1}] -> 5., Band[{2, 1}] -> 2., 
    Band[{1, 2}] -> 1.}, {m4, m4}];
d2udw2 = SparseArray[{Band[{1, 1}] -> 3., Band[{2, 1}] -> 0.3, 
    Band[{1, 2}] -> 2., Band[{3, 1}] -> 0.2}, {m4, m4}];

B0 = SparseArray[
    +(1/2 (DS^2).(DV)).KroneckerProduct[d2udx2, Idy, Idz, Idw]
     + (1/2 gamma^2 DV).KroneckerProduct[Idx, d2udy2, Idz, Idw]
     + (1/2 etad^2)*KroneckerProduct[Idx, Idy, d2udz2, Idw]
     + (1/2 etaf^2)*KroneckerProduct[Idx, Idy, Idz, d2udw2]
     + ((rosv*gamma)*(DS.DV)).KroneckerProduct[dudx, dudy, Idz, Idw]
     + ((rosd*etad)*(DS.Sqrt[DV])).KroneckerProduct[dudx, Idy, dudz, 
       Idw]
     + ((rovd*(gamma*etad))*Sqrt[DV]).KroneckerProduct[Idx, dudy, 
       dudz, Idw]
     + ((rosf*etaf)*(DS.Sqrt[DV])).KroneckerProduct[dudx, Idy, Idz, 
       dudw]
     + ((rovf*gamma*etaf)*(Sqrt[DV])).KroneckerProduct[Idx, dudy, Idz,
        dudw]
     + (rodf*etad*etaf)*KroneckerProduct[Idx, Idy, dudz, dudw]
     + ((DRD - DRF).DS).KroneckerProduct[dudx, Idy, Idz, Idw]
     + (κ (η*Id - DV)).KroneckerProduct[Idx, dudy, Idz, 
       Idw]
     - DRD.Id
    ]; // AbsoluteTiming
matcoef1 = (SparseArray[{{i_, i_} -> c1 - c2*Exp[-c3*t]}, {size, 
     size}, 0.]);
matcoef2 = (SparseArray[{{i_, i_} -> f1 - f2*Exp[-f3*t]}, {size, 
     size}, 0.]);
B[t_] := B0 + (lambdad (matcoef1 - DRD)).KroneckerProduct[Idx, Idy, 
      dudz, Idw] + ((lambdaf (matcoef2 - DRF)) - 
       rosf*etaf*Sqrt[DV]).KroneckerProduct[Idx, Idy, Idz, 
      dudw]; // AbsoluteTiming

B = B[t]; // AbsoluteTiming
MatrixPlot[B]

I know that it is correct, but is not at all efficient. For example, for the above case, only constructing $B$ takes around 27 seconds in my laptop, while I should choose much more inputs for $m_1,m_2,m_3,m_4$ in my real calculations.

I think the most time-consuming case is where at the end I involve matcoef1 and matcoef2 into calculations. Please note that I need to save $B$ since I will impose the boundary conditions later (based on a combination of its rows).

Does anyone know how we can speed up the process of constructing sparse matrices involving time-dependent entries for large sizes? I also tried to use Module or Compile, but I failed!

$\endgroup$
  • $\begingroup$ It is absolutely not clear to me why it should not be sufficient to compute B[t] only for numeric t. Using symbolic matrices of this size just feels wrong. $\endgroup$ – Henrik Schumacher Nov 16 '18 at 21:18
  • $\begingroup$ Thanks for your care. You're right, it seems strange, but I don't have any other options. For example, consider that I should add a multiplier of rows 1000-1010 into the rows 2050-2060! I mean I impose the boundary conditions of the PDE. And then, I will get a final B[t] and give it to NDSolve. Because of this, I need to compute this fast. However, even if I do not impose the boundaries, NDSolve will need B[t] symbolically to do the time stepping. Am I right? If you know any tips, please write. $\endgroup$ – Fazlollah Nov 16 '18 at 21:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.