<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Fri, Oct 23, 2015 at 3:01 AM, Afanasiev Michael <span dir="ltr"><<a href="mailto:michael.afanasiev@erdw.ethz.ch" target="_blank">michael.afanasiev@erdw.ethz.ch</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div style="word-wrap:break-word">
Hi Matthew,
<div><br>
</div>
<div>So I’m discretizing via a tensor product of Lagrange polynomials co-located at the Gauss-Lobatto-Legendre (GLL) points. The polynomial order might range from 4 to 10 or something like that. The current problem is solved on 2D hexes. </div>
<div>
<div><br>
</div>
<div>I had found DMPlexCreateSection, and followed plex/ex1to get things set up. You can see my attempt below. Basically I’ve got 4 fields (displacement, velocity, acceleration, and strain) over each element. Here I’ve tried to setup a tensor product
of 4th order (5 GLL points) Lagrange polynomials (line 11). This seemed to <i>
somewhat</i> achieve what I wanted — I created a global vector and wrote it to a vtk file with VecView, and the numbering seemed to make sense. You can also see my attempt at defining a boundary condition (it looked like DMPlexCreateFromExodus labeled side
sets as “Face Sets”, seems to have worked). </div>
<div><br>
</div>
<div>Does this seem to be headed in the right direction?</div></div></div></blockquote><div><br></div><div>Yes, however I have some questions. Starting out, I think the GLL points include the endpoints, so</div><div>that means for polynomial degree k that numDof[] should really have 4*1 dofs on each vertex,</div><div>4*(k-1) dofs on each edge, and 4*(k-1)^2 on each quad. It looks like below you have numDof[] for</div><div>a 1D mesh with DG element.</div><div><br></div><div>The "Face Sets" is the right label to use for boundary conditions. This will eliminate those variables</div><div>from the global system, but they will be present in the local spaces.</div><div><br></div><div>With elements like these, it is common (I think) to eliminate the cell unknown explicitly, since the</div><div>system is dense and not connected to other cells. For this, you would retain the vertex and edge</div><div>unknowns, but not the cell unknowns. I have not tried to do this myself, so I do not know if there</div><div>are any pitfalls.</div><div><br></div><div>You can see an example of a similar implementation specifically for the kind of spectral elements</div><div>you are considering here: <a href="https://github.com/jedbrown/dohp">https://github.com/jedbrown/dohp</a>. It would probably be useful to understand</div><div>what is done there as you implement.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div style="word-wrap:break-word"><div>
<div>Cheers,</div>
<div>Mike.</div>
<div><br>
</div>
<div>
<div>DM</div>
<div>mesh::createSection(const DM &dm)</div>
<div>{</div>
<div><br>
</div>
<div>01 // displacement, velocity, acceleration, strain</div>
<div>02 IS bcPointIs[1];</div>
<div>03 PetscInt numBc = 1;</div>
<div>04 PetscInt bcField[1];</div>
<div>05 PetscInt numFields = 4;</div>
<div>06 PetscInt dim; DMGetDimension(dm, &dim);</div>
<div>07 PetscInt numComp[numFields];</div>
<div>08 for (auto i=0; i<numFields; i++) {numComp[i] = dim;}</div>
<div>09 PetscInt numDof[numFields*(dim+1)];</div>
<div>10 for (auto i=0; i<numFields*(dim+1); i++) {numDof[i] = 0;}</div>
<div>11 for (auto i=0; i<numFields; i++) {numDof[i*(dim+1)+dim] = 5;}</div>
<div>12 bcField[0] = 0;</div>
<div>13 PetscSection section;</div>
<div>14 DMPlexGetStratumIS(dm, "Face Sets", 1, &bcPointIs[0]);</div>
<div>15 DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBc, bcField,</div>
<div>16 NULL, bcPointIs, NULL, §ion);</div>
<div>17 ISDestroy(&bcPointIs[0]);</div>
<div>18 PetscSectionSetFieldName(section, 0, "u");</div>
<div>19 PetscSectionSetFieldName(section, 1, "v");</div>
<div>20 PetscSectionSetFieldName(section, 2, "a");</div>
<div>21 PetscSectionSetFieldName(section, 3, "e");</div>
<div>22 DMSetDefaultSection(dm, section);</div>
<div>23 return dm;</div>
<div>}</div>
</div>
<div><span class=""><br>
<div>--<br>
Michael Afanasiev<br>
Ph.D. Candidate<br>
Computational Seismology<br>
Institut für Geophysik<br>
ETH Zürich<br>
<br>
Sonneggstrasse 5, NO H 39.2<br>
CH 8092 Zürich<br>
<a href="mailto:michael.afanasiev@erdw.ethz.ch" target="_blank">michael.afanasiev@erdw.ethz.ch</a><br>
</div>
<br>
</span><div><div class="h5"><div>
<blockquote type="cite">
<div>On Oct 22, 2015, at 2:32 AM, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:</div>
<br>
<div>
<div dir="ltr" style="font-family:Optima-Regular;font-size:12px;font-style:normal;font-variant:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px">
<div class="gmail_extra">
<div class="gmail_quote">On Wed, Oct 21, 2015 at 3:07 PM, Dave May<span> </span><span dir="ltr"><<a href="mailto:dave.may@erdw.ethz.ch" target="_blank">dave.may@erdw.ethz.ch</a>></span><span> </span>wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div dir="ltr">Hey Mike,<br>
<br>
<br>
<div class="gmail_extra"><br>
<div class="gmail_quote">On 21 October 2015 at 18:01, Afanasiev Michael<span> </span><span dir="ltr"><<a href="mailto:michael.afanasiev@erdw.ethz.ch" target="_blank">michael.afanasiev@erdw.ethz.ch</a>></span><span> </span>wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex">
<div style="word-wrap:break-word">Hey Dave,
<div><br>
</div>
<div>So I’ve got a couple of days where there’s nothing pressing to work on… was thinking of ripping out the parallel routines (ugly) in my wave propagation code and replacing them with Petsc DM routines. I can read in my exodusii mesh with DMPLEX,
and partition it with chaco, which gives me a nicely partitioned DM. This takes me like 5 lines of code. That’s amazing.</div>
<div><br>
</div>
<div>But here I’m stuck, and am having a whale of a time with the documentation. All I<span> </span><i>think</i><span> </span>I need is a way to modify the exodus-created DM, and
add to it the degrees of freedom that are introduced by my quadrature rule. This would be really neat. I can just treat each sub-domain as its own mesh, with its own global numbering. Then whenever necessary I can scatter stuff the the<span> </span><i>real</i> global
degrees of freedom with something like VecLocToGlob. Most of the things I like about the code could stay the same (element-wise, matrix-free nature), just these parallel broadcasts would be infinitely nicer.</div>
<div><br>
</div>
</div>
</blockquote>
<div><br>
</div>
<div>First off - I don't use DMPLEX.</div>
</div>
</div>
</div>
</blockquote>
<div><br>
</div>
<div>Dave is refreshingly candid about his shortcomings ;)</div>
<div> </div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div dir="ltr">
<div class="gmail_extra">
<div class="gmail_quote">
<div> </div>
</div>
</div>
</div>
</blockquote>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div dir="ltr">
<div class="gmail_extra">
<div class="gmail_quote">
<div></div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex">
<div style="word-wrap:break-word">
<div></div>
<div>But I just can’t figure out how to set this up. The main problem really boils down to: what’s the best way to add my quadrature points to an already-created DM, which was constructed with an exodus file? I guess I could do this after the file
is read, but before the partitioning. In this case though, what’s stopping the partitioner from cutting an element in half? It seems like it would be a lot cleaner to do this post-partitioning.</div>
<div><br>
</div>
</div>
</blockquote>
<div><br>
</div>
<div>Presumably what is read from exodus is just the vertices of the hexes, and what you want to do is define the function space (given by your GLL locations) on top of element geometry read in. Is that what you are asking about?</div>
</div>
</div>
</div>
</blockquote>
<div><br>
</div>
<div>So Dave is right. We read in topology and geometry from ExodusII. Then you define a function space on top. How</div>
<div>exactly are you discretizing? In order to create vectors, do local to global, etc. Petsc really only need to know the</div>
<div>amount of data associated with each mesh piece. You can define this with PetscSection. If you give me an idea</div>
<div>what you want I can help you write the code easily I think.</div>
<div><br>
</div>
<div> Thanks,</div>
<div><br>
</div>
<div> Matt</div>
<div> </div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div dir="ltr">
<div class="gmail_extra">
<div class="gmail_quote">
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex">
<div style="word-wrap:break-word">
<div></div>
<div>Any hints here?</div>
</div>
</blockquote>
<div><br>
</div>
<div>Actually I have no experience with this object.<br>
I would just send an email to<span> </span><br>
</div>
<div> <span> </span><a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a><br>
</div>
<div>asking for help.<br>
</div>
<div> <br>
</div>
<div>The developer of DMPLEX (Matt Knepley) will definitely answer within in 1 day.<br>
</div>
<div><br>
</div>
<div>Cheers,<br>
</div>
<div> Dave<br>
</div>
<div><br>
</div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex">
<div style="word-wrap:break-word">
<div><br>
</div>
<div>Best,</div>
<div>Mike.<span><font color="#888888"><br>
<div>--<br>
Michael Afanasiev<br>
Ph.D. Candidate<br>
Computational Seismology<br>
Institut für Geophysik<br>
ETH Zürich<br>
<br>
Sonneggstrasse 5, NO H 39.2<br>
CH 8092 Zürich<br>
<a href="mailto:michael.afanasiev@erdw.ethz.ch" target="_blank">michael.afanasiev@erdw.ethz.ch</a><br>
</div>
<br>
</font></span></div>
</div>
</blockquote>
</div>
<br>
</div>
</div>
</blockquote>
</div>
<br>
<br clear="all">
<div><br>
</div>
--<span> </span><br>
<div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener</div>
</div>
</div>
</div>
</blockquote>
</div>
<br>
</div></div></div>
</div>
</div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>