gaz3ll3 / femex Goto Github PK
View Code? Open in Web Editor NEWuse mex solve 2D second order linear PDE and transport equation
License: MIT License
use mex solve 2D second order linear PDE and transport equation
License: MIT License
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.
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 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.
sometimes the equation might not be linear, it is advisable to make it run with Newton's method if it is stable.
If it is a vector, then column vector should be default.
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.
the libraries loading is not working for 2016a.
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:
replacing integral in PDEcs.
priority minor.
in solver.cc
Try to build library with it. support static and dynamic,
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.
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.
though the .so file has been replaced.
workaround: put ones(10) * ones(10); clear ans;
in startup.m
In PDEcs module,
there should be some functions like:
umfpack
AMG
as preconditioner
merge 3d fmm into mainstream
It is quite straightforward to consider using the FEM code as forward solver and some preconditioner method/ krylov method for optimization, like BFGS/LBFGS.
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.
when scattering is large, use diffusive equation to acclelerate.
previous treecode is slow and not accurate enough, needs to update.
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.
use lumped mass to replace actual mass matrix to accelerate.
make it vectorized.
The lumped mass is directly extracted from original matrix, that will involve some thin layer elements outside, though it will not harm the result.
FMM is expected to be an more efficient version of treecode.
quadtree is not going to be in femex anymore.
Need to set preload library in order to run with un-modified matlab.
In meshing
serial running.
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.
add 3d case
currently fmm are running serially. Now there are two options to do it.
The concurrent queue should be a good choice.
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.
Tried with Matlab 2014b
.
Problem should be from OpenMP used with Boundary module. Now it is reset to serial version.
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.
building matrix is slow.
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.
most cases will have ~10^-3 relative error. not matter how we choose the surfaces.
In order to usequiver
to show gradients of solution on boundary.
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.
with 3595552
separate pml and domain.
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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.