You can subscribe to this list here.
2003 
_{Jan}

_{Feb}

_{Mar}

_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}
(2) 
_{Oct}
(2) 
_{Nov}
(27) 
_{Dec}
(31) 

2004 
_{Jan}
(6) 
_{Feb}
(15) 
_{Mar}
(33) 
_{Apr}
(10) 
_{May}
(46) 
_{Jun}
(11) 
_{Jul}
(21) 
_{Aug}
(15) 
_{Sep}
(13) 
_{Oct}
(23) 
_{Nov}
(1) 
_{Dec}
(8) 
2005 
_{Jan}
(27) 
_{Feb}
(57) 
_{Mar}
(86) 
_{Apr}
(23) 
_{May}
(37) 
_{Jun}
(34) 
_{Jul}
(24) 
_{Aug}
(17) 
_{Sep}
(50) 
_{Oct}
(24) 
_{Nov}
(10) 
_{Dec}
(60) 
2006 
_{Jan}
(47) 
_{Feb}
(46) 
_{Mar}
(127) 
_{Apr}
(19) 
_{May}
(26) 
_{Jun}
(62) 
_{Jul}
(47) 
_{Aug}
(51) 
_{Sep}
(61) 
_{Oct}
(42) 
_{Nov}
(50) 
_{Dec}
(33) 
2007 
_{Jan}
(60) 
_{Feb}
(55) 
_{Mar}
(77) 
_{Apr}
(102) 
_{May}
(82) 
_{Jun}
(102) 
_{Jul}
(169) 
_{Aug}
(117) 
_{Sep}
(80) 
_{Oct}
(37) 
_{Nov}
(51) 
_{Dec}
(43) 
2008 
_{Jan}
(71) 
_{Feb}
(94) 
_{Mar}
(98) 
_{Apr}
(125) 
_{May}
(54) 
_{Jun}
(119) 
_{Jul}
(60) 
_{Aug}
(111) 
_{Sep}
(118) 
_{Oct}
(125) 
_{Nov}
(119) 
_{Dec}
(94) 
2009 
_{Jan}
(109) 
_{Feb}
(38) 
_{Mar}
(93) 
_{Apr}
(88) 
_{May}
(29) 
_{Jun}
(57) 
_{Jul}
(53) 
_{Aug}
(48) 
_{Sep}
(68) 
_{Oct}
(151) 
_{Nov}
(23) 
_{Dec}
(35) 
2010 
_{Jan}
(84) 
_{Feb}
(60) 
_{Mar}
(184) 
_{Apr}
(112) 
_{May}
(60) 
_{Jun}
(90) 
_{Jul}
(23) 
_{Aug}
(70) 
_{Sep}
(119) 
_{Oct}
(27) 
_{Nov}
(47) 
_{Dec}
(54) 
2011 
_{Jan}
(22) 
_{Feb}
(19) 
_{Mar}
(92) 
_{Apr}
(93) 
_{May}
(35) 
_{Jun}
(91) 
_{Jul}
(32) 
_{Aug}
(61) 
_{Sep}
(7) 
_{Oct}
(69) 
_{Nov}
(81) 
_{Dec}
(23) 
2012 
_{Jan}
(64) 
_{Feb}
(95) 
_{Mar}
(35) 
_{Apr}
(36) 
_{May}
(63) 
_{Jun}
(98) 
_{Jul}
(70) 
_{Aug}
(171) 
_{Sep}
(149) 
_{Oct}
(64) 
_{Nov}
(67) 
_{Dec}
(126) 
2013 
_{Jan}
(108) 
_{Feb}
(104) 
_{Mar}
(171) 
_{Apr}
(133) 
_{May}
(108) 
_{Jun}
(100) 
_{Jul}
(93) 
_{Aug}
(126) 
_{Sep}
(74) 
_{Oct}
(59) 
_{Nov}
(145) 
_{Dec}
(93) 
2014 
_{Jan}
(38) 
_{Feb}
(45) 
_{Mar}
(26) 
_{Apr}
(41) 
_{May}
(125) 
_{Jun}
(70) 
_{Jul}
(61) 
_{Aug}
(66) 
_{Sep}
(60) 
_{Oct}
(110) 
_{Nov}
(27) 
_{Dec}
(30) 
2015 
_{Jan}
(43) 
_{Feb}
(67) 
_{Mar}
(56) 
_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

From: Roy Stogner <roystgnr@ic...>  20150324 19:56:59

On Tue, 24 Mar 2015, Sylvain Vallaghe wrote: > CXX libopt_lafpoptimizer.lo > 329 1 sincos{ic=SINCOS}( (F64) *((F64*) a_16545_V$33b3), (F64 > *)sin_16545_V$33b1, (F64 *)cos_16545_V$33b2 ); [CALL_CONVENTION_UNIX_ABI] > 333 1 sincosf{ic=SINCOSF}( (F32) *((F32*) a_16546_V$33b6), (F32 > *)sin_16546_V$33b4, (F32 *)cos_16546_V$33b5 ); [CALL_CONVENTION_UNIX_ABI] > fpoptimizer/grammar_data.cc(81): internal error: assertion failed: > lower_dynamic_init_aggregate_constant: bad aggr kind > (shared/cfe/edgcpfe/lower_init.c, line 6341) > > const ParamSpec_NumConstant <Value_t> > plist_n_container<Value_t>::plist_n[20] = > ^ > > compilation aborted for > /home/akselos/software/libmeshsrc/contrib/fparser/fpoptimizer.cc (code 4) > > > Any ideas of what the problem could be ? "internal error" suggests a bug in icpc; you should probably report it to them. As a possible workaround, you might try configuring with "withfparser=devel" to see if that avoids triggering the bug. As a much more likely workaround, configure with "disablefparser" if you don't need it. > Also general comments on how to use intel tools properly with libmesh are > welcome (intel compilers, intel mkl, intel mpi). My apologies; I don't think I've used anything more recent than icpc 13.1 and a contemporaneous MKL.  Roy 
From: Sylvain Vallaghe <sylvain.vallaghe@gm...>  20150324 19:48:23

Hi all, I'm trying to compile libmesh with the intel compilers version 15.0.2 on Ubuntu 14.04.1 LTS. I'm running into the following error during compilation: CXX libopt_lafpoptimizer.lo 329 1 sincos{ic=SINCOS}( (F64) *((F64*) a_16545_V$33b3), (F64 *)sin_16545_V$33b1, (F64 *)cos_16545_V$33b2 ); [CALL_CONVENTION_UNIX_ABI] 333 1 sincosf{ic=SINCOSF}( (F32) *((F32*) a_16546_V$33b6), (F32 *)sin_16546_V$33b4, (F32 *)cos_16546_V$33b5 ); [CALL_CONVENTION_UNIX_ABI] fpoptimizer/grammar_data.cc(81): internal error: assertion failed: lower_dynamic_init_aggregate_constant: bad aggr kind (shared/cfe/edgcpfe/lower_init.c, line 6341) const ParamSpec_NumConstant <Value_t> plist_n_container<Value_t>::plist_n[20] = ^ compilation aborted for /home/akselos/software/libmeshsrc/contrib/fparser/fpoptimizer.cc (code 4) Any ideas of what the problem could be ? Also general comments on how to use intel tools properly with libmesh are welcome (intel compilers, intel mkl, intel mpi). Thanks, Sylvain 
From: Roy Stogner <roystgnr@ic...>  20150320 21:31:22

Copying to libmeshusers, since this frequentlyasked question's answer *should* have showed up in a quick Google search but didn't for me. On Fri, 20 Mar 2015, Sahai, Amal wrote: > Thank you so much for all your help. Quick question  I was going > through the shock fitting utility and it talks about how it works > for a serial mesh. Do you think I would have to make changes to the > routine it uses for doing ray tracing if I were to use it for a > FINS run in parallel? Sorry if my question is vague, I am still in > the process of figuring out what is being done there and wasn't able > to think of something more specific. In lieu of apology, whenever you find time you can help me write the inaugural section, "SerialMesh vs ParallelMesh", for the currentlynonexistent libMesh FAQ. I really should be smart enough to *save* these things when I write them out so I have a URL to point to in the future. Here's the simple explanation for libMesh users: Both SerialMesh and ParallelMesh work in parallel. SerialMesh works in parallel in a serialized fashion: every processor knows everything about the geometry and topology and partitioning and DoF indexing of every element and node in the mesh. ParallelMesh works in parallel in a distributed fashion: every processor knows everything about it's own "local" elements and about the "semilocal" elements which touch its own elements, but for elements just past the "semilocal" layer all each processor knows is that a "remote" element exists, and for elements beyond those the local processor has no information whatsoever. Obviously the latter data structure is much cheaper on memory and in large problems can be cheaper on the CPU, but it makes writing "nonlocal" algorithms like ray tracing harder since you can only do them in a tricky communicationsynchronized way. And another tip for you and other FINS users: Check out the TransitionBase code. That's another example of "node tracing", and IIRC it's still only SerialMesh compatible, but at least it's an example of the ray calculations being done on their home processors via manual communication.  Roy 
From: Roy Stogner <roystgnr@ic...>  20150320 04:42:13

