node15.html

htmdoc for html coding

linear systems, amongst which are:
   <UL><LI> Direct methods, such as Gauss, Crout and Cholesky methods, where the matrix is stored in main or
     secondary memory ([<A HREF="node65.html#mod_5"><A NAME=tex2html39 HREF="../Guide5/welcome.html">MODULEF User Guide - 5</A></A>],<BR>[4] [<A HREF="node65.html#mod_24">24</A>].
    <LI> Iterative methods, such as relaxation and conjugate gradient with or without preconditioning  
     ([<A HREF="node65.html#mod_5"><A NAME=tex2html39 HREF="../Guide5/welcome.html">MODULEF User Guide - 5</A></A>], [<A HREF="node65.html#mod_102">102</A>]).                                                       
   </UL>
<P>
In all these cases the matrix corresponding to the finite element discretization, is sparse. In order to minimize
the space occupied in memory, an appropriate method of matrix storage is sought:
  <UL><LI> For direct methods (excluding hyper-matrices  [<A HREF="node65.html#mod_62">Van-Ingelandt-1985</A>]) skyline storage is used (see description 
   of the D.S. <b> MUA</b> in <BR>[4][<A HREF="node65.html#guide_2"><A NAME=tex2html35 HREF="../Guide2/welcome.html">MODULEF User Guide - 2</A></A>]).
   <LI> For iterative methods compact storage is used where only non-zero coefficients are stored (see description
   of the D.S. <b> AMAT</b> in <BR>[4][<A HREF="node65.html#guide_2"><A NAME=tex2html35 HREF="../Guide2/welcome.html">MODULEF User Guide - 2</A></A>]).
  </UL>
<P>
The non-linear problem can often be solved by regarding it as that of a series of linear problems, each
treated as described above.
<P>
The choice of solution method for linear systems is not obvious, and depends strongly on the nature of the problem
to be solved. In the case of simple systems, all methods produce good results. When there are many systems 
to solve, for example an iterative process, there are some points in favour of both types of solution methods:
  <UL><LI> Advantage of the direct method: If the matrix is constant, it is factorized only once, afterward the
   solution requires only one forward- and backward-substitution.                       
   <LI> Advantage of the iterative method: A good approximation to the solution is known at the previous
    iteration or at the previous time-step, and  the number of iterations necessary is therefore reduced.
  </UL>
<P>
In order to analyze the different sub-steps, consider the simple example of solving a linear system by a direct
Cholesky method in main memory (m.m.).
<P>
  <OL><LI> <b> Construction of pointers</b> <BR> <A NAME=steppoin>&#160;</A>
   This step constructs the pointers required for the skyline storage by looking at all the elements in the mesh, and 
   calculates the amount of memory to be occupied by the matrix. Computation will therefore only be started if 
   the memory space available is sufficient.
<P>
   This step is performed by module <b> PREPAC</b> [<A HREF="node65.html#mod_5"><A NAME=tex2html39 HREF="../Guide5/welcome.html">MODULEF User Guide - 5</A></A>]. The result is found in the D.S. <b> MUA</b>, 
   which does not contain the matrix but only the pointers.
<P>
  <LI> <b> Assembly</b> <BR> 
   Using the D.S. <b> MUA</b> created in the previous step, and the D.S. <b> TAE</b> containing the element matrices 
   and  arrays,
   the modules <b> ASSMUA</b> and <b> ASEMBV</b> respectively assembles the element matrices and the 
   right-hand-side vectors, i.e., computes the following equations (see equation (<A HREF="node15.html#eq7">1.9</A>)):
     <DIV ALIGN=center><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img108.gif"></DIV>  
     and
     <DIV ALIGN=center><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img109.gif"></DIV>
<P>
   The results of this step is found in the D.S. <b> MUA</b>, containing the system matrix, and the D.S. <b> B</b>,
   containing the right-hand-side vector.
<P>
  <LI> <b> Implementation of imposed boundary conditions</b> <BR> 
  Boundary conditions fall into the following two categories:   
  <OL><LI> those which are implemented in the variational formulation, for example the Neumann or Fourier conditions 
     in the classical thermal problem (see [<A HREF="node65.html#mod_100">100</A>]: Introduction), and
<P>
    <LI> those which are not taken into account in the variational formulation, for example the Dirichlet
     conditions.
  </OL>
<P>
  The purpose of this  step, which consists of two sub-steps, is the implementation of the second type of boundary 
  conditions.
<P>
    <OL><LI> <b> Construction of the D.S. BDCL</b> <BR> 
    The D.S. <b> BDCL</b> contains a description of the imposed boundary conditions (numbers  and values of 
    restricted degrees of freedom).
    This data structure is created by the modules <b> COBDC1</b> or <b> COBDCL</b> [<A HREF="node65.html#mod_18">18</A>].   The module
    <b> COBDCL</b> constructs the D.S. <b> BDCL</b>, degree of freedom by degree of freedom, 
    from the data cards. Module <b> COBDC1</b> constructs the same data structure by utilizing global notions, for
    example reference numbers and the names of the variational unknowns, and is simpler to use than <b> COBDCL</b>.
    Both modules treat boundary conditions defined as linear relations, which is particularly useful when periodic 
    conditions are encountered.
<P>
    <b> Remark:</b> 
     <OL><LI> If the problem under consideration contains boundary conditions defined as linear relations, this step 
     must be performed before step <A HREF="node15.html#steppoin">1</A>, as the matrix profile is modified.     
     <LI> Module <b> COBDC1</b> utilizes the D.S. <b> COOR</b>, which contains the nodal coordinates. This data
     structure can be generated by the module <b> CORNOE</b> [<A HREF="node65.html#mod_63">63</A>].                                               
    </OL>
<P>
    <LI> <b> Effective implementation of the boundary conditions</b> <BR> 
     This step is concerned with the derivation from <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img97.gif"> to <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img110.gif">, considering equations (<A HREF="node15.html#eq9">1.11</A>) and
     (<A HREF="node15.html#eq10">1.12</A>).
<P>
     There are several techniques for implementing boundary conditions, as described  
     in [<A HREF="node65.html#mod_5"><A NAME=tex2html39 HREF="../Guide5/welcome.html">MODULEF User Guide - 5</A></A>]: module <b> CLIMPC</b>, for skyline matrices. Conceptually, we can distinguish two classes:
     <OL><P>
<P>
       <LI> <b> The penalized method: </b>
<P>
       <P>
<P>
       The point of departure is the following &quot;computer&quot; equality:
       <DIV ALIGN=center><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img111.gif"></DIV>
       if <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img112.gif"> (for example: <IMG BORDER=0 ALIGN=BOTTOM ALT="" SRC="img113.gif"> and <IMG BORDER=0 ALIGN=BOTTOM ALT="" SRC="img114.gif">).  <BR>
<P>
       In practice we start with <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img97.gif"> and <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img115.gif">  of equation (<A HREF="node15.html#eq9">1.11</A>), and replace the 
       diagonal coefficients of <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img116.gif"> by <b>VTG</b> and the coefficients of <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img117.gif"> by <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img118.gif">, where
       <P><A NAME=eqU0>&#160;</A><IMG BORDER=0 ALIGN=BOTTOM ALT="" SRC="img119.gif"><P>  
       where
       <DL COMPACT><DT><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img120.gif">
<DD> is a vector of values of imposed degrees of freedom equal to zero
                     in this example (see (<A HREF="node15.html#eq10">1.12</A>)).
<P>
 </DL>
<P>
       <P>
<P>
       A variation on the above method consists of replacing the diagonal coefficients of <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img116.gif"> by:
       <DIV ALIGN=center><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img121.gif"></DIV>
       and the right-hand-side <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img117.gif"> by:                 
       <DIV ALIGN=center><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img122.gif"></DIV>
       where <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img123.gif">.  <BR>
<P>
       <P>
<P>
       Advantages and disadvantages:
<P>
       <UL><LI> In the first variation, the initial values of the  diagonal coefficients are lost.
         <LI> In the second variation, it is necessary to know the values of the diagonal coefficients in order
          to take the boundary conditions into account in the right-hand-side vector. It is therefore necessary to 
          do both incorporations of boundary conditions simultaneously, which could be inconvenient in an iterative 
          method where the matrix is constant.
        </UL> 
       <P>
       To summarize:
       <DL COMPACT><DT>Variation 1
<DD> incorporates the boundary conditions into the matrix and right-hand-side vector 
            separately <BR>[4] ([<A HREF="node65.html#mod_5"><A NAME=tex2html39 HREF="../Guide5/welcome.html">MODULEF User Guide - 5</A></A>]: <b> PROFIL</b> and <b> CLIMPC</b>), with
<P>
           <DIV ALIGN=center><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img124.gif"></DIV>
            and
             <DIV ALIGN=center><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img125.gif"></DIV>
<P>
         <DT>Variation 2
<DD> incorporates the boundary conditions simultaneously in <b> MUA</b> and <b> B</b>, with
<P>
           <DIV ALIGN=center><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img126.gif"></DIV> 
 </DL>
<P>
       <P>
<P>
      <LI> <b> The direct method:</b>
<P>
       <P>
<P>
       The idea is to derive equation (<A HREF="node15.html#eq10">1.12</A>) from (<A HREF="node15.html#eq9">1.11</A>) effectively. Note that in equation
       (<A HREF="node15.html#eq10">1.12</A>) an 
       important property, namely the symmetry of the matrix, is lost. The importance lies in the property that 
       for a symmetric matrix only half the matrix needs to be stored. We therefore want to derive (<A HREF="node15.html#eq12">1.15</A>),
       which is equivalent to equation (<A HREF="node15.html#eq10">1.12</A>), from equation (<A HREF="node15.html#eq9">1.11</A>), where (<A HREF="node15.html#eq12">1.15</A>) is symmetric.
<P>
       In the general case (imposed d.o.f. has a value of <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img120.gif">, zero or nonzero (<A HREF="node15.html#eqU0">1.14</A>)), if we consider the
       decomposition of <IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img30.gif">:
       <DIV ALIGN=center><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img127.gif"></DIV>
       where 
       <DL COMPACT><DT><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img128.gif">
<DD> is the vector of free degrees of freedom, and
         <DT><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img120.gif">
<DD> is the vector of the blocked degrees of freedom,
<P>
 </DL>
<P>
        then equation (<A HREF="node15.html#eq10">1.12</A>) can be written as follows:
        <DIV ALIGN=center><IMG BORDER=0 ALIGN=MIDDLE ALT="" SRC="img129.gif"></DIV>
<P>
        which is equivalent to:
<P>
        <P><A NAME=eq12>&#160;</A><IMG BORDER=0 ALIGN=BOTTOM ALT="" SRC="img130.gif"><P>
<P>
        <P>
<P>
        In practice, it is equation (<A HREF="node15.html#eq12">1.15</A>) that is calculated (using the notation in [<A HREF="node65.html#mod_5"><A NAME=tex2html39 HREF="../Guide5/welcome.html">MODULEF User Guide - 5</A></A>]:
         <b> PROFIL</b> and <b> CLIMPC</b>, with <b>NIVO = 3</b> and <b>VTG = 1</b>). <BR>
<P>
     </OL>
       <P>
      <b> Remark:</b> The different techniques above are independent of the matrix storage technique, and can
        therefore be used in other methods of storage.
<P>
    </OL>
<P>
    <P>
<P>
  <LI> <b> Cholesky factorization</b> <BR> 
    This step computes the triangular matrix <b>L</b> which satisfies the relationship:
    <DIV ALIGN=center><IMG BORDER=0 ALIGN=BOTTOM ALT="" SRC="img131.gif"></DIV>
    and is performed by module <b> CHOLPC</b> [<A HREF="node65.html#mod_5"><A NAME=tex2html39 HREF="../Guide5/welcome.html">MODULEF User Guide - 5</A></A>].
<P>
    The result is stored in the D.S. <b> MUA</b> which contains the factorized matrix.
<P>
  <LI> <b> Forward- and backward-substitution</b> <BR> 
    The forward- and back-substitution solves two linear systems of triangular matrices.  It
    is performed by module <b> DRCHPC</b> <BR>[4][<A HREF="node65.html#mod_5"><A NAME=tex2html39 HREF="../Guide5/welcome.html">MODULEF User Guide - 5</A></A>].
<P>
    The result is stored in the D.S. <b> B</b>.
<P>
  </OL>
<P>
    <b> Comments:</b> <BR> 
    The solution of the linear system is decomposed into several steps in order to simplify the computation as
    much as possible. For example:
    <UL><LI> If several systems sharing the same matrix are solved, the factorization is done only once, and  
          similarly,
    <LI> when assembling a stiffness and mass matrix, the preparation (calculation of pointers) is done 
          only once, etc.
    </UL>
<P>
<P>
<P>
<DT>Step 5:
<DD><b> Interpretation of results</b> <BR>
<P>
   The graphic visualization of results is particularly important. Indeed, it is difficult to interpret the 
   results by merely
   printing the output data structure <b> B</b>. In order to aid you with the interpretation of results, the following
   tools are available:
   <UL><LI> The computation of error norms and the presentation of the exact and computed solutions, if the analytical 
      solution is known (see [<A HREF="node65.html#mod_22">22</A>]: module <b> NORME</b>).
      <LI> Plotting of isovalues (in the case of the present example), stresses,  velocities, deformations, etc.
       [<A HREF="node65.html#mod_96"><A NAME=tex2html41 HREF="../Guide6/welcome.html">MODULEF User Guide - 6</A></A>].
   </UL>
<P>
   In the post-treatment, we can also utilize the modification modules for the data structures
   <b> NOPO</b> [<A HREF="node65.html#mod_104"><A NAME=tex2html37 HREF="../Guide3/welcome.html">MODULEF User Guide - 3</A></A>] and <b> B</b>. 
   For example, if the problem is symmetric, the computation need only be done on part 
   of the domain. However, it is desirable to present the results on the complete domain, by &quot;gluing&quot; together
   the sub-domains.
<P>
 </DL>
<P>
<P>
<P>
<H2><A NAME=SECTION04114000000000000000>1.1.4 Conclusions</A></H2>
<P>
<P>
<P>
The purpose of this example was to illustrate how find the suitable modules in the libraries, capable of solving the
problem under consideration. However, it does not show how to use a given module. The examples treated in the 
remainder of this report will illustrate this process in detail.
<P>
<P>
<P>
<HR SIZE=3 WIDTH="75&#37;"><IMG SRC="../icons/smallmod.gif" WIDTH=211 HEIGHT=50 ALIGN=BOTTOM
	ALT="Modulef"><A NAME=tex2html379 HREF="node14.html"><IMG BORDER=0 ALIGN=BOTTOM SRC="../icons/previous_motif.gif"
	ALT="previous"></A><A NAME=tex2html383 HREF="node14.html"><IMG BORDER=0 ALIGN=BOTTOM SRC="../icons/up_motif.gif"
	ALT="up"></A><A NAME=tex2html385 HREF="node16.html"><IMG BORDER=0 ALIGN=BOTTOM SRC="../icons/next_motif.gif"
	ALT="next"></A><A NAME=tex2html387 HREF="node2.html"><IMG BORDER=0 ALIGN=BOTTOM SRC="../icons/contents_motif.gif"
	ALT="contents"></A><A HREF="../Guide1-18/node15.html"><IMG BORDER=0 SRC="../icons/zoom18.gif" ALIGN=BOTTOM
	ALT="[BIG]"></A><A HREF="../Guide1-14/node15.html"><IMG BORDER=0 SRC="../icons/zoom14.gif" ALIGN=BOTTOM
	ALT="[Normal]"></A><A HREF="../Guide1-10/node15.html"><IMG BORDER=0 SRC="../icons/zoom10.gif" ALIGN=BOTTOM
	ALT="[small]"></A><BR>
<B> Next: </B> <A NAME=tex2html386 HREF="node16.html">2 Conversational tests</A>
<B>Up: </B> <A NAME=tex2html384 HREF="node14.html">1 Introduction</A>
<B> Prev: </B> <A NAME=tex2html380 HREF="node14.html">1 Introduction</A>
<B><A HREF="node2.html"
	>Contents</A></B>
<BR> <HR>
<P><ADDRESS>
</ADDRESS>
</BODY>
</HTML>