separate pml and domain.

Since our problem works on a geometric shape. When doing iteration for optimization, instead of doing BFGS or something else directly on the nodes/elements(fine nodes), probably we need a rough idea on doing optimization on several sub-domains. Then just like multi-level method, move our solution towards the optimal one.

Problem:

• Algorithm on how to divide nodes/elements(Aggregates)
• Integral difficulties, maybe only a few levels can be made. Decomp depends on element's shape, if our mesh is not structured, then we probably need a sorting on location of elements, and get a quad-tree like thing.
• After the last level, we leave the optimization with element-wise level(piecewise constant) solution, then throw it into node-based iteration would require 3x or more memory.
• According to above, the convergence rate is very important to us. Since coarse mesh gives a better rate but worse precision, slow convergence at the end would be the case for multilevel method.

the bbfmm part or kifmm use Eigen to provide most linear algebra operations, but the operations are limited to multiplications and divisions. We probably should use raw pointers instead.

MKL can run in serialized mode here, since all matrices are quite small. Not needed to run in parallel.

merge 3d fmm into mainstream

Since we only need double precision floating point operations, one way is to build the header with template initiated with double type. Then link all things with the library. Compile process will be great reduced.

The PML example does not have a source term added.

container unordered_set needs to be allocate memory in the process, which might cause problem.

a workaround is to split the all nodes into serveral buckets and each bucket will take a thread to work on.

Now the problem does not need to be solved, a 70+ seconds warming up is fine.

Need to set preload library in order to run with un-modified matlab.

separate mex and C++ codes into two parts.

Get Julia/Python interfaces for that. Since Matlab is not desirable sometimes, importing too many libraries.

Both can use cffi or similar wrapper or swig.

previous treecode is slow and not accurate enough, needs to update.

most cases will have ~10^-3 relative error. not matter how we choose the surfaces.

The Raylet structure defined as

typedef struct Raylet {
    int32_t elem;
    double first[2];
    double second[2];
} Raylet;

and sizeof gives 40 bytes(since there is padding __int64). first and second take 32 bytes.

divide and conquer is good, at last step, has to consider several cases, 9 cases.

x0/side ? x1/side
y0/side ? y1/side

the corner cases are not easy to handle since there are two intersections. In this case, we divide and conquer again. to reduce to the easier 5 cases.

The lumped mass is directly extracted from original matrix, that will involve some thin layer elements outside, though it will not harm the result.

• use OPENMP.
• optimize some parts. throw the intersecting boundary part away may help.

It is quite straightforward to consider using the FEM code as forward solver and some preconditioner method/ krylov method for optimization, like BFGS/LBFGS.

In meshing

Now using matrix to calculate the solution instead of iteration.

Set up: 75 s

Iteration takes around 5~6 s per iteration.
Matrix uses 6 s for building the matrix. And around 10 s to solve (or faster).

If matrix method is faster then iteration method will take less iterations.

Find a way to make the solving process faster, since umfpack gives accurate solution to 1e-13 or so, we only need accuracy at 1e-8. Some iteration methods might help. (AMG already tried, it requires sparse structure of the matrix, it will take more time and memory on amg.)

GMRES is the next to be tried.

Try to build library with it. support static and dynamic,

implement global regularization term, associate with a weight function.

roughly speaking, TV regularization term takes form in

\int \grad f \cdot \grad f dx

which is actually a stiffness matrix. however, it might only exists on boundary or part of the domain, depends on the problem. Now, we just take stiffness.

though the .so file has been replaced.

workaround: put ones(10) * ones(10); clear ans; in startup.m

Now can add OpenMP support to all other modules. However, meshing process is not necessay to be parallelized. And only some assembly process has to be faster.

make it vectorized.

sometimes the equation might not be linear, it is advisable to make it run with Newton's method if it is stable.

Since each nodes is stored twice in our model.

Now I think first intersection is enough.

typedef struct Raylet {
    int32_t elem;
    double first[2];
} Raylet;

which is 8 + 16 = 24 bytes(with padding). But have to set to guardian Raylet at last.

Now considering shift pyamg to C/C++ + mex combination, since Pyamg is totally serial, we could for sure avoid the annoying MPI problem with MATLAB.

building matrix is slow.

replacing integral in PDEcs.

priority minor.

• test performance of dot product with C++/C 's summing function. w/ OpenMP and w/ OpenMP.

with 3595552

in solver.cc

In PDEcs module,

there should be some functions like:

• Integral over boundary
  • Specifiy all boundary or part of boundary
• Integral over domain
  • Domain still cannot specify sub-domains
• Solver specification
  • Use umfpack
  • Use AMG as preconditioner
    • Usually, it will be slower if matrix size is not large, in Matlab, there is no way to run a huge matrix, anyway.
• Extract Neumann data
• Extract Dirichlet data (however, this is not necessary)
• Regularization terms
  • Total variation
  • Tikhonov
  • TSVD
  • Iterative
• Gradient for the regularization terms
• Objective function evaluation expenses estimate

when scattering is large, use diffusive equation to acclelerate.

use lumped mass to replace actual mass matrix to accelerate.

Need assorted solvers, now it uses MATLAB 's UMFPACK sparse solver, with reasonable performance.

ILUPACK's AMG solver is under integration into the project.

Considering using gmres or other Krylov methods.

Considering interfacing with Pardiso's solver, if MKL is present.

If it is a vector, then column vector should be default.

In order to usequiver to show gradients of solution on boundary.

For the solver, a refinement method is needed for meshing, along with upward projection operator and downward operator, i.e. mapping between grids.

Warning: this would take large memory for storing mesh mapping. Refining mesh now takes use of triangle's -r option, which reuse the previous nodes, but not the edges.

the libraries loading is not working for 2016a.

currently fmm are running serially. Now there are two options to do it.

• cilkplus. supported by gcc-5/g++-5 and icc/icpc, easy to program and spawn new task.
• tbb.
  • upward pass operates on source, bottom-up, on each level, the tasks are independent, which can run in parallel.
  • downward pass operates on tree, top-bottom, on each level, the tasks are also independent.

The concurrent queue should be a good choice.

• each task may add new tasks append on end
• when meeting current task, it should have all information needed from previous tasks.

add 3d case

since FEM does not work well. The good news is FDM will be much more easier to implement and easy to use. Again, the PML is used, otherwise, we need a larger domain instead of [0, 1] x [0, 1] domain to let the wave pass through without bounce. consider the bouncing path, we do not need it to be too large.

FMM is expected to be an more efficient version of treecode.

quadtree is not going to be in femex anymore.

serial running.

Tried with Matlab 2014b.

Problem should be from OpenMP used with Boundary module. Now it is reset to serial version.

Reference from stackoverflow