On Thu, 19 Mar 2015, Zhang wrote: > Another question is :Does the DirichletBoundary class also uses > nodewise assignment implicitly? DirichletBoundary calculates constraints by first doing interpolation at nodes, then (in 2D and 3D) projection on edges holding nodal dof coefficients fixed, then (in 3D) projection on faces holding nodal and edge dof coefficients fixed. If you set incompatible constraints on different boundary sides that share the same node or edge (e.g. a liddriven cavity problem) then which constraint is applied depends on which constraint you added to the system first. The *earliest* added constraint takes precedence.  Roy 
From: Zhang <zyzhang@nu...>  20150319 04:54:21

Dear Vikram, The usage of class template in Feel++ is quite good. I didn't find the default quadrature order in case drivencavity yet. Please leave me some time. BTW, I checked the Ke & Fe of the element at exactly the left corner. The weak integral cannot make sure a good diagdominance at the boundary nodes. This is the direct reason to the abnormal nonzero wall velocity. In Ke, the diagonal nonzeros are at most one order of mag. higher than the largest offdiag nonzeros. Such a ratio does not change with the penalty range(1.e9 ~ 1.E10). Here, qrule(FIFTH),qface.n_points()=3, phi_vel_face.size()=9, fe_vel (SECOND),fe_pres(ONE)). So I guess there may be improper match between the quadrature order and the velocity order for this cavity problem. Am I right? please tell me if you can. Another question is :Does the DirichletBoundary class also uses nodewise assignment implicitly? Cheers, Zhenyu > 原始邮件 > 发件人: "Vikram Garg" <vikram.v.garg@...> > 发送时间: 20150319 08:47:43 (星期四) > 收件人: grandrabbit <zyzhang_nuaa@...> > 抄送: libmeshusers <libmeshusers@...> > 主题: Re: [Libmeshusers] Proper way of imposing Weak Dirichlet boundary condition > > Thats interesting. Do both the libMesh and Feel++ implementations of the > weak boundary condition use similar quadrature rules ? > > For the NavierStokes and advectiondiffusion equations, there is evidence > that weak enforcement of boundary conditions is superior to strong > enforcement, see 'Weak imposition of Dirichlet boundary conditions in fluid > mechanics' , Bazilevs, Hughes, 2007, > https://www.ices.utexas.edu/media/reports/2005/0525.pdf > > To my knowledge, we havent really used weak boundary conditions (Nitsche's > method) to enforce boundary conditions in libMesh. This is because these > methods are more problem independent. We generally use either the direct > nodewise assignment of the bc or the penalty method, and these havent led > to any major convergence issues so far. > > Thanks. > > On Wed, Mar 18, 2015 at 7:30 PM, grandrabbit <zyzhang_nuaa@...> wrote: > > > Dear Vikram, > > Yes. I do agree with your remarks on the geometrical singularity, indeed. > > Of course, refinement can remedy this to some extents. But from the > > view of a practical large problem, between this weak implementation of > > Dirichlet BC and the direct nodewise value assiging way, which is better? > > May the latter way introduce other qeustions, such as convergence? > > So this is my first question here. > > > > At the same time, I compared the result discussed above with the results > > from a Feel++ code with triangular meshes of roughly same mesh resolution > > (gmsh.hsize=0.05). Here further weak terms are added for > > symmetry of the Jacobian matrix, such as terms of > > \int_{Gamma_D} (\mu\nabla vqI)\cdot\vec{n}\cdot (uu_bc) > > Since Feel++ encapsulated all details of BC implementation quite well, I > > can't > > tell you how it works yet at the level of source codes. > > I inserted these terms in my libmesh code as well. But The comparison > > showed that the results did not match the accuracy of the Feel++ codes yet. > > If you like, I can send it to you later. > > Zhenyu > > > > > >  原始邮件  > > *发件人:* "Vikram Garg";<vikram.v.garg@...>; > > *发送时间:* 2015年3月19日(星期四) 凌晨0:59 > > *收件人:* "grandrabbit"<zyzhang_nuaa@...>; > > *抄送:* "libmeshusers"<libmeshusers@...>; > > *主题:* Re: [Libmeshusers] Proper way of imposing Weak Dirichlet boundary > > condition > > > > What happens if you refine the mesh ? The spurious oscillation you are > > seeing could be because of the singularity at the corners where the > > boundary condition is not well defined. > > > > On Tue, Mar 17, 2015 at 7:59 PM, grandrabbit <zyzhang_nuaa@...> wrote: > > > >> Thank you for your timely response. > >> > >> I reset h_elem to 1.0. No big difference. In fact in this case the mesh > >> is uniform, with 20 cells along one dimension, h_elem is a constant, say, > >> 1.25e2. Since I set penalty to 1.e+8, the actually penalty is 8.e9, > >> similar to intro_ex3. Therefore I think the true reason is located > >> elsewhere . BTW, is there any example for > >> TransientNonlinearImplicitSystem? When I start the coding, I referred to > >> systems_of_equations_ex2, which uses TransientLinearImplicitSystem and runs > >> well. I changed to TransientNonlinearImplicitSystem and rewrite the > >> residual function based on that. > >> > >> Anyway, I list the results below as reference. Some screenshots attached > >> as well. Further suggestion's welcome. > >> > >> case 1(good one): DirichletBoundary class used, no spurious velocity near > >> the wall > >> penalty 1.e8, residual 1.e16~1.e17, nonlinear convergence 7.28e4 > >> u in [0.196,1.0]; p in [110, +111] > >> DirichletBoundary_class.png > >> > >> case 2(problem now): self coded penalty integral (see laminar2dnl.cpp), > >> nonzero velocity on the corner wall > >> penalty 1.e8, residual 1.e9~1.e10, nonlinear convergence 7.286e4 > >> u in [0.203,1.07]; p in [195, +197] > >> h_elem: 1.25e2 > >> h_elem_0.0125.png > >> h_elem_0.0125_corner.png > >> > >> case 3(problem as well): self coded penalty integral (see > >> laminar2dnl.cpp), nonzero velocity on the corner wall > >> penalty 1.e8, residual 1.e11~1.e12, nonlinear convergence 7.28e4 > >> u in [0.203,1.07]; p in [195, +197] > >> h_elem: 1.0 > >> h_elem_1.0.png > >> h_elem_1.0_corner.png > >> > >> Cheers, > >> > >> Zhenyu > >> 原始邮件 > >> 发件人: "Vikram Garg" <vikram.v.garg@...> > >> 发送时间: 20150318 01:45:12 (星期三) > >> 收件人: Zhang <zyzhang@...> > >> 抄送: "libmeshusers@..." < > >> libmeshusers@...> > >> 主题: Re: [Libmeshusers] Proper way of imposing Weak Dirichlet boundary > >> condition > >> > >> Zhenyu, can you tell us what happens if you remove the h_elem term in > >> your penalty formulation ? Usually, we just set the penalty to a very > >> large number, and do not involve the mesh size in setting it. See example > >> 3: http://libmesh.github.io/examples/introduction_ex3.html > >> > >> > >> > >> On Tue, Mar 17, 2015 at 8:09 AM, Zhang <zyzhang@...> wrote: > >> Dear Libmesher, > >> > >> I tried to write an incompressible NaiverStokes solver based Libmesh. > >> Following features includes: > >> coupled pressurevelocity method > >> TransientNonlinearImplicitSystem used with > >> compute_jacobian(..) and compute_residual(..) run separatedly > >> Liddriven cavity case for initial tests > >> BDF1 time marching > >> QUAD4 element, > >> vel_order=2, presorder=1, both LAGRANGE elements > >> Pressure pinned at node_id=0 > >> > >> I tried DirichletBoundary class to impose the velocity BC, the code runs > >> OK. > >> > >> Now the problem is : > >> I coded myself the common penalty method for term > >> \int_{\Gamma_D}dt*\gamma/h*u*v, > >> and applied the forms in boundary integral part of jacobian function > >> Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[j][qp_face] > >> *phi_face_[i][qp_face] > >> and > >> Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[i][qp_face] *(uu_bc) > >> at counterparts of residual function > >> the results show the boundary velocity not fully applied ,esp. at the > >> side walls near > >> the two top corners. Actually there are nonzero u normal to the local > >> side wall. > >> So I wonder anybody here will kindly show me some hints. > >> > >> Furthermore, I tried applied the weak Dirichlet conditions by add futher > >> boundary integrals like > >> dt\int_{\partial\Omega} (\mu\nabla \vec{ u}+pII)\cdot > >> \vec{n}\cdot\vec{v} + > >> dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot > >> \vec{n}\cdot\vec{u} > >> at left and > >> dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot > >> \vec{n}\cdot\vec{u_bc} > >> > >> The case is even worse. > >> > >> I also attached the code, and I am looking forward to any suggestions. > >> Many thanks. > >> > >> Zhenyu > >> > >> > >>  > >> Dive into the World of Parallel Programming The Go Parallel Website, > >> sponsored > >> by Intel and developed in partnership with Slashdot Media, is your hub > >> for all > >> things parallel software development, from weekly thought leadership > >> blogs to > >> news, videos, case studies, tutorials and more. Take a look and join the > >> conversation now. http://goparallel.sourceforge.net/ > >> _______________________________________________ > >> Libmeshusers mailing list > >> Libmeshusers@... > >> https://lists.sourceforge.net/lists/listinfo/libmeshusers > >> > >> > >> > >> > >> > >>  > >> Vikram GargPostdoctoral Associate > >> Predictive Engineering and Computational Science (PECOS) > >> The University of Texas at Austin > >> http://web.mit.edu/vikramvg/www/ > >> > >> > >> > >> http://vikramvgarg.wordpress.com/ > >> > >> http://www.runforindia.org/runners/vikramg > >> > >>  > >> Dive into the World of Parallel Programming The Go Parallel Website, > >> sponsored > >> by Intel and developed in partnership with Slashdot Media, is your hub > >> for all > >> things parallel software development, from weekly thought leadership > >> blogs to > >> news, videos, case studies, tutorials and more. Take a look and join the > >> conversation now. http://goparallel.sourceforge.net/ > >> _______________________________________________ > >> Libmeshusers mailing list > >> Libmeshusers@... > >> https://lists.sourceforge.net/lists/listinfo/libmeshusers > >> > >> > > > > > >  > > Vikram Garg > > Postdoctoral Associate > > Predictive Engineering and Computational Science (PECOS) > > The University of Texas at Austin > > http://web.mit.edu/vikramvg/www/ > > > > http://vikramvgarg.wordpress.com/ > > http://www.runforindia.org/runners/vikramg > > > > > >  > Vikram Garg > Postdoctoral Associate > Predictive Engineering and Computational Science (PECOS) > The University of Texas at Austin > http://web.mit.edu/vikramvg/www/ > > http://vikramvgarg.wordpress.com/ > http://www.runforindia.org/runners/vikramg >  > Dive into the World of Parallel Programming The Go Parallel Website, sponsored > by Intel and developed in partnership with Slashdot Media, is your hub for all > things parallel software development, from weekly thought leadership blogs to > news, videos, case studies, tutorials and more. Take a look and join the > conversation now. http://goparallel.sourceforge.net/ > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers 
From: Vikram Garg <vikram.garg@gm...>  20150319 00:48:11

