Recent Questions - Computational Science Stack Exchange most recent 30 from scicomp.stackexchange.com 2020-01-18T20:32:05Z https://scicomp.stackexchange.com/feeds https://creativecommons.org/licenses/by-sa/4.0/rdf https://scicomp.stackexchange.com/q/34240 0 Why Householder transformation can not be chosen to be an identity matrix? sunshine https://scicomp.stackexchange.com/users/32679 2020-01-18T10:03:30Z 2020-01-18T12:57:55Z <p>For Householder transformation, we know that <span class="math-container">$H = I-uu^T$</span>, where <span class="math-container">$\|u\|_2=\sqrt{2}$</span>. When it acts on any vector <span class="math-container">$x$</span>, <span class="math-container">$Hx$</span> and <span class="math-container">$x$</span> is symmetric with respect to <span class="math-container">$span(u)^T$</span>. But I have read a monography "Stewart. Matrix Algorithm I: Basic Decomposition, 1998, SIAM". It is written as follows on Page 257:</p> <blockquote> <p>Combining these two observations we get Algorithm 1.1 — a program to generate a Householder transformation. Note that when <span class="math-container">$x = 0$</span>, any <span class="math-container">$u$</span> will do. In this case the program housegen returns <span class="math-container">$u = 2e_1$</span>, where <span class="math-container">$e_1$</span> is the first column of an identity matrix. This choice does not make <span class="math-container">$H$</span> the identity, but the identity with its <span class="math-container">$(1,1)$</span>-element changed to <span class="math-container">$-1$</span>. (In fact, it is easy to see that a Householder transformation can never be the identity matrix, since it transforms <span class="math-container">$u$</span>into <span class="math-container">$-u$</span>.)</p> </blockquote> <p>My question is that I donot understand why the Householder matrixcannot be an identity matrix. Because when <span class="math-container">$x=0$</span>, it means that we donot need to do any transformation, i.e., we can take <span class="math-container">$H = I$</span> identity matrix. But the author said we cannot take identity matrix. Where do I misunderstand? Thanks very much.</p> https://scicomp.stackexchange.com/q/34237 1 Initial condition for Kuramoto-Sivashinsky Dave https://scicomp.stackexchange.com/users/33890 2020-01-17T20:52:51Z 2020-01-17T22:39:37Z <p>For a project in my advanced numerical method class I have to solve the 1D <a href="https://en.wikipedia.org/wiki/Kuramoto%E2%80%93Sivashinsky_equation" rel="nofollow noreferrer">Kuramoto-Sivashinsky</a> equation of which I know little. I just know that it was derived the equation to model the diffusive instabilities in a laminar flame front. It reads </p> <p><span class="math-container">$$u_t + u u_x + \lambda u_{xx} + \eta u_{xxxx} = 0.$$</span></p> <p>Obviously I need an initial condition <span class="math-container">$u_0$</span> but I don't know what to choose. I would like to simulate something with a physical meaning.</p> <p>I'm also looking for advice and documentation on this equation.</p> <p>Thank you in advance.</p> https://scicomp.stackexchange.com/q/34235 3 Solving saddle point problem having non-invertible top-left block with a PETSc nested matrix janou195 https://scicomp.stackexchange.com/users/11916 2020-01-17T18:06:30Z 2020-01-17T18:06:30Z <p>My system is a symmetric FE problem with lagrange multipliers:</p> <p><span class="math-container">$Z=\begin{pmatrix}A &amp; C^T \\ C &amp; 0\end{pmatrix}$</span></p> <p>The matrix <span class="math-container">$A$</span> is positive semi-definite, non-invertible. The whole matrix is invertible.</p> <p>I am working on the development of a finite element software where I want to solve such problems with PETSc in parallel. At this time I would like the <span class="math-container">$Z$</span> matrix to be a <code>MATNEST</code> in which the <span class="math-container">$A$</span> and <span class="math-container">$C$</span> matrices are assembled independently.</p> <p>With the representative simple code at the end of this message (1D laplacian with Lagrange multipliers), and using a <code>PCFIELDSPLIT</code> preconditioner, the solver diverges. If I artificially force the matrix <span class="math-container">$A$</span> to be invertible (adding the command-line option <code>-force_invertible</code>), then I get the right solution.</p> <p>So, what can I do to solve such a system with a nested matrix?</p> <p>Another question: is it possible to use direct methods (Mumps) with nested matrices in parallel?</p> <pre><code>#include "petscsys.h" /* framework routines */ #include "petscvec.h" /* vectors */ #include "petscmat.h" /* matrices */ #include "petscksp.h" #include &lt;vector&gt; #include &lt;string&gt; #include &lt;iostream&gt; #include &lt;numeric&gt; // Try to solve with Petsc a simple 1D laplacian problem (or a series of springs) // o-////-o-////-o-////-o-////-o-////-o-////-o-////-o-////-o ... // // The boundary conditions (Dirichlet) are imposed using Lagrange multipliers. // The system to be solved is of the form: // // Z x = y // --^-- -^- -^- // [A Ct] [u] = [b] // [C 0 ] [l] [d] // // A is positive semi-definite (1 eigenvalue is 0) // Z is indefinite invertible static char help[] = "Saddle point problem: scalar 1D laplacian with lagrange multipliers.\n\n"; int main(int argc,char **argv) { // Initialization ------------------------------------------------ // MPI_Init(NULL, NULL); PetscErrorCode ierr; ierr = PetscInitialize(&amp;argc,&amp;argv,NULL,help);CHKERRQ(ierr); int rank, nproc; MPI_Comm_rank(PETSC_COMM_WORLD, &amp;rank); MPI_Comm_size(PETSC_COMM_WORLD, &amp;nproc); PetscViewer viewer = PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD); PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_DENSE); PetscBool forceInvertible=PETSC_FALSE; PetscOptionsGetBool(NULL,NULL, "-force_invertible", &amp;forceInvertible, NULL); // Problem definition and dof spliting among processes ----------- // double k=1.; std::vector&lt;int&gt; globalDofs={0, 1, 2, 3, 4, 5, 6, 7, 8}; std::vector&lt;std::pair&lt;int, double&gt;&gt; globalDirichlet={{0, 0.}, {8, 1.}}; // first: global dof number, second: imposed value auto start=[&amp;](int rank){ int q=(globalDofs.size()-1)/nproc; int r=(globalDofs.size()-1)%nproc; return q*rank + std::min(rank, r); }; std::vector&lt;int&gt; dofs; for (int i=start(rank); i&lt;=start(rank+1); ++i) { dofs.push_back(globalDofs[i]); } std::vector&lt;int&gt; isLocal(dofs.size(), 1); if (rank!=0) isLocal=0; int nLocalDofs=std::accumulate(isLocal.begin(), isLocal.end(), 0); std::vector&lt;int&gt; dirichletIdx; for (int j=0; j&lt;globalDirichlet.size(); ++j) { for (int i=0; i&lt;dofs.size(); ++i) { if (isLocal[i] &amp;&amp; globalDirichlet[j].first==dofs[i]) { dirichletIdx.push_back(j); } } } for (int p=0; p&lt;nproc; ++p) { MPI_Barrier(PETSC_COMM_WORLD); if (rank==p) { std::cout &lt;&lt; "Rank: " &lt;&lt; rank &lt;&lt; std::endl; std::cout &lt;&lt; " dofs : "; for (int d: dofs) { std::cout &lt;&lt; d &lt;&lt; " "; } std::cout &lt;&lt; std::endl; std::cout &lt;&lt; " isLocal: "; for (int d: isLocal) { std::cout &lt;&lt; d &lt;&lt; " "; } std::cout &lt;&lt; std::endl; std::cout &lt;&lt; " nLocalDofs: " &lt;&lt; nLocalDofs; std::cout &lt;&lt; std::endl; std::cout &lt;&lt; " dirichlet: "; for (int i: dirichletIdx) { std::cout &lt;&lt; "{" &lt;&lt; globalDirichlet[i].first &lt;&lt; ", " &lt;&lt; globalDirichlet[i].second &lt;&lt; "} "; } std::cout &lt;&lt; std::endl; } } // Matrix A ------------------------------------------------------ // Mat A; MatCreate(PETSC_COMM_WORLD, &amp;A); MatSetType(A, MATMPIAIJ); MatSetSizes(A, nLocalDofs, nLocalDofs, PETSC_DETERMINE, PETSC_DETERMINE); MatMPIAIJSetPreallocation(A, 5, NULL, 1, NULL); auto setValue=[&amp;](int i, int j, double v){ if (forceInvertible &amp;&amp; (i==2 || j==2)) { MatSetValue(A, 2, 2, 1., ADD_VALUES); } else { MatSetValue(A, i, j, v, ADD_VALUES); } }; for (int i=0; i&lt;dofs.size()-1; ++i) { // k * [ -1 1 ] // [ 1 -1 ] setValue(dofs[i] , dofs[i] , -k); setValue(dofs[i+1], dofs[i+1], -k); setValue(dofs[i] , dofs[i+1], k); setValue(dofs[i+1], dofs[i] , k); } MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); MatView(A, viewer); // Matrix C ------------------------------------------------------ // Mat C; MatCreate(PETSC_COMM_WORLD, &amp;C); MatSetType(C, MATMPIAIJ); MatSetSizes(C, dirichletIdx.size(), nLocalDofs, PETSC_DETERMINE, PETSC_DETERMINE); MatMPIAIJSetPreallocation(C, 1, NULL, 0, NULL); for (int i: dirichletIdx) { MatSetValue(C, i, globalDirichlet[i].first, 1., ADD_VALUES); } MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY); // Zero Matrix --------------------------------------------------- // Mat N; MatCreate(PETSC_COMM_WORLD, &amp;N); MatSetType(N, MATMPIAIJ); MatSetSizes(N, dirichletIdx.size(), dirichletIdx.size(), PETSC_DETERMINE, PETSC_DETERMINE); MatMPIAIJSetPreallocation(N, 0, NULL, 0, NULL); MatAssemblyBegin(N, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(N, MAT_FINAL_ASSEMBLY); // Matrix Z ------------------------------------------------------ // Mat Z; Mat subZ; subZ=A; MatCreateTranspose(C, &amp;subZ); subZ=C; subZ=N; MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, subZ, &amp;Z); MatView(Z, viewer); // Vector b ------------------------------------------------------ // Vec b; VecCreate(PETSC_COMM_WORLD, &amp;b); VecSetType(b, VECMPI); VecSetSizes(b, nLocalDofs, PETSC_DECIDE); VecAssemblyBegin(b); VecAssemblyEnd(b); // Vector d ------------------------------------------------------ // Vec d; VecCreate(PETSC_COMM_WORLD, &amp;d); VecSetType(d, VECMPI); VecSetSizes(d, dirichletIdx.size(), PETSC_DECIDE); for (int i: dirichletIdx) { VecSetValue(d, i, globalDirichlet[i].second, ADD_VALUES); } VecAssemblyBegin(d); VecAssemblyEnd(d); // Vector y ------------------------------------------------------ // Vec y; Vec suby; suby=b; suby=d; VecCreateNest(PETSC_COMM_WORLD, 2, NULL, suby, &amp;y); VecView(y, viewer); // KSP ----------------------------------------------------------- // KSP ksp; KSPCreate(PETSC_COMM_WORLD,&amp;ksp); KSPSetOperators(ksp, Z, Z); KSPSetFromOptions(ksp); PC pc; KSPGetPC(ksp,&amp;pc); PCSetType(pc, PCFIELDSPLIT); PCFieldSplitSetDetectSaddlePoint(pc, PETSC_TRUE); IS isg; MatNestGetISs(Z, isg, NULL); PCFieldSplitSetIS(pc, "u", isg); PCFieldSplitSetIS(pc, "l", isg); // PC pc; // KSPSetType(ksp, KSPPREONLY); // KSPGetPC(ksp,&amp;pc); // PCSetType(pc, PCLU); // PCFactorSetMatSolverType(pc,MATSOLVERMUMPS); Vec x; ierr = VecDuplicate(y, &amp;x); KSPSolve(ksp, y, x); // View ---------------------------------------------------------- // MatView(Z, viewer); VecView(y, viewer); KSPView(ksp, viewer); VecView(x, viewer); // --------------------------------------------------------------- // ierr = PetscFinalize(); MPI_Finalize(); return ierr; } </code></pre> https://scicomp.stackexchange.com/q/34234 3 Hit-n-Run Monte Carlo on convex polytope Davide Papapicco https://scicomp.stackexchange.com/users/33881 2020-01-17T12:42:41Z 2020-01-17T12:42:41Z <p>So, I'm currently trying to implement a MCMC to uniformly sampling hyper-points from the polytope defined as <span class="math-container">$\mathbb{K}=\{x\in\mathbb{R}^{n}\;\;\text{s.t.}\;\; A\,x=b \}$</span> in the specific case where, given a generic linear transformation <span class="math-container">$A\in\mathbb{R}^{m\times n}$</span>, <span class="math-container">$b\equiv0\in\mathbb{R}^m$</span> and the boundary conditions <span class="math-container">$0\leq x\leq1$</span> hold.</p> <p>Now, although I was able to succesfully perform the simulation (I am using Julia), there are some things I am not quite sure about:</p> <ul> <li><p>given the fact that polytopes of this this kind tend to have star-like shapes in higher dimensions, two pre-processing steps are required, before launching the simulation: </p> <ol> <li><p>the first regards the so-called <em>blocked-flux adjustment</em> which consists, quote from the text of the exercise, find fluxes <span class="math-container">$i$</span> such that <span class="math-container">$\max_{x\in K} x_i = \min_{x\in K} x_i = z_i$</span> and remove such variables from the system, adjusting the vector <span class="math-container">$b$</span>. Can anyone please explain to me what the heck this means?</p></li> <li><p>the second consists in finding an optimal inner point in <span class="math-container">$\mathbb{K}$</span> as a starting point of the chain, which intuitevely enough has to be located far from the vertexes of the polytope. The text tells me this: <em>it can be done e.g. by computing</em> <span class="math-container">$\frac{1}{2n}\sum_{i=1}^n (x^{\min,i} + x^{\max,i})$</span> <em>where</em> <span class="math-container">$x^{\min,i} \in \arg\min_{x\in K} x_i$</span> and <span class="math-container">$x^{\max,i} \in \arg\max_{x\in K} x_i$</span>. Here I just do not understand the notation: I suppose I should compute a weighted average of the midpoints of each edge of the polytope but I can't see how this is related to the above formulation.</p></li> </ol></li> <li><p>as any coherent MCMC, the walk in the state space has to be s.t. it satisfies the <strong>detailed balance</strong>, which, for the present case, the text tells me that it should be the target distribution <span class="math-container">$$p(x) \propto \delta^m(Sx-b)\prod_i \theta(u_i-x_i)\theta(x_i-l_i)$$</span> again, I have no idea how this is obtained nor how to compute them.</p></li> </ul> https://scicomp.stackexchange.com/q/34233 1 Efficient ways to numerically evaluate matrix exponentials ScientificPythonNovice https://scicomp.stackexchange.com/users/33823 2020-01-17T12:40:37Z 2020-01-17T19:42:55Z <p>What are some computationally efficient ways to solve matrix exponentials, i.e. functions of the form : f(X)=<span class="math-container">$e^{X}$</span>, where X is a square matrix ?</p> <p>So far I have been able to diagonalise some matrices, and find the exponent of individual diagonal elements, but not all matrices I'm dealing with will be diagonalisable.</p> <p>I am using Python with SciPy/NumPy, so solutions that can be implemented here will be most useful. If not, general solutions/solutions from other platforms are welcome too.</p> <p>EDIT: (1) I need the exponential itself, not a solution using it. (2) The matrix X is dense, typically small (3x3 or 4x4), need not be symmetric or Hermitian.</p> https://scicomp.stackexchange.com/q/34232 0 HSS preconditioner with gmres [closed] mh_zh https://scicomp.stackexchange.com/users/33879 2020-01-17T12:23:01Z 2020-01-18T02:38:16Z <p>I have a question about HSS preconditioner with GMRES method. For implementing the HSS preconditioner with GMRES, we need to solve the linear system of the form (I + H)(I + S)z =r, for a given r at each iteration. I have written MATLAB code for this algorithm, but it has errors. I have not found the correct way to modify it. How to solve the error in this code?</p> <p>This is my code:</p> <pre><code>clc; clear; close all; restart = 30; maxit = 500; tol = 1e-6; p =64; h=1/(p+1); m=p^2; n=2*p^2; N=m+n; I=speye(p); F=(1/h)*gallery('tridiag',p,-1,1,0); T=(1/h^2)*gallery('tridiag',p,-1,2,-1); B=[kron(I,F),kron(F,I)]'; A=[kron(I,T)+kron(T,I) sparse(m,m);sparse(m,m) kron(I,T)+kron(T,I)]; AA=[A B;-B' sparse(m,m)]; u0=zeros(N,1); f=ones(n,1); g=ones(m,1); rhs=[f;g]; x0=zeros(n,1); y0=zeros(m,1); alpha=10^-4; %% Hss tic; M = @(x)hsspre_inexact(x,alpha,A,B,rhs,x0,y0,maxit); [x,flag,relres,iter,resvec]=gmres(AA,rhs,restart,tol,maxit,M); time=toc; iter = (iter(1)-1)*restart+iter(2); % solve the sub-linear system function p = hsspre_inexact(alpha,A,B,rhs,x0,y0,maxit) % P_hss = [alpha*In+A O ] * [alpha*In B'] % [ O alpha*Im] [-B alpha*Im] [m,n]=size(B); In = speye(n); Im = speye(m); for i = 1:maxit p = [x0;y0]; f = rhs(1:n); g = rhs((n+1):(n+m)); r1=f-A*x0-B'*y0; r2=g+B*y0; u1 = pcg(alpha*In+A,r1); u2 = (1/alpha)*(r2+B*u1); z2 = pcg((1/alpha)*B*B'+alpha*Im,r2); z1 = (1/alpha)*(u1-B'*u2); p = [z1;z2]; end end </code></pre> <p>This is an error:</p> <pre><code>Error in gmres (line 228) r = iterapp('mldivide',m1fun,m1type,m1fcnstr,r,varargin{:}); </code></pre> https://scicomp.stackexchange.com/q/34230 1 Using adolc for the sign function in c++ Smilia https://scicomp.stackexchange.com/users/29701 2020-01-17T10:10:55Z 2020-01-18T06:27:55Z <p>Here is an implementation of the sign function in C++ using Adolc librairy for automatic differentiation.</p> <pre><code>template&lt;class Tdouble&gt; Tdouble sgn(const Tdouble &amp;x) { Tdouble s_plus, s_minus, half(.5); // set s_plus to sign(x)/2, except for case x == 0, s_plus = -.5 condassign(s_plus, +x, -half, +half); // set s_minus to -sign(x)/2, except for case x == 0,s_minus = -.5 condassign(s_minus, -x, -half, +half); // set s to sign(x) return 0.5*(1-(s_plus - s_minus)); } </code></pre> <p>My question is: why do we need to compute s_minus and s_plus ? what are the advantages ?</p> <p>What if I use simply :</p> <pre><code>template&lt;class Tdouble&gt; Tdouble sgn(const Tdouble &amp;x) { Tdouble res; condassign(res,x,1,-1); return res; } </code></pre> https://scicomp.stackexchange.com/q/34226 1 Determinant of a matrix after removing or adding lines and columns Christophe https://scicomp.stackexchange.com/users/360 2020-01-17T06:23:15Z 2020-01-17T07:57:57Z <p>In quantum mechanics, the wavefunction of N electrons is given by a determinant. I am working on a Monte Carlo algorithm. At each Monte Carlo step, I need to add or remove an electron, which means adding or removing a line and a column to my matrix.</p> <p>I use Gauss elimination to compute an upper triangle matrix. Adding a line and a column is simple and can be done in a time growing as <span class="math-container">$O(n^2)$</span>. Removing the last line and column is also feasible simply in <span class="math-container">$O(n^2)$</span>. To the ith line, I saw no other possibilty than recomputing the pivots from i to n-1 which requires a time <span class="math-container">$O(n^3)$</span>.</p> <p>Does anybody see another possibily?</p> https://scicomp.stackexchange.com/q/34223 2 Inverting really big symmetric block matrix Caspar https://scicomp.stackexchange.com/users/33871 2020-01-17T00:29:00Z 2020-01-17T00:29:00Z <p>I have a really big symmetric 7.000.000 X 7.000.000 matrix that i would like to invert. The matrix is extremely sparse and it can be rearranged as to become a block matrix. The biggest blocks are around 1500 X 1500 but most blocks are much smaller like 4 X 4. What is an efficient way to find the inverse? The scipy.sparse.linalg.inv in python gives me a memory error.</p> https://scicomp.stackexchange.com/q/34222 4 Why the two Gram-Schmidt algorithms produce different results for qr factorization? sunshine https://scicomp.stackexchange.com/users/32679 2020-01-16T23:35:18Z 2020-01-16T23:40:54Z <p>For the qr factorization using classic Gram-Schmidt algorithm, I found the 2 different implementations below. The first one uses the <code>for</code> loop to compute the upper triangular matrix <code>R</code>, but the second one uses the matrix-vector multiplication. Since they are mathematically equivalent, but they produce different results. I cannot find why it is so, can anyone give me some indications? Furthermore, I often find that whenever I understand an algorithm theoretically, if I implemented it, I will have some trouble because of different numerical results. This is why I always ask questions about how, actually, to implement the algorithm with specific code. Because even if I know the algorithm, I still cannot generate the correct results. Thanks.</p> <pre class="lang-matlab prettyprint-override"><code>% test the classical Gram-Schmidt algorithm clc;clear;format compact; A = hilb(7); %% method 1 R(1,1) = norm(A(:,1)); Q(:,1) = A(:,1)/R(1,1); m = size(A,2); for j=2:m for i=1:j-1 R(i,j) = Q(:,i)'*A(:,j);% the difference between 2 methods from below end q_hat = A(:,j)-Q(:,1:j-1)*R(1:j-1,j); R(j,j) = norm(q_hat); Q(:,j) = q_hat/R(j,j); end %% method 2 R1(1,1) = norm(A(:,1)); Q1(:,1)=A(:,1)/R1(1,1); for j=2:m R1(1:j-1,j) = Q1(:,1:j-1)'*A(:,j); % the difference between 2 methods from above q_hat = A(:,j)-Q1(:,1:j-1)*R1(1:j-1,j); R1(j,j) = norm(q_hat); Q1(:,j) = q_hat/R(j,j); end %% compare norm(A-Q*R) norm(Q'*Q-eye(m)) norm(A-Q1*R1) norm(Q1'*Q1-eye(m)) </code></pre> <p>My results are as follows:</p> <pre class="lang-matlab prettyprint-override"><code>ans = 0 ans = 0.3369 ans = 8.9238e-10 ans = 0.1890 </code></pre> https://scicomp.stackexchange.com/q/34221 2 pdepe or Crank-Nicolson? How much is pdepe good? Na'omi https://scicomp.stackexchange.com/users/32414 2020-01-16T22:34:48Z 2020-01-17T00:46:03Z <p>I am beginner in MATLAB and similar. I sow and discussed with my professors doing simulations some times: they wrote down a lot of calculus, most of them using Crank-Nicolson Method and so implement them on MATLAB.</p> <p>Until now, I could not imagine that MATLAB has pdepe. Today I've learned about pdepe and try to code a PDE that my professor coded with Crank-Nicolson. The result was absolutely the same and I had less calculus.</p> <p>I am thinking about the difference between the two methods. In fact, I'd like to know when is not so good use pdepe (once I thought this is marvelous, very easy...!). What can be the benefits of other methods on simple researchs?</p> <p>Many thanks in advance.</p> https://scicomp.stackexchange.com/q/34219 0 How to compute over 1 billion particles? Артур Клочко https://scicomp.stackexchange.com/users/33864 2020-01-16T16:18:49Z 2020-01-16T20:48:55Z <p>I want to simulate human erythrocytes in capillaries. I calculated, that for 1 meter long and 1 mm in diameter capillary there are about 3 billions blood cells.</p> <p>Erythrocytes are actually discs, but let’s assume that it’s a simple cuboid. What do we have</p> <pre><code>Erythrocyte -&gt; vertex -&gt; double x double y double z </code></pre> <p>In c++ double weights 8 bytes. 1 vertex=24 bytes, 8 vertexes = 24 bytes = data size of an erythrocyte. 24*7 billions bytes = 168 GB. I have only 4 Gb RAM, unfortunately. What should I do? I have never been up in computing.</p> <p>Try to decrease calculations?</p> <p>Actually there is a wide distribution of vessels diameter, for example the largest is the aorta(as far as I know) - few cm, and the shortest could have diameter about single erythrocyte size ~ 10um, by the way 10um it’s only 100th part of mm, so there is no sense to decreasing vessel diameter. </p> <p>Decrease vessel’s width? Okay, let it be 0.3 m = 30 cm. Recalculate:</p> <p>Assume vessel is a cylinder, with diameter=0.001 m and height=0.3 m. The volume equals <span class="math-container">$2.35619449*10^{11} \mu m^3$</span>. Average erythrocyte volume <span class="math-container">$90 \mu m^3$</span>.</p> <p>We have 2.6 billions erythrocytes, and 2.6 billions * 24 bytes of data. I still need some additional power computing solutions. Don’t I?</p> https://scicomp.stackexchange.com/q/34216 4 Why the solid FEM problem can not be solved after constraining 3 degrees of freedom? Xu Hui https://scicomp.stackexchange.com/users/33756 2020-01-16T08:58:09Z 2020-01-16T20:21:22Z <p>I write a simple MATLAB code for solving solid FEM problem.</p> <p>The problem looks like that</p> <pre><code>(1) (2) x-------x | / | | / | | / | x-------x (3) (4) </code></pre> <p>Each node has two degrees of freedom(DOF), which stand for <span class="math-container">$x$</span> and <span class="math-container">$y$</span> displacements.</p> <h1>Boundary condition</h1> <p>I want only constraint <code>node 1</code>'s <span class="math-container">$x,y$</span> DOF and <code>node 2</code>'s $x DOF, that is constraining DOFs <code>1,2,3</code> to zero. </p> <h1>Load condition</h1> <p>All nodes have a node force equal to 1.</p> <h1>Result</h1> <p>I found that the global stiffness matrix has a rank 7 but not the rank 8, which leads to the problem can not be solved.</p> <h1>Problem</h1> <p>In the 2D solid FEM case, the 3 DOFs are needed to be constrained for limiting rigid body translation along <span class="math-container">$x,y$</span> axis, and the rotaion. But after I constrain 3 DOFs, the problem still can not be solved. Why?</p> <h1>Appendix</h1> <p>The matlab code looks like</p> <pre><code>clc clear close all %% model x_range=10; y_range=10; num_x_node=2; num_y_node=2; num_node=num_x_node*num_y_node; DOF=2; %% generate grid node_coordinate_table=calculateNodeCoordinateTableForRectangleDomain(num_x_node,num_y_node,x_range,y_range); element_node_table=calculateElementNodeTableForRectangleDomainTriangleGrid(num_x_node,num_y_node); %% calculate element matrix K=calcGlobalMatrix(node_coordinate_table,element_node_table,DOF); %% load P=ones(num_node*DOF,1); %% constrain is_fix_dof=zeros(num_node*DOF,1); is_fix_dof(1)=1; is_fix_dof(2)=1; is_fix_dof(3)=1; %% modify discrete equation by constrain (cross zero center one method) for i=1:1:num_node*DOF if is_fix_dof(i)==1 K(i,:)=0; K(:,i)=0; K(i,i)=1; P(i)=0; end end %% solution x=K\P %% plot grid figure set(gcf,'position',[50,50,1000,800]) hold on plotNode(node_coordinate_table); plotElement(element_node_table,node_coordinate_table); %% function function K=calcGlobalMatrix(node_coordinate_table,element_node_table,DOF) num_node=size(node_coordinate_table,1); K=zeros(DOF*num_node); element_matrix_array=calculateElementMatrixArray(node_coordinate_table,element_node_table); num_element=size(element_node_table,1); for i=1:1:num_element K0=element_matrix_array{i}; node_global_index_array=element_node_table(i,:); for node_i=1:1:3 for node_j=1:1:3 for i1=1:1:DOF for j1=1:1:DOF node_global_index_i=node_global_index_array(node_i); node_global_index_j=node_global_index_array(node_j); row=(node_i-1)*DOF+i1; col=(node_j-1)*DOF+j1; ROW=(node_global_index_i-1)*DOF+i1; COL=(node_global_index_j-1)*DOF+j1; K(ROW,COL)=K(ROW,COL)+K0(row,col); end end end end end end function [element_matrix_array]=calculateElementMatrixArray(node_coordinate_table,element_node_table) num_element=size(element_node_table,1); element_matrix_array=cell(num_element,1); for i=1:1:num_element node_a_id=element_node_table(i,1); node_b_id=element_node_table(i,2); node_c_id=element_node_table(i,3); x=zeros(3,1); x(1)=node_coordinate_table(node_a_id,1); x(2)=node_coordinate_table(node_b_id,1); x(3)=node_coordinate_table(node_c_id,1); y=zeros(3,1); y(1)=node_coordinate_table(node_a_id,2); y(2)=node_coordinate_table(node_b_id,2); y(3)=node_coordinate_table(node_c_id,2); element_matrix_array{i}=calculateElementMatrixByNodeCoordinate(x,y); end end function plotNode(node_coordinate_table) x=node_coordinate_table(:,1); y=node_coordinate_table(:,2); plot(x,y,'o'); num_node=size(node_coordinate_table,1); for i=1:1:num_node node_x=x(i); node_y=y(i); text(node_x,node_y,num2str(i)); end end function plotElement(element_node_table,node_coordinate) num_element=size(element_node_table,1); for i=1:1:num_element node_a_id=element_node_table(i,1); node_b_id=element_node_table(i,2); node_c_id=element_node_table(i,3); x_a=node_coordinate(node_a_id,1); y_a=node_coordinate(node_a_id,2); x_b=node_coordinate(node_b_id,1); y_b=node_coordinate(node_b_id,2); x_c=node_coordinate(node_c_id,1); y_c=node_coordinate(node_c_id,2); plot([x_a,x_b,x_c,x_a],[y_a,y_b,y_c,y_a],'k-'); x_o=(x_a+x_b+x_c)/3; y_o=(y_a+y_b+y_c)/3; text(x_o,y_o,[num2str(i)]);% -1 for plot element index end end function node_coordinate_table=calculateNodeCoordinateTableForRectangleDomain(num_x_node,num_y_node,x_range,y_range) num_node=num_x_node*num_y_node; num_x_element=num_x_node-1; num_y_element=num_y_node-1; node_coordinate_table=zeros(num_node,2); n=1; dx=x_range/num_x_element; dy=y_range/num_y_element; for j=1:1:num_y_node for i=1:1:num_x_node x = (i-1) * dx; y = y_range - (j-1) * dy; node_coordinate_table(n,1)=x; node_coordinate_table(n,2)=y; n=n+1; end end end function element_node_table=calculateElementNodeTableForRectangleDomainTriangleGrid(num_x_node,num_y_node) % element node table generate is order in matlab format, that is, id start from 1 num_x_element=num_x_node-1; num_y_element=num_y_node-1; num_element=num_x_element*num_y_element*2; element_node_table = zeros(num_element,3); n=1; for j=1:1:num_y_element for i=1:1:num_x_element % a--c % | /| % |/ | % b--d node_a_id=i+(j-1)*num_x_node; node_b_id=node_a_id+num_x_node; node_c_id=node_a_id+1; node_d_id=node_b_id+1; % upper element abc element_node_table(n,1)=node_a_id; element_node_table(n,2)=node_b_id; element_node_table(n,3)=node_c_id; n=n+1; % lower element cbd element_node_table(n,1)=node_c_id; element_node_table(n,2)=node_b_id; element_node_table(n,3)=node_d_id; n=n+1; end end end function K0=calculateElementMatrixByNodeCoordinate(x,y) % reference % finite element method, wang xucheng chapter 2 % model E_0=100; nu_0=0; Thickness=1; % calculate element matrix a=zeros(3,1); b=zeros(3,1); c=zeros(3,1); for i=1:1:3 if(i==1);j=2;m=3;end if(i==2);j=3;m=1;end if(i==3);j=1;m=2;end a(i)=x(j)*y(m)-x(m)*y(j); b(i)=y(j)-y(m); c(i)=-x(j)+x(m); end A=1/2*(b(i)*c(j)-b(j)*c(i)); K0=zeros(6,6); for r=1:1:3 for s=1:1:3 K0(2*r-1,2*s-1)=b(r)*b(s)+(1-nu_0)/2*c(r)*c(s); K0(2*r-1,2*s)=nu_0*b(r)*c(s)+(1-nu_0)/2*c(r)*b(s); K0(2*r,2*s-1)=nu_0*c(r)*b(s)+(1-nu_0)/2*b(r)*c(s); K0(2*r,2*s)=c(r)*c(s)+(1-nu_0)/2*b(r)*b(s); end end K0=K0*E_0*Thickness/(4*A*(1-nu_0^2)); end <span class="math-container">`</span> </code></pre> https://scicomp.stackexchange.com/q/34213 29 why is A*v+B*v faster than (A+B)*v? Sam C. https://scicomp.stackexchange.com/users/26743 2020-01-15T21:54:42Z 2020-01-16T15:42:04Z <p><span class="math-container">$A$</span> and <span class="math-container">$B$</span> are <span class="math-container">$n \times n$</span> matrices and <span class="math-container">$v$</span> is a vector with <span class="math-container">$n$</span> elements. <span class="math-container">$Av$</span> has <span class="math-container">$\approx 2n^2$</span> flops and <span class="math-container">$A+B$</span> has <span class="math-container">$n^2$</span> flops. Following this logic, <span class="math-container">$(A+B)v$</span> should be faster than <span class="math-container">$Av+Bv$</span>.</p> <p>Yet, when I run the following code in matlab</p> <pre class="lang-py prettyprint-override"><code>A = rand(2000,2000); B = rand(2000,2000); v = rand(2000,1); tic D=zeros(size(A)); D = A; for i =1:100 D = A + B; (D)*v; end toc tic for i =1:100 (A*v+B*v); end toc </code></pre> <p>The opposite is true. A<em>v+B</em>v is over twice as fast. Any explanations?</p> https://scicomp.stackexchange.com/q/34212 1 Is the similar subdivision of a delaunay mesh still delaunay? EMP https://scicomp.stackexchange.com/users/9446 2020-01-15T18:15:42Z 2020-01-18T13:52:42Z <p>I have a delaunay triangulation for a 2d box with say an airfoil inside. If I uniformly refine this mesh by subdividing each triangle in the mesh into 4 triangles by halving each edge, is the resultant mesh still delaunay? I have been assuming so, but I have no proof this is the case. </p> https://scicomp.stackexchange.com/q/34208 -1 How to compute the gradient of T with Armadillo library [closed] ztdep https://scicomp.stackexchange.com/users/22404 2020-01-15T10:27:38Z 2020-01-17T00:20:29Z <p>I am using the Armadillo library to solve a 3d heat conduction problem on 3d unstructured grid system,the gradient of the T field is determined by the least square method. I have created a matrix equation for each Tet cell which has the following dimension A(4*3)*X(3)=B(4)</p> <p>Could you please give me some help about how to solve it with the least square method .</p> https://scicomp.stackexchange.com/q/34209 0 Is it possble to do this complex symbolic calculation with Matlab? Xinxin Peng https://scicomp.stackexchange.com/users/28691 2020-01-15T10:22:37Z 2020-01-15T12:15:15Z <p>Sorry it's bit abrupt, but recently I am caught up in some symbolic calcualtion which is tedious and almost impossible with mere human hands, so just wondering is it possible to solve the double integral of what's written below using Matlab or some other platform?</p> <p><span class="math-container">$$\int_{0}^{a}dk\frac{1}{k}\int_{|k-q|}^{|k+q|}dk'\frac{1}{k'}(k^2+q^2-k'^2)(k'^2+q^2-k^2)(k'^2-(k+q)^2)(k'^2-(k-q)^2)[\frac{1}{v_F(k-k')-w}+\frac{1}{v_F(k-k')+w}].$$</span></p> <p>And to clarify, all the variables are positive real numbers. The first integral is about <span class="math-container">$k'$</span> with upper and lower limit being <span class="math-container">$|k+q|$</span> and |k-q| respectively, and second integral is about <span class="math-container">$k$</span> from <span class="math-container">$0$</span> to <span class="math-container">$a$</span>. The symbols other than <span class="math-container">$k$</span> and <span class="math-container">$k'$</span> can all be seen as real positive constants.</p> https://scicomp.stackexchange.com/q/34206 -1 Change value of dependent sweep varietals Love Eva https://scicomp.stackexchange.com/users/33832 2020-01-15T06:48:16Z 2020-01-15T06:48:16Z <p>I am using Comsol to model a really hard problem. I am using the sweep for one variable (width), and as the problem state (Length=2*width).</p> <p>When I use the sweep the value of Length does not update. Is there a way to update the value of Length as it sweep the values of Width?</p> https://scicomp.stackexchange.com/q/34202 12 Example where autodiff works but symbolic differentiation will not? Lucas Roberts https://scicomp.stackexchange.com/users/33830 2020-01-15T02:46:16Z 2020-01-16T03:36:02Z <p>According to the survey paper on autodiff (linked) <a href="https://arxiv.org/pdf/1502.05767.pdf" rel="noreferrer">Autodiff works on inputs that cannot be specified in closed form but can be described by a sequence of code</a>, each component of which is differentiable. Autodiff also works on code that can be symbolically differentiated but in this latter case the benefits are less obvious and more subtle. However, I haven't been able to come up with an example where autodiff works but symbolic diff will not work. </p> <p>My question is: does there exist a simple example of autodiff with a code input that will not work with symbolic diff? </p> <p>Note: I realize simple is somewhat arbitrarily defined so let's say in less than 20 lines of code so that the answer isn't too long to read. If 20 lines isn't enough then something like the minimal number of lines in code would work. </p> https://scicomp.stackexchange.com/q/34200 4 Computation of triple nested loops as a convolution product? HansimGlück https://scicomp.stackexchange.com/users/33824 2020-01-14T13:31:57Z 2020-01-17T09:56:28Z <p>I'm trying to compute efficiently the following <span class="math-container">\begin{equation} A_j = \sum_{l'=1}^{\infty}\sum_{k= 0}^{K-1} L_{l'}T_ke^{2\pi i \frac{k}{K}j}\epsilon_{l',k} \end{equation}</span> for <span class="math-container">$j = 0,1, \ldots, K-2,K-1$</span>, where <span class="math-container">$L_l$</span> (dimension <span class="math-container">$NL$</span>) and <span class="math-container">$T_k$</span> (dimension <span class="math-container">$W$</span>) are 1D vectors of complex numbers and <span class="math-container">$\epsilon$</span> is a 2D matrix (dimensions : <span class="math-container">$L \times W$</span>) of complex numbers. </p> <p>Numerically, I'm performing the summation over <span class="math-container">$l'$</span> from <span class="math-container">$1$</span> to <span class="math-container">$NL$</span>, where <span class="math-container">$N$</span> is an integer that has to be chosen sufficiently large (to satisfy conditions unimportant in this context), therefore yielding <span class="math-container">\begin{equation} A_j = \sum_{l'=1}^{NL}\sum_{k= 0}^{K-1} L_{l'}T_ke^{2\pi i \frac{k}{K}j}\epsilon_{l'\text{mod} L,k}, \end{equation}</span> a truncated summation over <span class="math-container">$l'$</span> and the index of <span class="math-container">$\epsilon_{l'\text{mod} L,k}$</span> being taken with the modulo.</p> <p>While a straightforward (but naive) implementation of <span class="math-container">$A_j$</span> in terms of three nested for loops (two for the above double sum and a third that ranges over all <span class="math-container">$j$</span>) works fine, it is terribly slow, especially considering that this implementation is to be repeated for every time step (and there are many of them).</p> <p>I'm thus looking for an efficient implementation of the above, which looks very much like a 2D separable convolution product. This would allow one to make use of FFT's to compute it quite efficiently, taking advantage that in the frequency domain, a convolution becomes a multiplication of complex numbers.</p> <p>Does anybody have any idea ?</p> https://scicomp.stackexchange.com/q/34199 1 Can Representation Theory be studied computationally / numerically? user33821 https://scicomp.stackexchange.com/users/33821 2020-01-14T08:08:34Z 2020-01-14T08:08:34Z <p>Can a subfield such as the representation theory of Lie algebras be studied computationally / numerically -- is there an interplay between the abstract and the concrete? I would be grateful for an example, e.g. a reference to a journal paper.</p> https://scicomp.stackexchange.com/q/34197 2 Solution method of nonlinear heat transfer analysis vydesaster https://scicomp.stackexchange.com/users/28358 2020-01-13T22:08:09Z 2020-01-14T02:34:49Z <p>The governing equation of transient heat transfer analysis is described as follows: <span class="math-container">$$C \frac{dT}{dt}+K T = Q$$</span> <br> When using backward difference scheme for the discretization of the time we get the following equation: <span class="math-container">$$(\frac{C}{\Delta t}+K)T_n = Q_n+\frac{C}{\Delta t}T_{n-1}$$</span> Where <span class="math-container">$C$</span> is the heat capacity matrix, <span class="math-container">$K$</span> is the conductivity matrix, <span class="math-container">$T_n$</span> is the temperature for time step <span class="math-container">$n$</span>, <span class="math-container">$T_{n-1}$</span> is the temperature for time step <span class="math-container">$n-1$</span> and <span class="math-container">$Q$</span> is a heat flux vector. <br> The term becomes nonlinear as soon as I implement the radiation boundary condition. <br> My question now is, which methods can be used to solve the above nonlinear problem? Any further reading would be helpful.</p> https://scicomp.stackexchange.com/q/34194 1 How to vectorize 2 nested for loop with one complex condition in the inner loop? Jones G https://scicomp.stackexchange.com/users/16639 2020-01-13T20:44:49Z 2020-01-13T21:45:56Z <p>Octave calculations is too slow specially when you deal with scientific calculations that can span a very very large matrix, even just iterating. I have to vectorize it to speed up the calculation.</p> <p>Vectorizing a 1xn matrix would just be easy by following</p> <pre><code>for i = 1:n if (a(i) &gt; 5) a(i) -= 20 endif endfor </code></pre> <p><a href="https://octave.org/doc/v4.0.1/Basic-Vectorization.html" rel="nofollow noreferrer">https://octave.org/doc/v4.0.1/Basic-Vectorization.html</a></p> <p>But I dont have any clue when the matrix becomes nxm. My sample code is shown below.</p> <pre><code>segs = 1000000; r_int = 100; r_ot = 2000; x = linspace (0, 1000, segs); y = linspace (0, -1000, segs); [xx, yy] = meshgrid (x, y); circ = xx.*yy; circ_matrix = zeros(segs,segs); #this needs vectorization for j = 1:segs for i = 1:segs if((r_int&lt;=circ(i,j)) &amp;&amp;(circ(i,j)&lt;=r_ot)) circ_matrix(i,j)=1; endif; endfor; endfor; </code></pre> <p>I have marked the code that needs vectorization</p> <p>Some things tried:</p> <pre><code>circ_matrix((r_int&lt;=circ) &amp;&amp;(circ&lt;=r_ot))=1; </code></pre> <p>but did not work</p> https://scicomp.stackexchange.com/q/34185 -1 Simulating an object in orbit joethemow https://scicomp.stackexchange.com/users/33808 2020-01-13T00:43:50Z 2020-01-15T22:10:40Z <p>This question is more oriented around suggestions for simulation tools and how to approach simulating an object in orbit. So high level I am trying to simulate the concept of a <a href="https://en.wikipedia.org/wiki/Skyhook_(structure)" rel="nofollow noreferrer">Sky Hook</a>. What are the best softwares to run simulations to validate the dynamics of the structure and visually show the structure in motion? My best guess is MatLab and maybe Unity? Any direction is helpful!</p> https://scicomp.stackexchange.com/q/34179 3 Levinson Recursion for Non Square Toeplitz Matrices Royi https://scicomp.stackexchange.com/users/7951 2020-01-11T21:02:34Z 2020-01-17T06:58:31Z <p>Given a rectangular <a href="https://en.wikipedia.org/wiki/Toeplitz_matrix" rel="nofollow noreferrer">Toeplitz Matrix</a> <span class="math-container">$ H $</span>, how could one solve:</p> <p><span class="math-container">$$y = H x$$</span></p> <p>For instance, <span class="math-container">$ H $</span> can be Linear Convolution Matrix of the filter <span class="math-container">$ h $</span>:</p> <p><span class="math-container">$$H = \begin{bmatrix} {h}_{1} &amp; 0 &amp; 0 &amp; \ldots &amp; &amp; 0 \\ {h}_{2} &amp; {h}_{1} &amp; 0 &amp; \ldots &amp; &amp; 0 \\ {h}_{m} &amp; \ldots &amp; {h}_{1} &amp; 0 &amp; \ldots &amp; 0 \\ 0 &amp; {h}_{m} &amp; \ldots &amp; {h}_{1} &amp; 0 &amp; \ldots &amp; 0 \\ \vdots &amp; &amp; \ddots &amp; &amp; &amp; \vdots \\ 0 &amp; \ldots &amp; {h}_{m} &amp; {h}_{m - 1} &amp; \ldots &amp; {h}_{1} \\ 0 &amp; \ldots &amp; 0 &amp; {h}_{m} &amp; \ldots &amp; {h}_{2} \\ \vdots &amp; &amp; \ddots &amp; &amp; &amp; \vdots \\ 0 &amp; \ldots &amp; &amp; 0 &amp; \ldots &amp; {h}_{m} \\ \end{bmatrix}$$</span></p> <p>Which is clearly not square.<br> In the general case it would be generated by:</p> <pre class="lang-matlab prettyprint-override"><code>numRows = 10; %&lt;! Or any other number numCols = 20; %&lt;! Or any other number vR = randn(numRows, 1); vC = randn(numCols, 1); mH = toeplitz(vC, vR); </code></pre> <p>Since the Matrix is not square, I'm after the least squares solution:</p> <p><span class="math-container">$$\arg \min_{x} {\left\| H x - y \right\|}_{2}^{2}$$</span></p> <p>I have implemented the <a href="https://github.com/RoyiAvital/Projects/tree/master/LevinsonRecursion" rel="nofollow noreferrer">Levinson Recursion</a> for the square case yet I'd like to know if there is an extension to it.</p> <p>At the moment the trick I use for this specific case it to partition the matrix and data for a subset of square Toeplitz Matrices and use the algorithm I implemented.<br> I wonder if there is a better way.</p> https://scicomp.stackexchange.com/q/33284 2 Givens rotation vs 2x2 Householder reflection gTcV https://scicomp.stackexchange.com/users/21342 2019-08-21T06:15:55Z 2020-01-18T11:02:40Z <p>The usual story of Givens rotations vs Householder reflections is that Householder reflections are better if you want to map a long vector to <span class="math-container">$e_1$</span>, while Givens is better if you want to map a 2-vector to <span class="math-container">$e_1$</span>. However, I can't find a reference which explains why Givens is better in the 2-vector case. Can anyone explain?</p> https://scicomp.stackexchange.com/q/30321 1 What is the difference between Abaqus and Calculix contact input? hamza boulahia https://scicomp.stackexchange.com/users/26742 2018-10-09T11:38:20Z 2020-01-17T05:53:52Z <p>I would like to say first that am new at using Calculix.</p> <p>I'm using Abaqus/CAE to create a cup deep drawing simulation and everything worked perfectly but my objective is to run the same exact simulation on Calculix and compare the results</p> <p>So I used the same geometry and mesh and converted all the keywords from the input file of Abaqus to match the format of the input file for Calculix (it is a simple simulation so no big changes are made).</p> <p>but when I run the CCX solver it stops at <span class="math-container">$t=0.09$</span> where the total time is <span class="math-container">$t=1$</span>, and the error message is: </p> <blockquote> <p>Too many cutbacks.</p> </blockquote> <p>I tried:</p> <ul> <li>refining the mesh of the slave part</li> <li>lowering the initial time increment</li> </ul> <p>I don't know where is the problem although it seems to me that is in the contact definition because every other thing defined is very simple (like mesh, boundary condition, and surface definition etc...) and no complicated parameters are needed.</p> <p>This is how I defined the contact:</p> <pre><code> ** contact *Contact Pair, Interaction=SI1, Type=Surface To Surface Slave_tole, Master_poin *Surface Interaction, Name=SI1 *Surface Behavior, Pressure-Overclosure=Linear 1290000 *Friction 0.25, 350000 </code></pre> <p>I used the same values as in the Abaqus simulation.</p> https://scicomp.stackexchange.com/q/21417 5 Eigen - store sparse matrix as binary physicsGuy https://scicomp.stackexchange.com/users/18362 2015-11-27T18:33:51Z 2020-01-15T09:09:42Z <p>I need to store large sparse matrices in Eigen. I cannot find anything in the library except the function below, in <a href="http://eigen.tuxfamily.org/dox/unsupported/MarketIO_8h_source.html" rel="noreferrer">Eigen/Unsupported</a>. The problem with saveMarket is, that it saves in text format. Due to the size of my matrices I need to store my sparse matrices as binaries. Is there an easy way to adjust the function below to store as a binary. And an easy way to reload the matrix?</p> <pre><code>template&lt;typename SparseMatrixType&gt; bool saveMarket(const SparseMatrixType&amp; mat, const std::string&amp; filename, int sym = 0) { typedef typename SparseMatrixType::Scalar Scalar; std::ofstream out(filename.c_str(),std::ios::out); if(!out) return false; out.flags(std::ios_base::scientific); out.precision(64); std::string header; internal::putMarketHeader&lt;Scalar&gt;(header, sym); out &lt;&lt; header &lt;&lt; std::endl; out &lt;&lt; mat.rows() &lt;&lt; " " &lt;&lt; mat.cols() &lt;&lt; " " &lt;&lt; mat.nonZeros() &lt;&lt; "\n"; int count = 0; for(int j=0; j&lt;mat.outerSize(); ++j) for(typename SparseMatrixType::InnerIterator it(mat,j); it; ++it) { ++ count; internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out); // out &lt;&lt; it.row()+1 &lt;&lt; " " &lt;&lt; it.col()+1 &lt;&lt; " " &lt;&lt; it.value() &lt;&lt; "\n"; } out.close(); return true; } </code></pre> https://scicomp.stackexchange.com/q/20644 9 Use of machine learning in computational fluid dynamics EngrStudent - Reinstate Monica https://scicomp.stackexchange.com/users/3979 2015-09-04T11:22:27Z 2020-01-15T19:55:16Z <p><strong>Background:</strong><br> I have only built one working numeric solution to 2d Navier-Stokes, for a course. It was a solution for lid-driven cavity flow. The course, however, discussed a number of schemas for spatial discretizations and time discretizations. I have also taken more symbol-manipulation coursework applied to NS.</p> <p>Some of the numeric approaches to handle conversion of the analytic/symbolic equation from PDE to finite difference include:</p> <ul> <li>Euler FTFS, FTCS, BTCS</li> <li>Lax</li> <li>Midpoint Leapfrog</li> <li>Lax-Wendroff</li> <li>MacCormack</li> <li>offset grid (spatial diffusion allows information to spread)</li> <li>TVD </li> </ul> <p>To me, at the time, these seemed like "insert-name finds a scheme and it happens to work". Many of these were from before the time of "plentiful silicon". They are all approximations. In the limit they. in theory, lead to the PDE's. </p> <p>While Direct Numerical Simulation (<a href="https://en.wikipedia.org/wiki/Direct_numerical_simulation" rel="nofollow noreferrer">DNS</a>) is fun, and Reynolds Averaged Navier-Stokes (<a href="https://en.wikipedia.org/wiki/Reynolds-averaged_Navier%E2%80%93Stokes_equations" rel="nofollow noreferrer">RANS</a>) is also fun, they are the two "endpoints" of the continuum between computationally tractable, and fully representing the phenomena. There are multiple families of approaches that live interior to these. </p> <p>I have had CFD professors say, in lecture, that most CFD solvers make pretty pictures, but for the most part, those pictures do not represent reality and that it can be very tough, and take lots of work, to get a solver solution that does represent reality. </p> <p>The sequence of development (as I understand it, not exhaustive) is:</p> <ol> <li>start with the governing equations -> PDE's</li> <li>determine your spatial and temporal discretization -> grid and FD rules</li> <li>apply to the domain including initial conditions and boundary conditions </li> <li>solve (lots of variations on matrix inversion)</li> <li><p>perform gross reality checks, fit to known solutions, etc..</p></li> <li><p>build some simpler physical models derived from analytic results</p></li> <li>test them, analyze, and evaluate</li> <li>iterate (jumping back to either step 6, 3, or 2)</li> </ol> <p>Thoughts:<br> I have recently been working with CART models, oblique trees, random forests, and gradient boosted trees. They follow more mathematically derived rules, and the math drives the shape of the tree. They work to make discretized forms well. </p> <p>Although these human-created numeric approaches work somewhat, there is extensive "voodoo" needed to connect their results to the physical phenomena they are meant to model. Often the simulation does not substantially replace real-world testing and verification. It is easy to use the wrong parameter, or not account for variation in geometry or application parameters experienced in the real world.</p> <p><strong>Questions:</strong> </p> <ul> <li>Has there been any approach to let the nature of the problem define<br> the appropriate discretization, spatial and temporal differencing scheme, initial conditions, or solution?</li> <li>Can a high definition solution coupled with the techniques of machine learning be used to make a differencing scheme that has much larger step sizes but retains convergence, accuracy, and such?</li> <li>All of these schemes are accessibly "humanly tractable to derive" - they have a handful of elements. Is there a differencing scheme with thousands of elements that does a better job? How is it derived?</li> </ul> <p>Note: I will follow up with the empirically intialized and empirically derived (as opposed to analytically) in a separate question. </p> <p><strong>UPDATE:</strong></p> <ol> <li><p>Use of deep learning to accelerate lattice Boltzmann flows. Gave ~9x speedup for their particular case</p> <p>Hennigh, O. (in press). Lat-Net: Compressed Lattice Boltzmann Flow Simulations using Deep Neural Networks. Retrieved from: <a href="https://arxiv.org/pdf/1705.09036.pdf" rel="nofollow noreferrer">https://arxiv.org/pdf/1705.09036.pdf</a></p> <p>Repo with code (I think):<br> <a href="https://github.com/loliverhennigh/Phy-Net" rel="nofollow noreferrer">https://github.com/loliverhennigh/Phy-Net</a></p></li> <li><p>About 2 orders of magnitude faster than GPU, 4 orders of magnitude, or ~O(10,000x) faster than CPU, and same hardware. </p> <p>Guo, X., Li, W. &amp; Ioiro, F. Convolutional Neural Networks for Steady Flow Approximation. Retrieved from: <a href="https://autodeskresearch.com/publications/convolutional-neural-networks-steady-flow-approximation" rel="nofollow noreferrer">https://autodeskresearch.com/publications/convolutional-neural-networks-steady-flow-approximation</a></p></li> <li><p>Others who have looked into the topic about 20 years ago:</p> <p>Muller, S., Milano, M. &amp; Koumoutsakos P. Application of machine learning algorithms to flow modeling and optimization. Center for Turbulence Research Annual Research Briefs 1999 Retrieved from: <a href="https://web.stanford.edu/group/ctr/ResBriefs99/petros.pdf" rel="nofollow noreferrer">https://web.stanford.edu/group/ctr/ResBriefs99/petros.pdf</a> </p></li> </ol> <p><strong>Update (2017):</strong><br> <a href="https://arxiv.org/abs/1712.06567" rel="nofollow noreferrer">This</a> characterises the use of non-gradient methods in deep learning, an arena which has been exclusively gradient based. While the direct implication of activity is in deep learning, it also suggests that GA can be used as an equivalent in solving a very hard, very deep, very complex problem at the level consistent with or superior to gradient descent based methods. </p> <p>Within the scope of this question, it might suggest that a larger-scale, machine-learning based attack might allow "templates" in time and space that substantially accelerate convergence of gradient-domain methods. The article goes as far as to say that sometimes going in the direction of gradient descent moves away from the solution. While in any problem with local optima or pathological trajectories (most high-value real-world problems have some of these) it is expected that the gradient isn't globally informative, it is still nice to have it quantified and validated empirically as it was in this paper and the ability to "jump the bound" without requiring "reduction of learning" as you get in momentum or under-relaxation.</p> <p><strong>Update (2019):</strong><br> It seems that google now has a contribution "how to find a better solver" piece of the AI puzzle. <a href="https://ai.googleblog.com/2018/08/introducing-new-framework-for-flexible.html" rel="nofollow noreferrer">link</a> This is a part of making the AI make the solver.</p> <p>**Update (2020): ** And now they are doing it, and doing it well...<br> <a href="https://arxiv.org/pdf/1911.08655.pdf" rel="nofollow noreferrer">https://arxiv.org/pdf/1911.08655.pdf</a></p> <p>It could be argued that they could then deconstruct their NN to determine the actual discretization. I particularly like figure 4.</p> https://scicomp.stackexchange.com/q/19685 10 How to sample points in hyperbolic space? doetoe https://scicomp.stackexchange.com/users/5521 2015-05-19T13:33:05Z 2020-01-15T03:14:08Z <p>Hyperbolic space in the <a href="https://en.wikipedia.org/wiki/Poincar%C3%A9_half-plane_model" rel="noreferrer">Poincaré upper half space model</a> looks like ordinary <span class="math-container">$\Bbb R^n$</span> but with the notion of angle and distance distorted in a relatively simple way. In Euclidean space I can sample a random point uniformly in a ball in several ways, e.g. by generating <span class="math-container">$n$</span> independent Gaussian samples to obtain a direction, and separately sample a radial coordinate <span class="math-container">$r$</span> by uniformly sampling <span class="math-container">$s$</span> from <span class="math-container">$\left[0, \frac1{n+1}R^{n+1}\right]$</span>, where <span class="math-container">$R$</span> is the radius, and setting <span class="math-container">$r = \left((n+1)s\right)^{\frac1{n+1}}$</span>. In the hyperbolic upper half plane a sphere happens to still be a sphere, only its centre will not be the centre in the Euclidean metric, so we could do the same. </p> <p>If we want to sample according to a non-uniform distribution, but still in an isotropic way, e.g. a Gaussian distribution, this doesn't seem so easy. In Euclidean space we could just generate a Gaussian sample for each coordinate (this only works for the Gaussian distribution), or equivalently generate a multidimensional Gaussian sample. Is there a direct way to convert this sample to a sample in hyperbolic space? </p> <p>An alternative approach could be to first generate a direction uniformly distributed direction (e.g. from <span class="math-container">$n\$</span> Gaussian samples) then a Gaussian sample for the radial component, and finally generate the image under the <a href="https://en.wikipedia.org/wiki/Exponential_map_(Riemannian_geometry)" rel="noreferrer">exponential map</a> in the specified direction for the specified length. A variation would be to just take the Euclidean Gaussian sample and map it under the exponential map.</p> <p>My questions:</p> <ul> <li>what would be a good and efficient way to obtain a Gaussian sample with given mean and standard deviation in hyperbolic space?</li> <li>do the ways I describe above provide the desired sampling?</li> <li>did anyone work out the formula's already?</li> <li>how does this generalize to other metrics and other probability distributions?</li> </ul> <p>Thanks in advance.</p> <p><strong>EDIT</strong></p> <p>I just realized that even in the case of uniform sampling these questions remain; even though a sphere is a sphere, a uniform distribution would not be described by a constant function on a ball.</p>