Thanks Yves! So I gave up on Matlab for using mumps, and went back to the basic configuration of getFEM (without mumps, metis or mpi). But unfortunately I cannot replicate what I did before. The solver just exited in the middle of the solving (after ~1.5 hours of solving), with the following warning message:
Gmres is blocked, exiting Level 2 Warning in ../../src/getfem/getfem_model_solvers.h, line 105: gmres did not converge! The solver is by default the iterative solver. Do you have any idea of what this means? Thanks a lot and sorry for keep bugging you for silly questions! On Fri, Jul 7, 2017 at 6:13 PM, Yves Renard <[email protected]> wrote: > Dear Andy, > > Unfortunately, the parallel version using mpi is not compatible with > Matlab (may be the openmp one is, but I never tried). You should install > the sequential version of Mumps (without Metis) or use the python interface > which is compatible with the parallel version of Mumps (but in that case, > mpirun as to be used to run the python script, see the example in the > documentation). > > Best regards, > > Yves. > > ----- Original Message ----- > From: "Yu (Andy) Huang" <[email protected]> > To: [email protected] > Sent: Friday, July 7, 2017 11:51:36 PM > Subject: Re: [Getfem-users] Fwd: Simulating electric field distribution > > Thanks Yves and Konstantinos! > > The iterative solver was used in my case. So I tried to re-compile getFEM > with mumps and parallel computing enabled, and installed again (before this > I installed METIS and Mumps), but for some unknown reason it made Matlab > crash once the gf_model_get(md, 'solve') is executed. > > So I realized I must have messed up everything. Then I deleted getFEM and > downloaded the latest version of getFEM (v5.2) and re-installed it. Now it > always crashes Matlab when tries to solve the model (a segmentation > violation error was thrown out by Matlab). I did a "make check" and found > the check_all.sh failed. Do you know what might go wrong with my > installation? Below is the commands I used for installation: > > ./configure --enable-matlab --enable-mumps --with-pic > > make > > cp interface/src/matlab/gf_matlab.mexa64 interface/src/matlab/gf_matlab > > sudo make install > > sudo cp /usr/local/getfem_toolbox/gf_matlab > /usr/local/getfem_toolbox/gf_matlab.mexa64 > > > Thanks again for your help! > > > > On Thu, Jun 22, 2017 at 8:05 AM, Konstantinos Poulios < > [email protected]> wrote: > > > Dear Andy, > > > > Have you checked which linear solver is used by your: > > > > gf_model_get(md, 'solve'); > > > > ? > > > > If you have compiled GetFEM with mumps, the default behavior would be to > > use mumps for such large systems. Can you confirm? > > > > BR > > Kostas > > > > > > On Thu, Jun 22, 2017 at 1:31 AM, Yu (Andy) Huang <[email protected] > > > > wrote: > > > >> Thanks for your generous help, Yves!! Now it worked, and it also worked > >> on the real model! > >> > >> It took more than 4 hours to solve the real model (with about 1 million > >> tetrahedra elements). Do you know any tricks to speed up the solving? > I'm > >> thinking about parallelize the solving and was reading the > documentation. > >> It says I need to compile the package again with parallel option > enabled? > >> Also do you have any empirical values for the options in calling > >> gf_model_get(model,'solve')? e.g. the "max_iter" and "max_res"? > >> > >> Thanks again for your help! > >> > >> > >> > >> ---------- Forwarded message ---------- > >> From: Yves Renard <[email protected]> > >> Date: Tue, Jun 20, 2017 at 3:10 AM > >> Subject: Re: [Getfem-users] Simulating electric field distribution > >> To: "Yu (Andy) Huang" <[email protected]>, [email protected] > >> > >> > >> Dear Andy, > >> > >> I just add some commentaries on your code. > >> > >> Yves. > >> > >> > >> Le 19/06/2017 à 21:27, Yu (Andy) Huang a écrit : > >> > >> Dear getFEM users, > >> > >> Thanks for you previous help on my silly questions! Now I manage to > >> simulate the electric field on a toy sphere with two layers, each layer > >> having a different electrical conductivity. I'm just not sure if I did > it > >> properly, because when I compare the results to those I got from Abaqus > (a > >> commercial FEM solver), I see some difference that I don't understand > (see > >> attached screenshots). I used the *same mesh*, with the* same boundary > >> condition* and *same conductivity *values, but the distribution and > >> absolute values of voltage and field are all different between getFEM > and > >> Abaqus. My major concern is the way I coded the boundary conditions. The > >> problem I'm solving is a Laplacian equation of electric potential u, > with > >> the following BC: > >> > >> 1) injecting 1 A/m^2 current density at the north pole of the sphere: - > >> *n*.*J *= 1 > >> 2) ground at the south pole: V = 0 > >> 3) insulation at all outer boundary: *n*.*J* = 0 > >> 4) continuity for inner boundary: *n*.(*J1* - *J2*) = 0 > >> > >> I only coded explicitly (1) and (2) so not sure if it's good enough. I > >> put part of my code below (Matlab code). Any advice on the code is much > >> appreciated! > >> > >> % =================================================== > >> mesh = gfMesh('import','gmsh', 'toySphere.msh'); > >> > >> mfu = gf_mesh_fem(mesh, 1); % scalar-field (electric potential) > >> mfE = gf_mesh_fem(mesh, 3); % 3d vector-field (electric field) > >> > >> gf_mesh_fem_set(mfu, 'fem', gf_fem('FEM_PK(3,1)')); % P1 Lagrange > >> gf_mesh_fem_set(mfE, 'fem', gf_fem('FEM_PK(3,1)')); % P1 Lagrange > >> > >> mim = gf_mesh_im(mesh, gf_integ('IM_TETRAHEDRON(1)')); % integration > >> method > >> > >> The integration method is used for both the volumic term and the > boundary > >> conditions, so it has to be of order two > >> -> mim = gf_mesh_im(mesh, gf_integ('IM_TETRAHEDRON(2)')); > >> > >> > >> md=gf_model('real'); > >> gf_model_set(md, 'add fem variable', 'u', mfu); > >> > >> rid = gf_mesh_get(mesh,'regions'); > >> reg1 = gf_mesh_get(mesh, 'region', rid(1)); > >> reg2 = gf_mesh_get(mesh, 'region', rid(2)); > >> region1 = 1; > >> region2 = 2; > >> gf_mesh_set(mesh, 'region', region1, reg1); > >> gf_mesh_set(mesh, 'region', region2, reg2); > >> > >> Did you check that your regions were ok (for instance with gf_plot_mesh > ) > >> ? > >> > >> > >> % governing equation and conductivities > >> sigma = [0.276;0.126]; % conductivity values for the two layers in the > >> sphere > >> gf_model_set(md, 'add linear generic assembly brick', mim, > >> [num2str(sigma(1)) '*(Grad_u.Grad_Test_u)'],region1); > >> gf_model_set(md, 'add linear generic assembly brick', mim, > >> [num2str(sigma(2)) '*(Grad_u.Grad_Test_u)'],region2); > >> > >> fb1 = gf_mesh_get(mesh, 'outer faces with direction', [0 0 1], 0.1, > >> cvid(indAnode)); > >> fb2 = gf_mesh_get(mesh, 'outer faces with direction', [0 0 -1], 0.1, > >> cvid(indCathode)); > >> *% the boundary condition is injecting 1 A/m^2 current density at the > >> north pole of the sphere, with the south pole as ground* > >> *% here indAnode and indCathode is the index of the element > corresponding > >> to the north and south pole* > >> > >> > >> anode_area = 3; > >> cathode_area = 4; > >> gf_mesh_set(mesh, 'region', anode_area, fb1); > >> gf_mesh_set(mesh, 'region', cathode_area, fb2); > >> > >> > >> Did you check that your boundaries were ok (may be this is your first > >> graph ) ? > >> > >> > >> gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u', > >> mfu, cathode_area); > >> > >> gf_model_set(md, 'add initialized data','Jn', ones(gf_mesh_get(mesh, > >> 'nbpts'),1)); > >> % here is the line that I suspect mostly. I have to pass 'Jn' a vector, > >> otherwise it won't solve. > >> gf_model_set(md, 'add source term brick', mim, 'u', [num2str(sigma(1)) > >> '*(-Grad_u.Normal)'], anode_area, 'Jn'); > >> > >> The source term isthe right hand side of -J.n = f, so in your case just > >> 1. So you do not need any "gf_model_set(md, 'add initialized data','Jn', > >> ones(gf_mesh_get(mesh, 'nbpts'),1));" and you can simply add > >> > >> gf_model_set(md, 'add source term brick', mim, 'u', '1', anode_area); > >> > >> > >> > >> % Neumann BC of electric potential > >> > >> % solve > >> gf_model_get(md, 'solve'); > >> > >> % extracted solution > >> u = gf_model_get(md, 'variable', 'u'); > >> E = gf_model_get(md, 'interpolation', '-Grad_u', mfu); % electric field > >> > >> % display > >> figure; gf_plot(mfu, u, 'mesh','on','cvlst', get(mesh, 'outer > >> faces')); colormap(jet); colorbar > >> figure; gf_plot(mfE, E, 'mesh','on','norm','on','cvlst', get(mesh, > >> 'outer faces')); colormap(jet); colorbar > >> %=========================================================== > >> =================== > >> > >> > >> > >> > >> On Thu, Jun 8, 2017 at 10:55 PM, Yu (Andy) Huang < > [email protected]> > >> wrote: > >> > >>> Dear getFEM users, > >>> > >>> I'm entirely new to getFEM, and I'm trying to simulate the electric > >>> field distribution in the human brain when direct electric current is > >>> applied on the scalp surface. I know it's just a Laplacian equation of > the > >>> electric potential, and I managed to simulate the voltage distribution > on a > >>> toy (a cube). > >>> > >>> Now my question is: how do I simulate the electric field? should I add > >>> another variable of electric field? or can I just get the field from > the > >>> voltage solution? I tried both but without any luck. I added the > electric > >>> field as a new variable but did not figure out how to properly add > boundary > >>> condition using gf_model_set(). If calculating field from voltage, I > didn't > >>> find out which function to use to establish a relation between the > field > >>> variable and voltage variable. > >>> > >>> Any suggestion is appreciated! The examples in the documentation are > >>> generally mechanical problems, and there are very limited online > resources, > >>> so I really get stuck here. > >>> > >>> Thanks a lot! > >>> > >>> -- > >>> Yu (Andy) Huang, Ph.D. > >>> Postdoc fellow at Dept. of Biomedical Engineering, City College of New > >>> York > >>> Center for Discovery and Innovation, Rm. 3.320, > >>> <http://neuralengr.com/directions-to-cdi-center-for- > discovery-and-innovation/> > >>> 85 St Nicholas Terrace, New York, NY 10027 > >>> <http://neuralengr.com/directions-to-cdi-center-for- > discovery-and-innovation/> > >>> Tel: 1-646-509-8798 <%28646%29%20509-8798> > >>> Email: [email protected] > >>> *[email protected]* <[email protected]> > >>> http://www.parralab.org/people/yu-andy-huang/ > >>> > >> > >> > >> > >> -- > >> Yu (Andy) Huang, Ph.D. > >> Postdoc fellow at Dept. of Biomedical Engineering, City College of New > >> York > >> Center for Discovery and Innovation, Rm. 3.320, > >> <http://neuralengr.com/directions-to-cdi-center-for- > discovery-and-innovation/> > >> 85 St Nicholas Terrace, New York, NY 10027 > >> <http://neuralengr.com/directions-to-cdi-center-for- > discovery-and-innovation/> > >> Tel: 1-646-509-8798 <%28646%29%20509-8798> > >> Email: [email protected] > >> *[email protected]* <[email protected]> > >> http://www.parralab.org/people/yu-andy-huang/ > >> > >> > >> -- > >> > >> Yves Renard ([email protected]) tel : (33) > 04.72.43.87.08 > >> Pole de Mathematiques, INSA-Lyon fax : (33) 04.72.43.85.29 > >> 20, rue Albert Einstein > >> 69621 Villeurbanne Cedex, FRANCE > >> http://math.univ-lyon1.fr/~renard > >> > >> --------- > >> > >> > >> > >> > >> -- > >> Yu (Andy) Huang, Ph.D. > >> Postdoc fellow at Dept. of Biomedical Engineering, City College of New > >> York > >> Center for Discovery and Innovation, Rm. 3.320, > >> <http://neuralengr.com/directions-to-cdi-center-for- > discovery-and-innovation/> > >> 85 St Nicholas Terrace, New York, NY 10027 > >> <http://neuralengr.com/directions-to-cdi-center-for- > discovery-and-innovation/> > >> Tel: 1-646-509-8798 <%28646%29%20509-8798> > >> Email: [email protected] > >> *[email protected]* <[email protected]> > >> http://www.parralab.org/people/yu-andy-huang/ > >> <http://neuralengr.com/members/yu-%28andy%29-huang> > >> > > > > > > > -- > Yu (Andy) Huang, Ph.D. > Postdoc fellow at Dept. of Biomedical Engineering, City College of New York > Center for Discovery and Innovation, Rm. 3.320, > <http://neuralengr.com/directions-to-cdi-center-for- > discovery-and-innovation/> > 85 St Nicholas Terrace, New York, NY 10027 > <http://neuralengr.com/directions-to-cdi-center-for- > discovery-and-innovation/> > Tel: 1-646-509-8798 <(646)%20509-8798> > Email: [email protected] > *[email protected]* <[email protected]> > http://www.parralab.org/people/yu-andy-huang/ > <http://neuralengr.com/members/yu-%28andy%29-huang> > -- Yu (Andy) Huang, Ph.D. Postdoc fellow at Dept. of Biomedical Engineering, City College of New York Center for Discovery and Innovation, Rm. 3.320, <http://neuralengr.com/directions-to-cdi-center-for-discovery-and-innovation/> 85 St Nicholas Terrace, New York, NY 10027 <http://neuralengr.com/directions-to-cdi-center-for-discovery-and-innovation/> Tel: 1-646-509-8798 Email: [email protected] *[email protected]* <[email protected]> http://www.parralab.org/people/yu-andy-huang/ <http://neuralengr.com/members/yu-%28andy%29-huang>