Thats interesting. Do both the libMesh and Feel++ implementations of the weak boundary condition use similar quadrature rules ? For the NavierStokes and advectiondiffusion equations, there is evidence that weak enforcement of boundary conditions is superior to strong enforcement, see 'Weak imposition of Dirichlet boundary conditions in fluid mechanics' , Bazilevs, Hughes, 2007, https://www.ices.utexas.edu/media/reports/2005/0525.pdf To my knowledge, we havent really used weak boundary conditions (Nitsche's method) to enforce boundary conditions in libMesh. This is because these methods are more problem independent. We generally use either the direct nodewise assignment of the bc or the penalty method, and these havent led to any major convergence issues so far. Thanks. On Wed, Mar 18, 2015 at 7:30 PM, grandrabbit <zyzhang_nuaa@...> wrote: > Dear Vikram, > Yes. I do agree with your remarks on the geometrical singularity, indeed. > Of course, refinement can remedy this to some extents. But from the > view of a practical large problem, between this weak implementation of > Dirichlet BC and the direct nodewise value assiging way, which is better? > May the latter way introduce other qeustions, such as convergence? > So this is my first question here. > > At the same time, I compared the result discussed above with the results > from a Feel++ code with triangular meshes of roughly same mesh resolution > (gmsh.hsize=0.05). Here further weak terms are added for > symmetry of the Jacobian matrix, such as terms of > \int_{Gamma_D} (\mu\nabla vqI)\cdot\vec{n}\cdot (uu_bc) > Since Feel++ encapsulated all details of BC implementation quite well, I > can't > tell you how it works yet at the level of source codes. > I inserted these terms in my libmesh code as well. But The comparison > showed that the results did not match the accuracy of the Feel++ codes yet. > If you like, I can send it to you later. > Zhenyu > > >  原始邮件  > *发件人:* "Vikram Garg";<vikram.v.garg@...>; > *发送时间:* 2015年3月19日(星期四) 凌晨0:59 > *收件人:* "grandrabbit"<zyzhang_nuaa@...>; > *抄送:* "libmeshusers"<libmeshusers@...>; > *主题:* Re: [Libmeshusers] Proper way of imposing Weak Dirichlet boundary > condition > > What happens if you refine the mesh ? The spurious oscillation you are > seeing could be because of the singularity at the corners where the > boundary condition is not well defined. > > On Tue, Mar 17, 2015 at 7:59 PM, grandrabbit <zyzhang_nuaa@...> wrote: > >> Thank you for your timely response. >> >> I reset h_elem to 1.0. No big difference. In fact in this case the mesh >> is uniform, with 20 cells along one dimension, h_elem is a constant, say, >> 1.25e2. Since I set penalty to 1.e+8, the actually penalty is 8.e9, >> similar to intro_ex3. Therefore I think the true reason is located >> elsewhere . BTW, is there any example for >> TransientNonlinearImplicitSystem? When I start the coding, I referred to >> systems_of_equations_ex2, which uses TransientLinearImplicitSystem and runs >> well. I changed to TransientNonlinearImplicitSystem and rewrite the >> residual function based on that. >> >> Anyway, I list the results below as reference. Some screenshots attached >> as well. Further suggestion's welcome. >> >> case 1(good one): DirichletBoundary class used, no spurious velocity near >> the wall >> penalty 1.e8, residual 1.e16~1.e17, nonlinear convergence 7.28e4 >> u in [0.196,1.0]; p in [110, +111] >> DirichletBoundary_class.png >> >> case 2(problem now): self coded penalty integral (see laminar2dnl.cpp), >> nonzero velocity on the corner wall >> penalty 1.e8, residual 1.e9~1.e10, nonlinear convergence 7.286e4 >> u in [0.203,1.07]; p in [195, +197] >> h_elem: 1.25e2 >> h_elem_0.0125.png >> h_elem_0.0125_corner.png >> >> case 3(problem as well): self coded penalty integral (see >> laminar2dnl.cpp), nonzero velocity on the corner wall >> penalty 1.e8, residual 1.e11~1.e12, nonlinear convergence 7.28e4 >> u in [0.203,1.07]; p in [195, +197] >> h_elem: 1.0 >> h_elem_1.0.png >> h_elem_1.0_corner.png >> >> Cheers, >> >> Zhenyu >> 原始邮件 >> 发件人: "Vikram Garg" <vikram.v.garg@...> >> 发送时间: 20150318 01:45:12 (星期三) >> 收件人: Zhang <zyzhang@...> >> 抄送: "libmeshusers@..." < >> libmeshusers@...> >> 主题: Re: [Libmeshusers] Proper way of imposing Weak Dirichlet boundary >> condition >> >> Zhenyu, can you tell us what happens if you remove the h_elem term in >> your penalty formulation ? Usually, we just set the penalty to a very >> large number, and do not involve the mesh size in setting it. See example >> 3: http://libmesh.github.io/examples/introduction_ex3.html >> >> >> >> On Tue, Mar 17, 2015 at 8:09 AM, Zhang <zyzhang@...> wrote: >> Dear Libmesher, >> >> I tried to write an incompressible NaiverStokes solver based Libmesh. >> Following features includes: >> coupled pressurevelocity method >> TransientNonlinearImplicitSystem used with >> compute_jacobian(..) and compute_residual(..) run separatedly >> Liddriven cavity case for initial tests >> BDF1 time marching >> QUAD4 element, >> vel_order=2, presorder=1, both LAGRANGE elements >> Pressure pinned at node_id=0 >> >> I tried DirichletBoundary class to impose the velocity BC, the code runs >> OK. >> >> Now the problem is : >> I coded myself the common penalty method for term >> \int_{\Gamma_D}dt*\gamma/h*u*v, >> and applied the forms in boundary integral part of jacobian function >> Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[j][qp_face] >> *phi_face_[i][qp_face] >> and >> Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[i][qp_face] *(uu_bc) >> at counterparts of residual function >> the results show the boundary velocity not fully applied ,esp. at the >> side walls near >> the two top corners. Actually there are nonzero u normal to the local >> side wall. >> So I wonder anybody here will kindly show me some hints. >> >> Furthermore, I tried applied the weak Dirichlet conditions by add futher >> boundary integrals like >> dt\int_{\partial\Omega} (\mu\nabla \vec{ u}+pII)\cdot >> \vec{n}\cdot\vec{v} + >> dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot >> \vec{n}\cdot\vec{u} >> at left and >> dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot >> \vec{n}\cdot\vec{u_bc} >> >> The case is even worse. >> >> I also attached the code, and I am looking forward to any suggestions. >> Many thanks. >> >> Zhenyu >> >> >>  >> Dive into the World of Parallel Programming The Go Parallel Website, >> sponsored >> by Intel and developed in partnership with Slashdot Media, is your hub >> for all >> things parallel software development, from weekly thought leadership >> blogs to >> news, videos, case studies, tutorials and more. Take a look and join the >> conversation now. http://goparallel.sourceforge.net/ >> _______________________________________________ >> Libmeshusers mailing list >> Libmeshusers@... >> https://lists.sourceforge.net/lists/listinfo/libmeshusers >> >> >> >> >> >>  >> Vikram GargPostdoctoral Associate >> Predictive Engineering and Computational Science (PECOS) >> The University of Texas at Austin >> http://web.mit.edu/vikramvg/www/ >> >> >> >> http://vikramvgarg.wordpress.com/ >> >> http://www.runforindia.org/runners/vikramg >> >>  >> Dive into the World of Parallel Programming The Go Parallel Website, >> sponsored >> by Intel and developed in partnership with Slashdot Media, is your hub >> for all >> things parallel software development, from weekly thought leadership >> blogs to >> news, videos, case studies, tutorials and more. Take a look and join the >> conversation now. http://goparallel.sourceforge.net/ >> _______________________________________________ >> Libmeshusers mailing list >> Libmeshusers@... >> https://lists.sourceforge.net/lists/listinfo/libmeshusers >> >> > > >  > Vikram Garg > Postdoctoral Associate > Predictive Engineering and Computational Science (PECOS) > The University of Texas at Austin > http://web.mit.edu/vikramvg/www/ > > http://vikramvgarg.wordpress.com/ > http://www.runforindia.org/runners/vikramg >  Vikram Garg Postdoctoral Associate Predictive Engineering and Computational Science (PECOS) The University of Texas at Austin http://web.mit.edu/vikramvg/www/ http://vikramvgarg.wordpress.com/ http://www.runforindia.org/runners/vikramg 
From: grandrabbit <zyzhang_nuaa@qq...>  20150319 00:30:58

