<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class=""><br class=""><div><br class=""><blockquote type="cite" class=""><div class="">On Jan 19, 2019, at 12:58 AM, Valerio Barnabei <<a href="mailto:valerio.barnabei@gmail.com" class="">valerio.barnabei@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="">The only reason I may think about a slower code, is that
MatSetValues for MATIS has indirect memory access when mapping
global-to-local.</div></blockquote><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="">Are you preallocating the matrix properly? I
mean, are you calling MatISSetPreallocation with the same input used for
MatMPIAIJSetPreallocation?</div></blockquote><div class=""> </div><div class="">You're right. I'm preallocating a different number of nnz. I forgot how much preallocation influences performance. <br class=""></div><div class="">My doubt is: if i created my MATIS using MatCreateIS(PETSC_COMM_WORLD, 1, PETSC_DECIDE, PETSC_DECIDE,
ndofs
, ndofs, map, map, &K), how can i preallocate *_nnz for local submatrix? Do you suggest creating MATIS by passing the local size instead of the global size?<br class=""></div></div></div></div></blockquote><div><br class=""></div>The local size is not the size of the local FEM problems. It is the size of the local part of the global vector used in MatMult operations</div><div>If you want to do preallocation on the local FEM matrices, you do</div><div><br class=""></div><div>MatCreateIS(…&A)</div><div>MatISGetLocalMat(A,&B)</div><div>MatSeqAIJSetPreallocation(B,local_nonzeros…)</div><div><br class=""></div><div><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class=""><br class=""></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="">
Same loop, same entries. Only local indices in place of the global ones
</div></blockquote><div class=""><br class=""></div><div class="">And the performance is supposed to be the same, right?<br class=""></div></div></div></div></blockquote><div><br class=""></div>Same. Quoting Donald Knuth: <span style="font-variant-ligatures: normal; orphans: 2; widows: 2; background-color: rgb(255, 255, 255);" class="">Premature </span><span style="font-variant-ligatures: normal; orphans: 2; widows: 2; background-color: rgb(255, 255, 255);" class="">optimization</span><span style="font-variant-ligatures: normal; orphans: 2; widows: 2; background-color: rgb(255, 255, 255);" class=""> is the root of all </span><span style="font-variant-ligatures: normal; orphans: 2; widows: 2; background-color: rgb(255, 255, 255);" class="">evil</span><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class=""><br class=""></div><div class="">A couple more questions:</div><div class="">-I tried to use PCGAMG, and the convergence is really fast. Thank you for your suggestion. I'm planning to add an advection term to my equation, which preconditioner would you suggest?</div></div></div></div></blockquote><div><br class=""></div>Fast multigrid convergence with strong advecting terms is always tricky.</div><div><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class="">-I'm also trying to use
PCBDDC (as you suggested in an earlier mail), but i got errors, basically because I'm not sure about the options i'm supposed to use. If you think this PC could suite my case (considering the fact I will add an advection term soon), can I get any help to tweak it and make it work? Attached you'll find my error log of a run with -pc_type bddc. Is there any specific requirement (that i could have completely missed) in the domain decomposition strategy to use this particular PC?<br class=""></div></div></div></div></blockquote><div><br class=""></div>From the log, I see a bunch of errors coming from MatSetValues. The local-to-global map you have created is wrong.</div><div>I have an old branch where I create the l2g map brute force by doing a first pass of the values insertion loop (without actually inserting) and and then taking a unique list of global indices visited. It is not the most recommended way to create the map, but at least you have a comparison to understand where your code is wrong</div><div><a href="https://bitbucket.org/petsc/petsc/src/e1fa67b2d58bc4882f68a386d1578f90057aefbb/src/ksp/ksp/examples/tutorials/ex56.c?at=stefano_zampini/feature-pcgamg-useamat#lines-125" class="">https://bitbucket.org/petsc/petsc/src/e1fa67b2d58bc4882f68a386d1578f90057aefbb/src/ksp/ksp/examples/tutorials/ex56.c?at=stefano_zampini/feature-pcgamg-useamat#lines-125</a></div><div><br class=""></div><div>The only error with BDDC is not actually from BDDC, but from the local solvers, since your local matrices have a zero on the diagonal, and the default LU solver from PETSc does not tolerate this (no pivoting).<br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class=""><br class=""></div><div class="">Thank you for your help and your time.</div><div class=""><br class=""></div><div class="">Valerio<br class=""></div><div class=""><br class=""></div><div class=""><br class=""></div></div></div><br class=""><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Il giorno ven 18 gen 2019 alle ore 22:18 Stefano Zampini <<a href="mailto:stefano.zampini@gmail.com" class="">stefano.zampini@gmail.com</a>> ha scritto:<br class=""></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="overflow-wrap: break-word;" class=""><br class=""><div class=""><br class=""><blockquote type="cite" class=""><div class="">On Jan 18, 2019, at 6:47 PM, Valerio Barnabei <<a href="mailto:valerio.barnabei@gmail.com" target="_blank" class="">valerio.barnabei@gmail.com</a>> wrote:</div><br class="gmail-m_-7878666695508655580Apple-interchange-newline"><div class=""><div dir="ltr" class=""><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class="">Hello (again),<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class="">I'm trying to figure out how to change my program in order to use MATIS
type and PCBDD. <span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class="">I managed to build my MATIS matrix successfully, but my code is now slower.
I think the main issue is that i'm using MatSetValues inside the element loop,
and passing as values the e-th local element matrix.</span></div></div></div></blockquote><div class=""><br class=""></div>The only reason I may think about a slower code, is that MatSetValues for MATIS has indirect memory access when mapping global-to-local.</div><div class="">Are you preallocating the matrix properly? I mean, are you calling MatISSetPreallocation with the same input used for MatMPIAIJSetPreallocation?</div><div class=""><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class=""><span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class="">In a sort of pseudocode, where B is the shape functions matrix and D is a “diffusion”
matrix:<span class=""></span></span></div><p class="MsoNormal" style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:"Calibri",sans-serif"><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class=""><span class=""> </span></span></p><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">for e in localElements<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;text-indent:35.4pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">for qp in quadraturePoints{<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">
<span class=""> </span>for n in localNodes {<br class="">
<span class=""> </span>for m in localNodes {<span class=""> </span><br class="">
<span class=""> </span>ke += weights[qp]*detJ[qp]*B[n]^T*
D*B[m];<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt 106.2pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">}<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt 35.4pt;text-indent:35.4pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">}<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt 35.4pt;text-indent:35.4pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">For n in localNodes{<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt 70.8pt;text-indent:35.4pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">For m in localNodes{<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt 35.4pt;text-indent:35.4pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class=""><span class=""> </span>fe[n] -= ke[n* localNodes +m]*ge[e][m];//ge
is a simple array that lets me easily take into account Dirichlet conditions, for my simple case<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;text-indent:35.4pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">}<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;text-indent:35.4pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">for n in localNodes{<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;text-indent:35.4pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">//<span class="">
</span>idx [n]=LM [e][n];<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;text-indent:35.4pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">//LM is the DOF array (rows are Local elements,
columns are local nodes, values are global dofs, i.e. row/col indices for
global assembly of K, negative values are prescribed DOFs<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;text-indent:35.4pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">}<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">MatSetValues(K, localNodes,
(PetscInt*) &idx [0], localNodes, (PetscInt*) &idx [0], (PetscScalar*)
&ke[0], ADD_VALUES);<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">//same for
VecSetValues<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif;color:rgb(47,85,151)" lang="EN-GB" class="">}<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class="">Just to clarify: the array I used to create the map is a unique array for each
process, without the negative values of the prescribed DOFs; the array I’m
using in matsetvalues is each time the portion of the array used for the map relative
to the current element. In addition to that: the code is working, it gives the
same result of an older version where I used MATAIJ, but the portion of the
element assembly loop is slower, and so is the solving process (compared using
PCNONE on both, obtained same iterations to converge, same norm_2 of solution
vector).<span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class="">I guess using MatSetValues outside of the loop, using the same large array I
used to create the map for the matIS is supposed to be better, but I cannot see
a way to use it with this elemental assembly logic.</span></div></div></div></blockquote><div class=""><br class=""></div><div class="">MatSetValues (and MatSetValuesLocal) expects 2d data for the entries. Calling these once per element is the best you can do.</div><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class=""><span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class="">Stefano suggested to use MatSetValuesLocal: in that case what is supposed
to be idx? (I’m using global indices now).</span></div></div></div></blockquote><div class=""><br class=""></div>You should use the corresponding local indices. I mean, if your local -> global map is (0,1,2) -> (31,44,77), then, in place of 44, you should use 1 as idx for MatSetValuesLocal</div><div class=""><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class=""><span class=""></span></span></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class="">If I use MatSetValuesLocal, can I stick with this logic of adding values
inside the element loop?</span></div></div></div></blockquote><div class=""><br class=""></div>Same loop, same entries. Only local indices in place of the global ones</div><div class=""><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class="">Is there any other way you suggest to assign values
for this case?Any suggestion would be really appreciated.<br class=""></span></div><p class="MsoNormal" style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:"Calibri",sans-serif"><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class=""><span class=""></span></span></p><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class="">I’m not sure I explained my issue properly, I can give any further
information needed.</span></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class=""><br class=""></div><div style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:Calibri,sans-serif" class="">Valerio<br class=""></div><p class="MsoNormal" style="margin:0cm 0cm 0.0001pt;line-height:normal;font-size:11pt;font-family:"Calibri",sans-serif"><span style="font-size:12pt;font-family:"Times New Roman",serif" lang="EN-GB" class=""><span class=""></span></span></p><div style="margin:0cm 0cm 8pt;line-height:107%;font-size:11pt;font-family:Calibri,sans-serif" class=""><span lang="EN-GB" class=""><span class=""> </span></span><br class="gmail-m_-7878666695508655580webkit-block-placeholder"></div>
</div><br class=""><div class="gmail_quote"><div dir="ltr" class="gmail-m_-7878666695508655580gmail_attr">Il giorno ven 18 gen 2019 alle ore 15:55 Valerio Barnabei <<a href="mailto:valerio.barnabei@gmail.com" target="_blank" class="">valerio.barnabei@gmail.com</a>> ha scritto:<br class=""></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr" class=""><div class="">The mail was still incomplete. Sorry, i will write it again.</div><div class="">Apologies.<br class=""></div></div><br class=""><div class="gmail_quote"><div dir="ltr" class="gmail-m_-7878666695508655580gmail-m_-3579635411615542843gmail_attr">Il giorno ven 18 gen 2019 alle ore 15:53 Valerio Barnabei <<a href="mailto:valerio.barnabei@gmail.com" target="_blank" class="">valerio.barnabei@gmail.com</a>> ha scritto:<br class=""></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr" class=""><div dir="ltr" class=""><div dir="ltr" class=""><div dir="ltr" class=""><div dir="ltr" class=""><div dir="ltr" class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class="">Hello,</div><div class="">I'm trying to figure out how to change my program in order to use MATIS type and
PCBDD. <br class=""></div><div class="">I managed to build my MATIS matrix successfully, but my code is now slower. I think the main issue is that i'm using MatSetValues inside the element loop, and passing as values the e-th local element matrix. In a sort of pseudocode:</div><div class=""><br class=""></div><div class="">for e in localElements:{</div><div class="">for (int qp=0;qp<ngp_tot;qp++)<br class=""></div><div class=""> for(int A=0;A<Nnodesloc;A++){<br class=""> for(int B=0;B<Nnodesloc;B++){<br class=""> //Bat dot temp (NB: ba_n0=bat_0n)<br class=""> double prod=temp[0][A][qp]*BaT[B][0][qp]+temp[1][A][qp]*BaT[B][1][qp];<br class=""> //ke[A][B]+=weights[qp]*detJ[qp]*prod;<br class=""> ke[A*Nnodesloc+B]+=weights[qp]*detJ[qp]*prod;<br class=""></div><div class=""><br class=""></div><div class=""><br class=""></div></div><div dir="ltr" class=""><div class=""><div class="">?<br class=""></div><div class="">Actually i have lots of doubts:<br class=""></div><div class="">-when using
ISLocalToGlobalMappingCreate(), what's the blocksize?<br class=""></div>
</div></div><br class=""><div class="gmail_quote"><div dir="ltr" class="gmail-m_-7878666695508655580gmail-m_-3579635411615542843gmail-m_2365709834022617641gmail-m_3958850699414818303gmail-m_5556672655828192919gmail_attr">Il giorno lun 14 gen 2019 alle ore 15:08 Valerio Barnabei <<a href="mailto:valerio.barnabei@gmail.com" target="_blank" class="">valerio.barnabei@gmail.com</a>> ha scritto:<br class=""></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr" class=""><div class="">Thank you Stefano for your quick answer, I will try to change PC as you suggest, and I will read something about MATIS type. I will let you know if I still
have
doubts.<br class=""></div><div class=""><br class=""></div><div class="">Valerio<br class=""></div></div><br class=""><div class="gmail_quote"><div dir="ltr" class="">Il giorno dom 13 gen 2019 alle ore 16:08 Stefano Zampini <<a href="mailto:stefano.zampini@gmail.com" target="_blank" class="">stefano.zampini@gmail.com</a>> ha scritto:<br class=""></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr" class="">What problem are you trying to solve?<div class="">For standard elliptic problems (e.g. poisson or elasticity) you should use a preconditioner that scales, for example PCGAMG or BoomerAMG from HYPRE (PCHYPRE, if you configured PETSc with support for HYPRE via --download-hypre)<div class="">Alternatively, since you are using FEM with local Assembly, you can use PCBDDC, but this requires few extra calls</div><div class=""><br class=""></div><div class="">- call MatSetType(K,MATIS)</div><div class="">- provide an ISLocalToGlobalMapping of local to global degrees of freedom before you set preallocation via MatSetLocalToGlobalMapping (if you provide this call, then you can use MatSetValuesLocal for both AIJ amd MATIS types)</div></div><div class=""><br class=""></div><div class=""><br class=""></div></div><br class=""><div class="gmail_quote"><div dir="ltr" class="">Il giorno dom 13 gen 2019 alle ore 17:59 Valerio Barnabei via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank" class="">petsc-users@mcs.anl.gov</a>> ha scritto:<br class=""></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr" class="">
<div class="">Hello everybody,<br class=""></div><div class=""><br class=""></div><div class="">I have some doubts about the parallel solution of a simple FEM ksp problem. I searched in the mailing list but the information I found are not helping me.<br class=""></div><div class="">I
need to control manually the parallel assemble of the main matrix, but I
would like to use PetSC for the solution of the final linear system
K*U=F. <br class=""></div><div class="">The algorithm I am using, loads an
already decomposed FEM domain, and each process contributes to the
global matrix assembly using its portion of domain decomposition (made
in pre-processing). In the element loop, at each element I load the
element submatrix in the global matrix via MatSetValue function. Then,
out of the element loop, I write<br class=""></div><div class="">/*</div><div class="">...matrix preallocation<br class=""></div><div class="">...calculate element local stiffness matrix (local in a FEM sense and local in a domain decomposition sense)<br class=""></div><div class="">MatSetValue (ADD)</div><div class="">*/<br class=""></div><div class="">and then <br class=""></div><div class=""><span style="color:rgb(0,0,255)" class="">ierr=MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br class="">ierr=MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br class="">ierr=VecAssemblyBegin(F);CHKERRQ(ierr);<br class="">ierr=VecAssemblyEnd(F);CHKERRQ(ierr);</span><br class=""></div><div class=""><br class=""></div><div class="">This part scales perfectly and is sufficiently fast.</div><div class="">
My K matrix is a sort of banded diagonal matrix: my preallocation is not perfect, but it works at this stage of my work.
</div><div class="">Finally I solve the system with the ksp logic:</div><div class=""><br class=""></div><div class=""><span style="color:rgb(0,0,255)" class="">ierr=KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);<br class="">ierr=KSPSetOperators(ksp,K,K);CHKERRQ(ierr);<br class="">ierr=KSPGetPC(ksp,&pc);CHKERRQ(ierr);<br class="">ierr=PCSetType(pc,PCBJACOBI);CHKERRQ(ierr); //or PCNONE, or PCJACOBI, the behaviour is almost identical, the iteration number lowers but is still really high<br class="">ierr=KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);<br class="">ierr=KSPSetFromOptions(ksp);CHKERRQ(ierr);<br class="">ierr=KSPSolve(ksp,F,sol);<br class=""></span></div><div class=""><br class=""></div><div class="">When
I run the code I am noticing some strange behaviours. Indeed, the
numerical solution appear right, and the code scales if I make a simple
strong-scaling test. However, the ksp_log and log_view (see below, logs for a 160k quad linear element uniform 2dmesh) shows me a huge number of iterations, and the
time spent in "solve" function is very high. Furthermore, the more I
increase the problem size keeping a constant load per process (weak
scaling) the more iterations seems to be needed for convergence. In general, when increasing the problem size (i.e. 16M elements) the convergence is not achieved before the maximum iteration number is achieved. <br class=""></div><div class=""><br class=""></div><div class="">Question 1: Am I doing it right? Or am I making some trouble with the MPI communication management? (i'm aware that my domain decomposition does not match the petsc matrix decomposition, I'm expecting it to cripple my algorithm a bit)<br class=""></div><div class="">Question 2: Why the solver takes so much time and so many iterations? the problem is really simple, and the solution appears correct when plotted. Am I ill conditioning the system?</div><div class="">Thanks in advance for your help. I can add any further information if needed.</div><div class=""><br class=""></div><div class="">Valerio<br class=""></div>
</div>
</blockquote></div><br clear="all" class=""><div class=""><br class=""></div>-- <br class=""><div dir="ltr" class="gmail-m_-7878666695508655580gmail-m_-3579635411615542843gmail-m_2365709834022617641gmail-m_3958850699414818303gmail-m_5556672655828192919gmail-m_-8752206085825601650gmail-m_-5858931685459960259gmail_signature">Stefano</div>
</blockquote></div>
</blockquote></div>
</div></div></div></div></div></div></div>
</blockquote></div>
</blockquote></div>
</div></blockquote></div><br class=""></div></blockquote></div>
<span id="cid:f_jr2kzd0b0"><log.out></span></div></blockquote></div><br class=""></body></html>