Dear Vikram, Yes. I do agree with your remarks on the geometrical singularity, indeed. Of course, refinement can remedy this to some extents. But from the view of a practical large problem, between this weak implementation of Dirichlet BC and the direct nodewise value assiging way, which is better? May the latter way introduce other qeustions, such as convergence? So this is my first question here. At the same time, I compared the result discussed above with the results from a Feel++ code with triangular meshes of roughly same mesh resolution (gmsh.hsize=0.05). Here further weak terms are added for symmetry of the Jacobian matrix, such as terms of \int_{Gamma_D} (\mu\nabla vqI)\cdot\vec{n}\cdot (uu_bc) Since Feel++ encapsulated all details of BC implementation quite well, I can't tell you how it works yet at the level of source codes. I inserted these terms in my libmesh code as well. But The comparison showed that the results did not match the accuracy of the Feel++ codes yet. If you like, I can send it to you later. Zhenyu  原始邮件  发件人: "Vikram Garg";<vikram.v.garg@...>; 发送时间: 2015年3月19日(星期四) 凌晨0:59 收件人: "grandrabbit"<zyzhang_nuaa@...>; 抄送: "libmeshusers"<libmeshusers@...>; 主题: Re: [Libmeshusers] Proper way of imposing Weak Dirichlet boundary condition What happens if you refine the mesh ? The spurious oscillation you are seeing could be because of the singularity at the corners where the boundary condition is not well defined. On Tue, Mar 17, 2015 at 7:59 PM, grandrabbit <zyzhang_nuaa@...> wrote: Thank you for your timely response. I reset h_elem to 1.0. No big difference. In fact in this case the mesh is uniform, with 20 cells along one dimension, h_elem is a constant, say, 1.25e2. Since I set penalty to 1.e+8, the actually penalty is 8.e9, similar to intro_ex3. Therefore I think the true reason is located elsewhere . BTW, is there any example for TransientNonlinearImplicitSystem? When I start the coding, I referred to systems_of_equations_ex2, which uses TransientLinearImplicitSystem and runs well. I changed to TransientNonlinearImplicitSystem and rewrite the residual function based on that. Anyway, I list the results below as reference. Some screenshots attached as well. Further suggestion's welcome. case 1(good one): DirichletBoundary class used, no spurious velocity near the wall penalty 1.e8, residual 1.e16~1.e17, nonlinear convergence 7.28e4 u in [0.196,1.0]; p in [110, +111] DirichletBoundary_class.png case 2(problem now): self coded penalty integral (see laminar2dnl.cpp), nonzero velocity on the corner wall penalty 1.e8, residual 1.e9~1.e10, nonlinear convergence 7.286e4 u in [0.203,1.07]; p in [195, +197] h_elem: 1.25e2 h_elem_0.0125.png h_elem_0.0125_corner.png case 3(problem as well): self coded penalty integral (see laminar2dnl.cpp), nonzero velocity on the corner wall penalty 1.e8, residual 1.e11~1.e12, nonlinear convergence 7.28e4 u in [0.203,1.07]; p in [195, +197] h_elem: 1.0 h_elem_1.0.png h_elem_1.0_corner.png Cheers, Zhenyu 原始邮件 发件人: "Vikram Garg" <vikram.v.garg@...> 发送时间: 20150318 01:45:12 (星期三) 收件人: Zhang <zyzhang@...> 抄送: "libmeshusers@..." <libmeshusers@...> 主题: Re: [Libmeshusers] Proper way of imposing Weak Dirichlet boundary condition Zhenyu, can you tell us what happens if you remove the h_elem term in your penalty formulation ? Usually, we just set the penalty to a very large number, and do not involve the mesh size in setting it. See example 3: http://libmesh.github.io/examples/introduction_ex3.html On Tue, Mar 17, 2015 at 8:09 AM, Zhang <zyzhang@...> wrote: Dear Libmesher, I tried to write an incompressible NaiverStokes solver based Libmesh. Following features includes: coupled pressurevelocity method TransientNonlinearImplicitSystem used with compute_jacobian(..) and compute_residual(..) run separatedly Liddriven cavity case for initial tests BDF1 time marching QUAD4 element, vel_order=2, presorder=1, both LAGRANGE elements Pressure pinned at node_id=0 I tried DirichletBoundary class to impose the velocity BC, the code runs OK. Now the problem is : I coded myself the common penalty method for term \int_{\Gamma_D}dt*\gamma/h*u*v, and applied the forms in boundary integral part of jacobian function Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[j][qp_face] *phi_face_[i][qp_face] and Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[i][qp_face] *(uu_bc) at counterparts of residual function the results show the boundary velocity not fully applied ,esp. at the side walls near the two top corners. Actually there are nonzero u normal to the local side wall. So I wonder anybody here will kindly show me some hints. Furthermore, I tried applied the weak Dirichlet conditions by add futher boundary integrals like dt\int_{\partial\Omega} (\mu\nabla \vec{ u}+pII)\cdot \vec{n}\cdot\vec{v} + dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot \vec{n}\cdot\vec{u} at left and dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot \vec{n}\cdot\vec{u_bc} The case is even worse. I also attached the code, and I am looking forward to any suggestions. Many thanks. Zhenyu  Dive into the World of Parallel Programming The Go Parallel Website, sponsored by Intel and developed in partnership with Slashdot Media, is your hub for all things parallel software development, from weekly thought leadership blogs to news, videos, case studies, tutorials and more. Take a look and join the conversation now. http://goparallel.sourceforge.net/ _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers  Vikram GargPostdoctoral Associate Predictive Engineering and Computational Science (PECOS) The University of Texas at Austin http://web.mit.edu/vikramvg/www/ http://vikramvgarg.wordpress.com/ http://www.runforindia.org/runners/vikramg  Dive into the World of Parallel Programming The Go Parallel Website, sponsored by Intel and developed in partnership with Slashdot Media, is your hub for all things parallel software development, from weekly thought leadership blogs to news, videos, case studies, tutorials and more. Take a look and join the conversation now. http://goparallel.sourceforge.net/ _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers  Vikram Garg Postdoctoral Associate Predictive Engineering and Computational Science (PECOS) The University of Texas at Austin http://web.mit.edu/vikramvg/www/ http://vikramvgarg.wordpress.com/ http://www.runforindia.org/runners/vikramg 
From: Vikram Garg <vikram.garg@gm...>  20150318 16:59:29

What happens if you refine the mesh ? The spurious oscillation you are seeing could be because of the singularity at the corners where the boundary condition is not well defined. On Tue, Mar 17, 2015 at 7:59 PM, grandrabbit <zyzhang_nuaa@...> wrote: > Thank you for your timely response. > > I reset h_elem to 1.0. No big difference. In fact in this case the mesh is > uniform, with 20 cells along one dimension, h_elem is a constant, say, > 1.25e2. Since I set penalty to 1.e+8, the actually penalty is 8.e9, > similar to intro_ex3. Therefore I think the true reason is located > elsewhere . BTW, is there any example for > TransientNonlinearImplicitSystem? When I start the coding, I referred to > systems_of_equations_ex2, which uses TransientLinearImplicitSystem and runs > well. I changed to TransientNonlinearImplicitSystem and rewrite the > residual function based on that. > > Anyway, I list the results below as reference. Some screenshots attached > as well. Further suggestion's welcome. > > case 1(good one): DirichletBoundary class used, no spurious velocity near > the wall > penalty 1.e8, residual 1.e16~1.e17, nonlinear convergence 7.28e4 > u in [0.196,1.0]; p in [110, +111] > DirichletBoundary_class.png > > case 2(problem now): self coded penalty integral (see laminar2dnl.cpp), > nonzero velocity on the corner wall > penalty 1.e8, residual 1.e9~1.e10, nonlinear convergence 7.286e4 > u in [0.203,1.07]; p in [195, +197] > h_elem: 1.25e2 > h_elem_0.0125.png > h_elem_0.0125_corner.png > > case 3(problem as well): self coded penalty integral (see > laminar2dnl.cpp), nonzero velocity on the corner wall > penalty 1.e8, residual 1.e11~1.e12, nonlinear convergence 7.28e4 > u in [0.203,1.07]; p in [195, +197] > h_elem: 1.0 > h_elem_1.0.png > h_elem_1.0_corner.png > > Cheers, > > Zhenyu > 原始邮件 > 发件人: "Vikram Garg" <vikram.v.garg@...> > 发送时间: 20150318 01:45:12 (星期三) > 收件人: Zhang <zyzhang@...> > 抄送: "libmeshusers@..." < > libmeshusers@...> > 主题: Re: [Libmeshusers] Proper way of imposing Weak Dirichlet boundary > condition > > Zhenyu, can you tell us what happens if you remove the h_elem term in > your penalty formulation ? Usually, we just set the penalty to a very > large number, and do not involve the mesh size in setting it. See example > 3: http://libmesh.github.io/examples/introduction_ex3.html > > > > On Tue, Mar 17, 2015 at 8:09 AM, Zhang <zyzhang@...> wrote: > Dear Libmesher, > > I tried to write an incompressible NaiverStokes solver based Libmesh. > Following features includes: > coupled pressurevelocity method > TransientNonlinearImplicitSystem used with > compute_jacobian(..) and compute_residual(..) run separatedly > Liddriven cavity case for initial tests > BDF1 time marching > QUAD4 element, > vel_order=2, presorder=1, both LAGRANGE elements > Pressure pinned at node_id=0 > > I tried DirichletBoundary class to impose the velocity BC, the code runs > OK. > > Now the problem is : > I coded myself the common penalty method for term > \int_{\Gamma_D}dt*\gamma/h*u*v, > and applied the forms in boundary integral part of jacobian function > Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[j][qp_face] > *phi_face_[i][qp_face] > and > Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[i][qp_face] *(uu_bc) > at counterparts of residual function > the results show the boundary velocity not fully applied ,esp. at the > side walls near > the two top corners. Actually there are nonzero u normal to the local > side wall. > So I wonder anybody here will kindly show me some hints. > > Furthermore, I tried applied the weak Dirichlet conditions by add futher > boundary integrals like > dt\int_{\partial\Omega} (\mu\nabla \vec{ u}+pII)\cdot > \vec{n}\cdot\vec{v} + > dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot \vec{n}\cdot\vec{u} > at left and > dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot > \vec{n}\cdot\vec{u_bc} > > The case is even worse. > > I also attached the code, and I am looking forward to any suggestions. > Many thanks. > > Zhenyu > > >  > Dive into the World of Parallel Programming The Go Parallel Website, > sponsored > by Intel and developed in partnership with Slashdot Media, is your hub > for all > things parallel software development, from weekly thought leadership > blogs to > news, videos, case studies, tutorials and more. Take a look and join the > conversation now. http://goparallel.sourceforge.net/ > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > > > > > >  > Vikram GargPostdoctoral Associate > Predictive Engineering and Computational Science (PECOS) > The University of Texas at Austin > http://web.mit.edu/vikramvg/www/ > > > > http://vikramvgarg.wordpress.com/ > > http://www.runforindia.org/runners/vikramg > >  > Dive into the World of Parallel Programming The Go Parallel Website, > sponsored > by Intel and developed in partnership with Slashdot Media, is your hub for > all > things parallel software development, from weekly thought leadership blogs > to > news, videos, case studies, tutorials and more. Take a look and join the > conversation now. http://goparallel.sourceforge.net/ > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > >  Vikram Garg Postdoctoral Associate Predictive Engineering and Computational Science (PECOS) The University of Texas at Austin http://web.mit.edu/vikramvg/www/ http://vikramvgarg.wordpress.com/ http://www.runforindia.org/runners/vikramg 
From: grandrabbit <zyzhang_nuaa@qq...>  20150318 00:59:33

Thank you for your timely response. I reset h_elem to 1.0. No big difference. In fact in this case the mesh is uniform, with 20 cells along one dimension, h_elem is a constant, say, 1.25e2. Since I set penalty to 1.e+8, the actually penalty is 8.e9, similar to intro_ex3. Therefore I think the true reason is located elsewhere . BTW, is there any example for TransientNonlinearImplicitSystem? When I start the coding, I referred to systems_of_equations_ex2, which uses TransientLinearImplicitSystem and runs well. I changed to TransientNonlinearImplicitSystem and rewrite the residual function based on that. Anyway, I list the results below as reference. Some screenshots attached as well. Further suggestion's welcome. case 1(good one): DirichletBoundary class used, no spurious velocity near the wall penalty 1.e8, residual 1.e16~1.e17, nonlinear convergence 7.28e4 u in [0.196,1.0]; p in [110, +111] DirichletBoundary_class.png case 2(problem now): self coded penalty integral (see laminar2dnl.cpp), nonzero velocity on the corner wall penalty 1.e8, residual 1.e9~1.e10, nonlinear convergence 7.286e4 u in [0.203,1.07]; p in [195, +197] h_elem: 1.25e2 h_elem_0.0125.png h_elem_0.0125_corner.png case 3(problem as well): self coded penalty integral (see laminar2dnl.cpp), nonzero velocity on the corner wall penalty 1.e8, residual 1.e11~1.e12, nonlinear convergence 7.28e4 u in [0.203,1.07]; p in [195, +197] h_elem: 1.0 h_elem_1.0.png h_elem_1.0_corner.png Cheers, Zhenyu 原始邮件 发件人: "Vikram Garg" <vikram.v.garg@...> 发送时间: 20150318 01:45:12 (星期三) 收件人: Zhang <zyzhang@...> 抄送: "libmeshusers@..." <libmeshusers@...> 主题: Re: [Libmeshusers] Proper way of imposing Weak Dirichlet boundary condition Zhenyu, can you tell us what happens if you remove the h_elem term in your penalty formulation ? Usually, we just set the penalty to a very large number, and do not involve the mesh size in setting it. See example 3: http://libmesh.github.io/examples/introduction_ex3.html On Tue, Mar 17, 2015 at 8:09 AM, Zhang <zyzhang@...> wrote: Dear Libmesher, I tried to write an incompressible NaiverStokes solver based Libmesh. Following features includes: coupled pressurevelocity method TransientNonlinearImplicitSystem used with compute_jacobian(..) and compute_residual(..) run separatedly Liddriven cavity case for initial tests BDF1 time marching QUAD4 element, vel_order=2, presorder=1, both LAGRANGE elements Pressure pinned at node_id=0 I tried DirichletBoundary class to impose the velocity BC, the code runs OK. Now the problem is : I coded myself the common penalty method for term \int_{\Gamma_D}dt*\gamma/h*u*v, and applied the forms in boundary integral part of jacobian function Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[j][qp_face] *phi_face_[i][qp_face] and Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[i][qp_face] *(uu_bc) at counterparts of residual function the results show the boundary velocity not fully applied ,esp. at the side walls near the two top corners. Actually there are nonzero u normal to the local side wall. So I wonder anybody here will kindly show me some hints. Furthermore, I tried applied the weak Dirichlet conditions by add futher boundary integrals like dt\int_{\partial\Omega} (\mu\nabla \vec{ u}+pII)\cdot \vec{n}\cdot\vec{v} + dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot \vec{n}\cdot\vec{u} at left and dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot \vec{n}\cdot\vec{u_bc} The case is even worse. I also attached the code, and I am looking forward to any suggestions. Many thanks. Zhenyu  Dive into the World of Parallel Programming The Go Parallel Website, sponsored by Intel and developed in partnership with Slashdot Media, is your hub for all things parallel software development, from weekly thought leadership blogs to news, videos, case studies, tutorials and more. Take a look and join the conversation now. http://goparallel.sourceforge.net/ _______________________________________________ Libmeshusers mailing list Libmeshusers@... https://lists.sourceforge.net/lists/listinfo/libmeshusers  Vikram GargPostdoctoral Associate Predictive Engineering and Computational Science (PECOS) The University of Texas at Austin http://web.mit.edu/vikramvg/www/ http://vikramvgarg.wordpress.com/ http://www.runforindia.org/runners/vikramg 
From: Vikram Garg <vikram.garg@gm...>  20150317 17:45:41

Zhenyu, can you tell us what happens if you remove the h_elem term in your penalty formulation ? Usually, we just set the penalty to a very large number, and do not involve the mesh size in setting it. See example 3: http://libmesh.github.io/examples/introduction_ex3.html On Tue, Mar 17, 2015 at 8:09 AM, Zhang <zyzhang@...> wrote: > Dear Libmesher, > > I tried to write an incompressible NaiverStokes solver based Libmesh. > Following features includes: > coupled pressurevelocity method > TransientNonlinearImplicitSystem used with > compute_jacobian(..) and compute_residual(..) run separatedly > Liddriven cavity case for initial tests > BDF1 time marching > QUAD4 element, > vel_order=2, presorder=1, both LAGRANGE elements > Pressure pinned at node_id=0 > > I tried DirichletBoundary class to impose the velocity BC, the code runs > OK. > > Now the problem is : > I coded myself the common penalty method for term > \int_{\Gamma_D}dt*\gamma/h*u*v, > and applied the forms in boundary integral part of jacobian function > Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[j][qp_face] > *phi_face_[i][qp_face] > and > Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[i][qp_face] *(uu_bc) > at counterparts of residual function > the results show the boundary velocity not fully applied ,esp. at the side > walls near > the two top corners. Actually there are nonzero u normal to the local side > wall. > So I wonder anybody here will kindly show me some hints. > > Furthermore, I tried applied the weak Dirichlet conditions by add futher > boundary integrals like > dt\int_{\partial\Omega} (\mu\nabla \vec{ u}+pII)\cdot \vec{n}\cdot\vec{v} > + > dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot \vec{n}\cdot\vec{u} > at left and > dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot > \vec{n}\cdot\vec{u_bc} > > The case is even worse. > > I also attached the code, and I am looking forward to any suggestions. > Many thanks. > > Zhenyu > > >  > Dive into the World of Parallel Programming The Go Parallel Website, > sponsored > by Intel and developed in partnership with Slashdot Media, is your hub for > all > things parallel software development, from weekly thought leadership blogs > to > news, videos, case studies, tutorials and more. Take a look and join the > conversation now. http://goparallel.sourceforge.net/ > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > >  Vikram Garg Postdoctoral Associate Predictive Engineering and Computational Science (PECOS) The University of Texas at Austin http://web.mit.edu/vikramvg/www/ http://vikramvgarg.wordpress.com/ http://www.runforindia.org/runners/vikramg 
From: Zhang <zyzhang@nu...>  20150317 13:39:02

Dear Libmesher, I tried to write an incompressible NaiverStokes solver based Libmesh. Following features includes: coupled pressurevelocity method TransientNonlinearImplicitSystem used with compute_jacobian(..) and compute_residual(..) run separatedly Liddriven cavity case for initial tests BDF1 time marching QUAD4 element, vel_order=2, presorder=1, both LAGRANGE elements Pressure pinned at node_id=0 I tried DirichletBoundary class to impose the velocity BC, the code runs OK. Now the problem is : I coded myself the common penalty method for term \int_{\Gamma_D}dt*\gamma/h*u*v, and applied the forms in boundary integral part of jacobian function Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[j][qp_face] *phi_face_[i][qp_face] and Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[i][qp_face] *(uu_bc) at counterparts of residual function the results show the boundary velocity not fully applied ,esp. at the side walls near the two top corners. Actually there are nonzero u normal to the local side wall. So I wonder anybody here will kindly show me some hints. Furthermore, I tried applied the weak Dirichlet conditions by add futher boundary integrals like dt\int_{\partial\Omega} (\mu\nabla \vec{ u}+pII)\cdot \vec{n}\cdot\vec{v} + dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot \vec{n}\cdot\vec{u} at left and dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot \vec{n}\cdot\vec{u_bc} The case is even worse. I also attached the code, and I am looking forward to any suggestions. Many thanks. Zhenyu 
From: grandrabbit <zyzhang_nuaa@qq...>  20150317 13:35:45

Dear Libmesher, I tried to write an incompressible NaiverStokes solver based Libmesh. Following features includes: coupled pressurevelocity method TransientNonlinearImplicitSystem used with compute_jacobian(..) and compute_residual(..) run separatedly Liddriven cavity case for initial tests BDF1 time marching QUAD4 element, vel_order=2, presorder=1, both LAGRANGE elements Pressure pinned at node_id=0 I tried DirichletBoundary class to impose the velocity BC, the code runs OK. Now the problem is : I coded myself the common penalty method for term \int_{\Gamma_D}dt*\gamma/h*u*v, and applied the forms in boundary integral part of jacobian function Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[j][qp_face] *phi_face_[i][qp_face] and Jxw_face[qp_face]*dt*penalty/h_elem*phi_face_[i][qp_face] *(uu_bc) at counterparts of residual function the results show the boundary velocity not fully applied ,esp. at the side walls near the two top corners. Actually there are nonzero u normal to the local side wall. So I wonder anybody here will kindly show me some hints. Furthermore, I tried applied the weak Dirichlet conditions by add futher boundary integrals like dt\int_{\partial\Omega} (\mu\nabla \vec{ u}+pII)\cdot \vec{n}\cdot\vec{v} + dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot \vec{n}\cdot\vec{u} at left and dt\int_{\partial\Omega} (\mu\nabla \vec{ v}+qII)\cdot \vec{n}\cdot\vec{u_bc} at right hand side But the case is even worse. I also attached the code here, and I am looking forward to any suggestions. Many thanks. Zhenyu 
From: John Peterson <jwpeterson@gm...>  20150313 18:35:36

On Fri, Mar 13, 2015 at 12:28 PM, Manav Bhatia <bhatiamanav@...> wrote: > Thanks for the correction, John. I had meant to say "float". > > At my end, I have temporarily changed all "floats" to "double" in the > metis/parmetis files in libmesh, but this is a temporary fix. > > It would be better to source it from the respective headers. But that > raises another point: as of now, it seems like the metis.h in contrib has a > hardcoded value of real width to be 32. I understand that this saves memory > space, but would there be a motivation to keep this scalar type to be > consistent with the rest of the library? > The floating point stuff is bad... converting double to float certainly has loss of precision issues, but I'm also worried about the case when PETSc and libmesh are configured with with64bitindices and withdofidbytes=8, respectively, and we continue using int when passing data through the libmesh/MetisParmetis interfaces.  John 
From: Manav Bhatia <bhatiamanav@gm...>  20150313 18:28:43

Thanks for the correction, John. I had meant to say "float". At my end, I have temporarily changed all "floats" to "double" in the metis/parmetis files in libmesh, but this is a temporary fix. It would be better to source it from the respective headers. But that raises another point: as of now, it seems like the metis.h in contrib has a hardcoded value of real width to be 32. I understand that this saves memory space, but would there be a motivation to keep this scalar type to be consistent with the rest of the library? Manav On Fri, Mar 13, 2015 at 1:18 PM, John Peterson <jwpeterson@...> wrote: > > > On Fri, Mar 13, 2015 at 11:54 AM, Manav Bhatia <bhatiamanav@...> > wrote: > >> Correction: I had meant to write "real as the scalar type" (instead of >> "scalar as the real type") in my previous message. >> > > Hmm... I don't see a single instance of Real in either metis_partitioner.h > or parmetis_partitioner.h. Perhaps you are referring to: > > std::vector<float> _tpwgts; > std::vector<float> _ubvec; > > which use float? I agree that should probably be changed, along with > using Metis' idx_t for indexing instead of hardcoding int. > >  > John > 
From: John Peterson <jwpeterson@gm...>  20150313 18:19:25

On Fri, Mar 13, 2015 at 11:54 AM, Manav Bhatia <bhatiamanav@...> wrote: > Correction: I had meant to write "real as the scalar type" (instead of > "scalar as the real type") in my previous message. > Hmm... I don't see a single instance of Real in either metis_partitioner.h or parmetis_partitioner.h. Perhaps you are referring to: std::vector<float> _tpwgts; std::vector<float> _ubvec; which use float? I agree that should probably be changed, along with using Metis' idx_t for indexing instead of hardcoding int.  John 
From: Manav Bhatia <bhatiamanav@gm...>  20150313 17:54:08

Correction: I had meant to write "real as the scalar type" (instead of "scalar as the real type") in my previous message. Manav On Fri, Mar 13, 2015 at 12:44 PM, Manav Bhatia <bhatiamanav@...> wrote: > Is there a motivation to keep "scalar" as the real type in > parmetis_partitioner.h? > > Why not change this to the real type in metis.h? > > Manav > > > On Fri, Mar 13, 2015 at 12:27 PM, Manav Bhatia <bhatiamanav@...> > wrote: > >> not a problem... >> enjoy your conference and have a safe trip! >> >> Manav >> >> >> On Fri, Mar 13, 2015 at 12:24 PM, Paul T. Bauman <ptbauman@...> >> wrote: >> >>> Thanks for the report. We'll get it sorted soon. I'm about to leave for >>> SIAM CSE, as are other developers I anticipate so it may be a few days. >>> >>> On Fri, Mar 13, 2015 at 1:19 PM, Manav Bhatia <bhatiamanav@...> >>> wrote: >>> >>>> Found it!!! >>>> >>>> So, I did specify "PETSc" as the source of metis during compilation of >>>> libmesh (as I have been doing successfully for a while on my mac). However, >>>> for some reason, the compile options included the libmesh/contrib/metis and >>>> libmesh/contrib/parmetis as search paths. >>>> >>>> As a result, the compiler was using metis.h from libmesh/contrib, but >>>> using libmetis.so and libparmetis.so library from petsc. The >>>> libmesh/contrib has a real width of 32, and uses "float" as scalar in >>>> parmetis_partitioner, while the petsc metis library was build with real >>>> width type of 64. >>>> >>>> I do not know why contrib stuff is still being used, but this mismatch >>>> was causing the error. >>>> >>>> On Fri, Mar 13, 2015 at 6:24 AM, Paul T. Bauman <ptbauman@...> >>>> wrote: >>>> >>>>> On Fri, Mar 13, 2015 at 12:52 AM, Manav Bhatia <bhatiamanav@...> >>>>> wrote: >>>>> >>>>>> Parmetis is quitting on me with segmentation faults. The code >>>>>> works just fine on other machines (both mac and some big intel based >>>>>> clusters). Before I spend more time debugging it, I wanted to get a quick >>>>>> check from others if they have had issues if parmetis. >>>>>> >>>>> >>>>> Make sure you use the configure option withmetis=PETSc >>>>> >>>>> The libMesh and PETSc ParMetis installations can (and often do) >>>>> collide and that option will tell libMesh to use PETSc's. >>>>> >>>> >>>> >>> >> > 
From: Manav Bhatia <bhatiamanav@gm...>  20150313 17:44:39

Is there a motivation to keep "scalar" as the real type in parmetis_partitioner.h? Why not change this to the real type in metis.h? Manav On Fri, Mar 13, 2015 at 12:27 PM, Manav Bhatia <bhatiamanav@...> wrote: > not a problem... > enjoy your conference and have a safe trip! > > Manav > > > On Fri, Mar 13, 2015 at 12:24 PM, Paul T. Bauman <ptbauman@...> > wrote: > >> Thanks for the report. We'll get it sorted soon. I'm about to leave for >> SIAM CSE, as are other developers I anticipate so it may be a few days. >> >> On Fri, Mar 13, 2015 at 1:19 PM, Manav Bhatia <bhatiamanav@...> >> wrote: >> >>> Found it!!! >>> >>> So, I did specify "PETSc" as the source of metis during compilation of >>> libmesh (as I have been doing successfully for a while on my mac). However, >>> for some reason, the compile options included the libmesh/contrib/metis and >>> libmesh/contrib/parmetis as search paths. >>> >>> As a result, the compiler was using metis.h from libmesh/contrib, but >>> using libmetis.so and libparmetis.so library from petsc. The >>> libmesh/contrib has a real width of 32, and uses "float" as scalar in >>> parmetis_partitioner, while the petsc metis library was build with real >>> width type of 64. >>> >>> I do not know why contrib stuff is still being used, but this mismatch >>> was causing the error. >>> >>> On Fri, Mar 13, 2015 at 6:24 AM, Paul T. Bauman <ptbauman@...> >>> wrote: >>> >>>> On Fri, Mar 13, 2015 at 12:52 AM, Manav Bhatia <bhatiamanav@...> >>>> wrote: >>>> >>>>> Parmetis is quitting on me with segmentation faults. The code works >>>>> just fine on other machines (both mac and some big intel based clusters). >>>>> Before I spend more time debugging it, I wanted to get a quick check from >>>>> others if they have had issues if parmetis. >>>>> >>>> >>>> Make sure you use the configure option withmetis=PETSc >>>> >>>> The libMesh and PETSc ParMetis installations can (and often do) collide >>>> and that option will tell libMesh to use PETSc's. >>>> >>> >>> >> > 
From: Manav Bhatia <bhatiamanav@gm...>  20150313 17:27:43

not a problem... enjoy your conference and have a safe trip! Manav On Fri, Mar 13, 2015 at 12:24 PM, Paul T. Bauman <ptbauman@...> wrote: > Thanks for the report. We'll get it sorted soon. I'm about to leave for > SIAM CSE, as are other developers I anticipate so it may be a few days. > > On Fri, Mar 13, 2015 at 1:19 PM, Manav Bhatia <bhatiamanav@...> > wrote: > >> Found it!!! >> >> So, I did specify "PETSc" as the source of metis during compilation of >> libmesh (as I have been doing successfully for a while on my mac). However, >> for some reason, the compile options included the libmesh/contrib/metis and >> libmesh/contrib/parmetis as search paths. >> >> As a result, the compiler was using metis.h from libmesh/contrib, but >> using libmetis.so and libparmetis.so library from petsc. The >> libmesh/contrib has a real width of 32, and uses "float" as scalar in >> parmetis_partitioner, while the petsc metis library was build with real >> width type of 64. >> >> I do not know why contrib stuff is still being used, but this mismatch >> was causing the error. >> >> On Fri, Mar 13, 2015 at 6:24 AM, Paul T. Bauman <ptbauman@...> >> wrote: >> >>> On Fri, Mar 13, 2015 at 12:52 AM, Manav Bhatia <bhatiamanav@...> >>> wrote: >>> >>>> Parmetis is quitting on me with segmentation faults. The code works >>>> just fine on other machines (both mac and some big intel based clusters). >>>> Before I spend more time debugging it, I wanted to get a quick check from >>>> others if they have had issues if parmetis. >>>> >>> >>> Make sure you use the configure option withmetis=PETSc >>> >>> The libMesh and PETSc ParMetis installations can (and often do) collide >>> and that option will tell libMesh to use PETSc's. >>> >> >> > 
From: Paul T. Bauman <ptbauman@gm...>  20150313 17:24:32

Thanks for the report. We'll get it sorted soon. I'm about to leave for SIAM CSE, as are other developers I anticipate so it may be a few days. On Fri, Mar 13, 2015 at 1:19 PM, Manav Bhatia <bhatiamanav@...> wrote: > Found it!!! > > So, I did specify "PETSc" as the source of metis during compilation of > libmesh (as I have been doing successfully for a while on my mac). However, > for some reason, the compile options included the libmesh/contrib/metis and > libmesh/contrib/parmetis as search paths. > > As a result, the compiler was using metis.h from libmesh/contrib, but > using libmetis.so and libparmetis.so library from petsc. The > libmesh/contrib has a real width of 32, and uses "float" as scalar in > parmetis_partitioner, while the petsc metis library was build with real > width type of 64. > > I do not know why contrib stuff is still being used, but this mismatch was > causing the error. > > On Fri, Mar 13, 2015 at 6:24 AM, Paul T. Bauman <ptbauman@...> > wrote: > >> On Fri, Mar 13, 2015 at 12:52 AM, Manav Bhatia <bhatiamanav@...> >> wrote: >> >>> Parmetis is quitting on me with segmentation faults. The code works >>> just fine on other machines (both mac and some big intel based clusters). >>> Before I spend more time debugging it, I wanted to get a quick check from >>> others if they have had issues if parmetis. >>> >> >> Make sure you use the configure option withmetis=PETSc >> >> The libMesh and PETSc ParMetis installations can (and often do) collide >> and that option will tell libMesh to use PETSc's. >> > > 
From: Manav Bhatia <bhatiamanav@gm...>  20150313 17:19:44

Found it!!! So, I did specify "PETSc" as the source of metis during compilation of libmesh (as I have been doing successfully for a while on my mac). However, for some reason, the compile options included the libmesh/contrib/metis and libmesh/contrib/parmetis as search paths. As a result, the compiler was using metis.h from libmesh/contrib, but using libmetis.so and libparmetis.so library from petsc. The libmesh/contrib has a real width of 32, and uses "float" as scalar in parmetis_partitioner, while the petsc metis library was build with real width type of 64. I do not know why contrib stuff is still being used, but this mismatch was causing the error. On Fri, Mar 13, 2015 at 6:24 AM, Paul T. Bauman <ptbauman@...> wrote: > On Fri, Mar 13, 2015 at 12:52 AM, Manav Bhatia <bhatiamanav@...> > wrote: > >> Parmetis is quitting on me with segmentation faults. The code works >> just fine on other machines (both mac and some big intel based clusters). >> Before I spend more time debugging it, I wanted to get a quick check from >> others if they have had issues if parmetis. >> > > Make sure you use the configure option withmetis=PETSc > > The libMesh and PETSc ParMetis installations can (and often do) collide > and that option will tell libMesh to use PETSc's. > 
From: Dmitry Karpeyev <karpeev@mc...>  20150313 17:06:35

Apologies in advance for multipleposting. There is a postdoctoral position available at the University of Chicago that may be of interest to the members of this list. The focus is on the development of new and existing PETSc solvers required by a range of applications. Specifically, we are interested in the development of effective solvers and preconditioners for variational inequalities with general constraints, improvements to linear/nonlinear multigrid preconditioners, nonlinear splitting and domaindecomposition methods. Relevant applications include nonlinear elasticity with frictional contact, stronglycoupled multiphysics models, problems with moving meshes, mortar finiteelement formulations, etc. Knowledge of PETSc, MPI programming and numerical algorithms is required. Familiarity with modern highperformance hardware as well as with FEM libraries libMesh and MOOSE is a big plus. Please forward this information to your colleagues that might be interested in this opportunity. Feel free to contact me for more details. Thanks! Dmitry. Dmitry Karpeyev Computation Institute University of Chicago and Argonne National Laboratory 
From: Manav Bhatia <bhatiamanav@gm...>  20150313 16:52:32

Yes, I have been doing that. U am getting an error from Parmetis saying The sum of tpwgts for constraint #0 is not 1.0 I searched Parmetis' website and did not find anything conclusive to point me in a direction. There was someone who did a reinstall of mpi to get around this, but that is hardly a solution in my case. Manav > On Mar 13, 2015, at 6:24 AM, Paul T. Bauman <ptbauman@...> wrote: > >> On Fri, Mar 13, 2015 at 12:52 AM, Manav Bhatia <bhatiamanav@...> wrote: >> Parmetis is quitting on me with segmentation faults. The code works just fine on other machines (both mac and some big intel based clusters). Before I spend more time debugging it, I wanted to get a quick check from others if they have had issues if parmetis. > > Make sure you use the configure option withmetis=PETSc > > The libMesh and PETSc ParMetis installations can (and often do) collide and that option will tell libMesh to use PETSc's. 
From: Paul T. Bauman <ptbauman@gm...>  20150313 11:24:41

On Fri, Mar 13, 2015 at 12:52 AM, Manav Bhatia <bhatiamanav@...> wrote: > Parmetis is quitting on me with segmentation faults. The code works > just fine on other machines (both mac and some big intel based clusters). > Before I spend more time debugging it, I wanted to get a quick check from > others if they have had issues if parmetis. > Make sure you use the configure option withmetis=PETSc The libMesh and PETSc ParMetis installations can (and often do) collide and that option will tell libMesh to use PETSc's. 
From: Manav Bhatia <bhatiamanav@gm...>  20150313 04:52:11

Hi, I am now running my code linked to libMesh that uses parmetis from petsc version 3.5.3. The specific version of parmetis is 4.0.2. My compiler is gcc 4.4.7 with openmpi1.8.3. Parmetis is quitting on me with segmentation faults. The code works just fine on other machines (both mac and some big intel based clusters). Before I spend more time debugging it, I wanted to get a quick check from others if they have had issues if parmetis. Thanks, Manav 
From: Lei Zhao <leizhao987@gm...>  20150312 17:23:34

On Thu, Mar 12, 2015 at 11:50 AM, Cody Permann <codypermann@...> wrote: > Try to stay out of parallel_implementation.h. Generally you don't need to > go down that far unless you are adding new type specializations. The > signatures you need for most communication patterns are in parallel.h: > > /** > * Blockingsend to one processor with datadefined type. > */ > template <typename T> > void send (const unsigned int dest_processor_id, > T &buf, > const MessageTag &tag=no_tag) const; > > /** > * Blockingreceive from one processor with datadefined type. > */ > template <typename T> > Status receive (const unsigned int dest_processor_id, > T &buf, > const MessageTag &tag=any_tag) const; > > > So sending/recieving looks more like this: > _communicator.send(rank, some_vector); > _communicator.receive(rank, some_vector): > I tried and it works. Thanks! I think I should add some "tag" for safety, right? But when I try the following, the vector was sent, but never received. Why? _communicator.send(rank, send_vector, 101); _communicator.send(rank, receive_vector, 101); > BTW  When you are using MOOSE, you'll have a builtin Communicator in > every object. It's simply "_comm". > You mean simply use _comm instead of _communicator? What are the differences between them? Thank you, Lei > Cody > > > > > On Thu, Mar 12, 2015 at 10:40 AM Lei Zhao <leizhao987@...> wrote: > >> Hi All, >> >> I am working on some parallelism with Communicator. But I don't understand >> exactly the syntax of _communicator.send and _communicator.receive. This >> is >> the definition of send for sending a vector in parallel_implementation.h. >> >> *template <typename T>* >> *inline void Communicator::send (const unsigned int dest_processor_id,* >> * std::vector<T> &buf,* >> * const DataType &type,* >> * Request &req,* >> * const MessageTag &tag) const* >> *{* >> * START_LOG("send()", "Parallel");* >> >> *#ifndef NDEBUG* >> * // Only catch the return value when asserts are active.* >> * const int ierr =* >> *#endif* >> * ((this>send_mode() == SYNCHRONOUS) ?* >> * MPI_Issend : MPI_Isend) (buf.empty() ? NULL : &buf[0],* >> * cast_int<int>(buf.size()),* >> * type,* >> * dest_processor_id,* >> * tag.value(),* >> * this>get(),* >> * req.get());* >> >> * libmesh_assert (ierr == MPI_SUCCESS);* >> >> * STOP_LOG("send()", "Parallel");* >> *}* >> >> >> Suppose I want to send a vector<double> from processor 0 to processor 1. I >> tried the following: >> >> if (_communicator.rank() == 0) >> { >> std::vector<double> sender; >> sender.push_back(1.1); >> sender.push_back(2.2); >> Request & req; >> * _communicator.send(1, &sender, MPI::DOUBLE, req, 10.231);* >> } >> else if (_communicator.rank() == 1) >> { >> std::vector<double> receiver(2); >> Request & req; >> * _communicator.receive(0, &receiver, MPI::DOUBLE, req, 10.231);* >> } >> >> >> I encountered the compilation errors as follows. Does anyone know how to >> setup the correct syntax for this case? >> >> >> *'Request' was not declared in this scope'MPI' has not been declared* >> >> >> Thank you very much, >> Lei >> >> >>  >> Lei Zhao >> Research Assistant >> University of WisconsinMadison >> 206 Materials Science & Engineering >> 1509 University Avenue >> Madison, WI 53706, US >> Email: lzhao47@... >> Tel.: (608)3328385 >>  >>  >> Dive into the World of Parallel Programming The Go Parallel Website, >> sponsored >> by Intel and developed in partnership with Slashdot Media, is your hub >> for all >> things parallel software development, from weekly thought leadership >> blogs to >> news, videos, case studies, tutorials and more. Take a look and join the >> conversation now. http://goparallel.sourceforge.net/ >> _______________________________________________ >> Libmeshusers mailing list >> Libmeshusers@... >> https://lists.sourceforge.net/lists/listinfo/libmeshusers >> > 