text:
# Compact storage of K

As we have seen in POISS and ELAST, the standard f.e. program needs to assemble the global stiffness matrix K from the element matrices written to scratch file by subroutine stiff. But storing the whole n \times n matrix K would be impractical. Luckily, we can make use of the fact that K will be banded, as you have seen in the graphical demonstration of the assembly process in PREFEL, with zeros in entries outside the band (which lies either side of the diagonal of the matrix). In many applications K will also be symmetric, so that we only need to store the lower (or upper) half of the band, together with the diagonal itself.

## Symmetric band storage

The simplest example of band storage is in POISS.FOR, since in this application there is only one degree of freedom per node; hence, the equation corresponding to the i th node is in the i th row of the stiffness matrix equation. The stiffness matrix will thus have NPOIN rows.

The size of the band is defined by the bandwidth NBW. This defined as the smallest number for which: K_{i,j}=0 for all i,j satisfying |i-j| > {NBW} . If K is symmetric, we will store the lower half of the band, i.e. all entries K_{i,j} for which i gt;= j . The halfbandwidth NHBW is defined by NBW = 2*NHBW + 1. Since we only get a nonzero entry in K_{i,j} if nodes with global numbering i and j are contained in the same element, the halfbandwidth NHBW can be calculated as the maximum over all elements l , of the difference in node numbers for the nodes in l . In POISS.FOR, this calculation is done at the start of subroutine assmb. We will therefore use NPOIN rows of GSTIF (one row for each node in the mesh) and NHBW+1 columns. We denote the latter value by NBAND, and dimension GSTIF at the start of POISS.FOR as
DIMENSION GSTIF(MPOIN,MBAND)
where the parameters MPOIN and MBAND are set in the PARAMETER statement.

Then, on row i we store only k_{i,i- {NHBW}},...,k_{i,i} , for i=1,... ,NTOTV. We store these entries in an array GSTIF with NTOTV rows (where NTOTV is the total number of degrees of freedom) and (NHBW+1) columns. The mapping from K to GSTIF is: n K_{i,j} --> {GSTIF(I,J-I+NHBW+1)} for j=i- {NHBW},...,i .

Then the diagonal of the matrix (entries where i=j ) will be stored in the final column of the array. You can see the entries of the element stiffness matrices being placed in GSTIF, in subroutine assmb of POISS.FOR.

The mapping is slightly more complicated in ELAST.FOR, where there are two degrees of freedom per node. Here, the u and v d.o.f.s of node n (i.e. displacements in the x and y directions at node n ) are global degrees of freedom numbers i_n+1 and i_n+2 , where i_n=2(n-1) . You can see this mapping being used at the start of subroutine assmb. After making this adjustment, the local d.o.f.s are mapped to GSTIF by the formula in the previous section. At the start of assmb, the halfbandwidth is found by finding the position of the first non-zero entry in row I: this is stored in the array NMINVC. Thus, if NMINVC(I)=MI, it means that k_{i,mi} \neq 0 but k_{i,j}=0 {for} j=1,...,mi-1 . Then the halfbandwidth NHBW is given by \max_{i=1, {NTOTV}}(i-mi) .

It will be seen that the size of the bandwidth, and hence the amount of storage required, will depend on the order in which the nodes of the mesh have been numbered. For this reason it is very important to use the element renumbering facility in PREFEL when producing the .dat file for input to the main engine', to obtain a sensible numbering of the elements --- and then to re-number the nodes accordingly. Note that the bandwidth does not depend on the element numbering directly.

## Non-symmetric band storage

When the matrix K is nonsymmetric, the whole band must be stored. This is the case in THERM.FOR. Here, the array GSTIF has 3N rows (as there are 3 degrees of freedom per node) and NBW columns (NBW is the bandwidth). In this program we introduce the use of arrays NVAR and NDOF, both of size NPOIN. NDOF(N) stores the number of degrees of freedom at node N. NVAR(N) tells us which rows of the global matrix K these degrees of freedom correspond to: if NVAR(N)=IN, then the i'th d.o.f. of node N corresponds to row IN+i. NVAR is formed from NDOF by NVAR(N) = \sum_{i=1}^{N-1} NDOF(i). These arrays are used in other main engines', especially where there are variable degrees of freedom per node.

The mapping from K to GSTIF is again n K_{i,j} --> {GSTIF(I,J-I+NHBW+1)}. Now, however, there are NBW columns, and the diagonal entries of K lie down the (NHBW+1)'th column.

## Skyline storage

While band storage is considerably better than storing the whole matrix, it still leaves a large number of zero entries within the band. A more efficient storage system is skyline storage, where we decide separately for each row where we need to start the storage.

This storage technique is used in the soil consolidation program CONSL.FOR. As with band storage, the amount of storage used depends on the numbering of the nodes in the mesh. We again use the array NMINVC(I) to tell us the column where the first non-zero entry on row I of K lies. Then (as K will by symmetric in CONSL) we store only the entries from this point up to the diagonal entry, in GSTIF. GSTIF now becomes a one-dimensional array, and we use a pointer array LPOINT(I) to tell us where the I'th row of K begins. That is, the diagonal entry K_{i,i} will be stored as GSTIF(LPI), where LPI=LPOINT(I), for I=1 to NTOTV. The dimension of GSTIF is then LPOINT(NTOTV), which we denote by NSTIF. The parameter MSTIF is used to dimension GSTIF at the start of the program. The arrays NMINVC and LPOINT are set up by a subroutine mask in CONSL.FOR. LPOINT(I) can be found from LPOINT(I-1) by
LPOINT(I) = LPOINT(I-1) + I - NMINVC(I) + 1
since there are I-NMINVC(I)+1 entries to store on row I. We start with LPOINT(1)=1.

The mapping from the lower diagonal of K to GSTIF is now n K_{i,j} --> {GSTIF(LPOINT(I) - I + J)} if NMINVC(I) lt;= J lt;= I. This mapping is performed by a function gpos; if K_{i,j} lies outside the envelope of K which is being stored, gpos returns a value of zero. Otherwise, n K_{i,j} --> {GSTIF(GPOS(LPOINT,NTOTV,I,J))}.

# Equation solution algorithms

As with the algorithms for compact storage of K described above, various equation solution algorithms are used in the main engines' to solve n Gu = f, and will be introduced here, carefully graded from the simple to the more complex. We denote the number of rows and columns by n.

## Gaussian elimination

The simplest matrix solution algorithm is the well-known Gaussian elimination, in which elementary row operations are applied systematically to G -- and simultaneously to the right-hand-side vector -- to get zeroes in the entries below the diagonal; the resulting upper triangular system can then be solved by back-substitution. This algorithm, which does not require that the matrix is symmetric, is used in the main engine' THERM.FOR. The subroutine solve calls the subroutines reduce and backsu in turn, to perform the reduction (to upper triangular form) and back-substitution phases.

This algorithm is represented in the following fragments of pseudo-code. first we show the reduction phase, in which the lower triangular matrix L is formed:

C...Reduction:
do 20 i=1,n-1
pivot = G(i,i)
do 10 k=i+1,n
fact = G(k,i)/pivot
do 5 j=i+1,n
G(k,j) = G(k,j) - fact*G(i,j)
5     continue
f(k) = f(k) - fact*f(i)
10   continue
20 continue

Here, G(i,j) indicates the (i,j)'th entry of G; no modification has been made to represent the compact storage of G. If this pseudo-code is compared with the code in subroutine reduce of THERM.FOR, it will be seen how the mappings to the band storage of G in GSTIF (as described in the previous section for non-symmetric band storage) are implemented. Also, the start and end values in the do loops are modified, because we only need to perform these loops for entries within the band. These factors in translation from the pseudo-code given in the Manual, to the actual code in the main engines', should be borne in mind for when studying the following algorithms.

The solution is completed by the back-substitution phase, performed in subroutine backsu of THERM.FOR:

c...back-substitution phase
u(n) = f(n)/G(n,n)
do 10 i=n-1,1, step -1
pivot = G(i,i)
do 5 j=i+1,n
f(i) = f(i) - G(i,j)*u(j)
5   continue
u(i) = f(i)/pivot
10 continue

## Symmetric Choleski decomposition

We now consider an algorithm closely related to Gaussian elimination, and used in the Basic-level programs POISS and ELAST. A Choleski decomposition of a positive definite matrix K consists of a lower triangular matrix L (with 1's down the diagonal) and an upper triangular matrix U , such that K=LU . If K is symmetric, we modify this to finding a lower triangular matrix L such that K=LL^T . (Now we no longer have 1's down the diagonal of K ) Then our equation to be solved becomes L(L^Tu) = f . Once the decomposition has been performed, we solve the equation Ly=f by forward substitution, to find an intermediate vector y ; this becomes the right-hand-side in the final stage, which is to solve L^Tu=y by back-substitution. This final stage is identical to the back-substitution stage of Gaussian elimination.

The entries in L must be found systematically. The first one which can be found is l_{1,1} , since it will be seen by forming the matrix LL^T and comparing it with K , that k_{1,1} = {l_{1,1}}^2 . Once l_{1,1} is known, we could calculate all the entries below this in the first column. However, we will work in a row-based manner, and find l_{2,1} , use this to find l_{2,2} , before moving on to the 3rd row to find l_{3,1},l_{3,2},l_{3,3} --- then move on to the 4th row, and so on. This row-based approach matches the row-based skyline storage algorithm described earlier.

This Choleski decomposition algorithm is used in the Basic-level programs POISS and ELAST, and can be found in subroutine chodec. The pseudo-code for this decomposition, applied to the matrix G, is as follows:

c                       T
c...decompose G into L.L
i = 1
aii = G(1,1)
clii = sqrt(aii)
G(1,1) = clii
c...decompose each row in turn
do 20 i=2,n
c...work along the row
do 10 j=1,i-1
cljj = G(j,j)
aij = G(i,j)
sum = 0.0
do 5 k=1,j-1
gik = G(i,k)
gjk = G(j,k)
sum = sum + gik*gjk
5     continue
clij = (aij - sum)/cljj
G(i,j) = clij
10   continue
c...finally find the diagonal entry on the row
aii = G(i,i)
sum = 0.0
do 15 k=1,i-1
clik = G(i,k)
sum = sum + clik*clik
15   continue
clii2 = aii - sum
clii = sqrt(clii2)
G(i,i) = clii
20 continue

Note that as each entry in L is found, it is written to G, so that the corresponding entry of the stiffness matrix gets overwritten. At the end of the algorithm, G contains the entries of L in its lower triangle.

In implementing this algorithm in a main engine', we must use a mapping to convert G(i,j) to a storage space in GSTIF, as described in the earlier section on storage techniques.

The system Gu=f can now be rewritten LL^Tu=f , and is solved in two stages:

• solve Ly=f using forward substitution;
• solve L^Tu=y using back-substitution.
The pseudo-code for this, on which subroutine chosub in POISS.FOR is based, is:

c...forward substitution: solve Ly = f
y(1) = f(1)/G(1,1)
do 10 i=2,n
bi = f(i)
sum = 0.0
do 5 k=1,i-1
gik = G(i,k)
ak = y(k)
sum = sum + gik*ak
5   continue
clii = G(i,i)
y(i) = (bi - sum)/clii
10 continue
c                                   T
c...back-substitution phase: solve L x = y
u(n) = y(n)/G(n,n)
do 20 i=n-1,1, step -1
yi = y(i)
sum = 0.0
do 15 k=i+1,n
gki = G(k,i)
ak = u(k)
sum = sum + gki*ak
15   continue
clii = G(i,i)
u(i) = (yi - sum)/clii
20 continue

In program CONSL.FOR, we use the same solution algorithm, but modify its Fortran coding to take account of the skyline storage technique, described above. When row-based skyline storage is used with a Choleski decomposition, the algorithm given above for the final back-substitution phase (solve L^Tu=y ) is rather inefficient. This is because it works with rows of L^T , which are columns of L , the entries of which are scattered throughout the row-based storage array GSTIF. It can easily be modified, however, so that once an entry u_i is found, this is eliminated from all equations from i-1 up to row 1, which involves working up the column of L^T . The modified algorithm, replacing the last phase of the previous pseudo-code, is:

c...back-substitution phase: solve L x = y
c...We can do this in a column-based manner!
do 20 i=n,1, step -1
yi = y(i)
clii = G(i,i)
xi = yi/clii
u(i) = xi
c...once we have found x_i, transfer the term involving this, to
c...the r.h.s. in all earlier equations
if (i.eq.1) goto 20
do 15 j=i-1,1, step -1
clij = G(i,j)
y(j) = y(j) - clij*xi
15   continue
20 continue

This form of the routine is also used in PLADV, where skyline storage is also employed.

## Frontal solution algorithm

The storage of a global stiffness matrix, even using a compact method such as skyline storage, still requires a large amount of space for complex finite element meshes. Researchers have therefore looked for ways of solving Gu=f without having to assemble G fully. One way of doing this is to use iterative solution algorithms, which will be discussed in the next section. However, we first look at a direct solution method which is widely used in commercial f.e. programs. This is called the frontal method. A full description is beyond the scope of this manual (see Chapter 10 for recommended texts), but essentially it is a form of Gaussian elimination in which the element assembly process is combined with the solution. First, the mesh topology is examined to identify the element at which each degree of freedom last appears. Then, entries in element stiffness matrices are read in element-by-element, and assigned space in an array of active variables. Once the last appearance of a d.o.f. has occurred, that variable is eliminated, and the space thus freed up is used for a new variable. Hence, there will be a front' of active variables passing through the matrix. When this elimination phase is completed, the back-substitution phase can occur.

The coding for this algorithm, for the case of a symmetric matrix, can be seen in subroutine front of the advanced elasticity program ELADV.FOR. I am indebted to Dr David Naylor for permission to use this subroutine. A number of global one-dimensional arrays have to be used, but the algorithm is still more space-efficient than skyline storage, and so ELADV should be used for elasticity analyses with complex meshes. There are comment lines in the subroutine, giving some guidance as to how the algorithm works.

The algorithm can be adapted to the case of a non-symmetric matrix G , in which case approximately twice as much storage will be needed. My own adaptation of front for the nonsymmetric case, is the subroutine afront in the elasto-viscoplasticity program VPLAS.FOR. This program handles non-associated plastic flow, which results in non-symmetric stiffness matrices.

The amount of storage needed depends on the numbering of the elements, not the nodes, as the algorithm is dependent on the order in which the elements are assembled into the front. The subroutine makes a preliminary run-through of the elements to check that sufficient storage is provided; this is defined by the parameter MFRON at the start of the program. The entries of G are allocated spaces in the one-dimensional array GSTIF, which is of dimension MFRON. As the assembly process is handled within front, there is no assmb subroutine within ELADV and VPLAS.

## Iterative solution methods

The standard iterative method for solving Gu=f , G symmetric positive definite, is the preconditioned conjugate gradients (PCG) algorithm. This may be summarised as follows:

an {Initialization}
v^{(0)} = 0
r^{(0)} = f
k = 0
z^{(0)} = C^{-1}r^{(0)}
p^{(0)} = z^{(0)}
{Begin iteration}
\alpha_k = \frac{(z^{(k)},r^{(k)})}{(Gp^{(k)},p^{(k)}}
v^{(k+1)} = v^{(k)} + \alpha_k p^{(k)}
r^{(k+1)} = r^{(k)} - \alpha_k Gp^{(k)}
z^{(k+1)} = C^{-1} r^{(k+1)}
ta_k = \frac{(z^{(k+1)},r^{(k+1)})}{(z^{(k)},r^{(k)}}
p^{(k+1)} = z^{(k+1)} + ta_k p^{(k)}
k = k + 1
{Next iteration} a

This form of the PCG algorithm is taken from Bartelt(1989) -- see Chapter 10 for details pf recommended texts. Here, (u,v) denotes the inner product of vectors u and v . The matrix C is the preconditioning matrix. The algorithm is terminated when the norm of the residual vector r^{(k+1)} is less than TOLER times the norm of f ; then the solution u=v^{(k+1)} .

The simplest form of preconditioning is diagonal preconditioning, in which C = {diag}(G) . This algorithm is implemented in the Basic-level program FRAME.FOR. The solution is performed in subroutine pcgsol, the pseudo-code for which is:

c...solve Kv=g using diagk as preconditioner
c...multiplications K.v are performed element-by-element,
c...reading the element matrices of K from the scratch file
c
c...initial guess;
do 2 i=1,n
v(i) = 0.0
r(i) = g(i)
z(i) = 0.0
2 continue
c
c...z = {Kdiag}^-1 * r
do 15 i=1,n
z(i) = r(i)/diagk(i)
15 continue
c
do 16 i=1,n
p(i) = z(i)
16 continue
ztr = 0.0
do 17 i=1,n
zi = z(i)
ri = r(i)
ztr = ztr + zi*ri
17 continue
call norm2(mtotv,ntotv,r,rnorm0)
itscg = 0
write(6,*)'PGSOL: iteration',itscg,',norm = ',rnorm0
c
25 itscg = itscg + 1
c
c...create y which is G*p
call ematgv( p , y )
pty = 0.0
do 28 i=1,n
pty = pty + p(i)*y(i)
28 continue
alpha = ztr/pty
do 30 i=1,n
v(i) = v(i) + alpha*p(i)
30 continue
c
do 40 i=1,n
r(i) = r(i) - alpha*y(i)
40 continue
c
call norm2(mtotv,ntotv,r,rnorm)
c
write(6,*)'       iteration',itscg,',norm = ',rnorm
c
if(rnorm.lt.cgtol) goto 200
c
c...z = {Kdiag}^-1 * r
do 45 i=1,n
z(i) = r(i)/diagk(i)
45 continue
c
ztr0 = ztr
ztr = 0.0
do 60 i=1,ntotv
ztr = ztr + z(i)*r(i)
60 continue
beta = ztr/ztr0
c
do 130 i=1,ntotv
p(i) = z(i) + beta*p(i)
130 continue
c
goto 25
c
200 continue
write(6,*)'Solution converged after ',itscg,' iterations.'

With the diagonal preconditioner (stored in the array diagk), the formation of C^{-1}v is achieved by simply dividing the elements of v by the corresponding elements of diagk. The array diagk is assembled from the element stiffness matrices in the subroutine assdk. The global matrix-vector multiplication Gv is performed element-by-element, reading the element stiffness matrices from the scratch file, in subroutine ematgv; thus, there is no assembly or storage of the global matrix G .

Diagonal preconditioning works surprisingly well, considering its simplicity. If a more sophisticated preconditioner is required, the most popular currently is to perform an Incomplete Choleski (IC) factorization of G . This is essentially a symmetric Choleski decomposition of G into LL^T , as described earlier, but avoiding some or all of the fill-in' which occurs within the envelope of G . That is, we may set L(i,j)=0 instead of L(i,j)=l_{i,j} if G(i,j)=0 . To compensate this, the calculated l_{i,j} is instead added to the diagonal term L(i,i) .

Both diagonal and IC preconditioning are implemented in the advanced plasticity program PLADV.FOR. The form of IC preconditioning there is adaptive: the fill-in is avoided if the magnitude of l_{i,j} is less than a user-specified multiple ( filtol) of the diagonal term. PLADV.FOR uses the skyline storage technique for G . The subroutine inchol in PLADV.FOR performs this IC decomposition. The pseudo-code below should be compared with that for subroutine chodec:

c...Incomplete Choleski decomposition
small = 1.e-8
c...count the amount of fill-in which occurs
kfill = 0
khole = 0
i = 1
a11 = G(1,1)
cl11 = sqrt(a11)
G(1,1) = cl11
c...decompose each row in turn
do 20 i=2,n
aii = G(i,i)
c...we will not fill in if the term is smaller than filtol
filtol = factor*abs(aii)
c...work along the row
do 10 j=1,i-1
cljj = G(j,j)
aij = G(i,j)
sum = 0.0
do 5 k=1,j-1
sum = sum + G(i,k)*G(j,k)
5     continue
clij = (aij - sum)/cljj
c...here is the 'Incomplete' part
gij = G(i,j)
if (gij.eq.0.0) then
khole = khole + 1
if (abs(clij).lt.filtol) then
G(i,i) = G(i,i) + clij
goto 10
endif
c...fill-in will occur
kfill = kfill + 1
endif
G(i,j) = clij
10   continue
c...finally find the diagonal entry on the row
aii = G(i,i)
sum = 0.0
do 15 k=1,i-1
clik = gstif(i,k)
if (clik.eq.0.0) goto 15
sum = sum + clik*clik
15   continue
clii2 = aii - sum
clii = sqrt(clii2)
G(i,i) = clii
20 continue
write(6,*) 'No. of holes found in G matrix = ',khole
write(6,*) 'No. of holes filled-in = ',kfill

As might be expected, there is a trade-off between the amount of fill-in which is prevented, and the effectiveness of the resulting preconditioner. The example datafile pltun4.dat, used with the solution options in PLADV, is discussed in Chapter 9.

Note that this is a very simple form of IC algorithm. There have been many improved algorithms proposed in recent years, and users wishing to apply this technique in practice are advised to search specialist journals for the lastest developments.

The reader will have noticed that additional tolerance parameters must be specified when an iterative solver is used. In PLADV, the iterative solvers are chosen by setting the parameter ISOLA to 2 or 3:

• ISOLA = 0 for direct solution using Choleski ( L.L^T ) decomposition and skyline storage;
• ISOLA = 1 for direct solution using L.D.L^T decomposition and skyline storage;
• ISOLA = 2 for iterative solution using conjugate gradients with diagonal preconditioning;
• ISOLA = 3 for Incomplete Choleski preconditioning with conjugate gradient solution.
If ISOLA = 2 or 3 is chosen, the convergence tolerance cgtol must be input; a value of 10^{-6} is suitable. (An interesting approach is to make the tolerance adaptive: use a greater tolerance in the early stages of the plasticity equilibrium iteration, to avoide unnecessary PCG iterations for a solution which is itself only an approximation, and a smaller tolerance for the final iterations). Then, if ISOLA=3, the fill-in tolerance factor should be input. If factor = 0.0 is chosen, then fill-in will always occur, and so Choleski decomposition will be complete'; in this case, the preconditioner of G will be G itself, and the PCG iteration will converge immediately. See the example of pltun4.dat in the next Chapter, for experience in choosing factor.

# Handling of fixed and specified degrees of freedom

The correct way to handle fixed or specified degrees of freedom in a finite element solution, is to delete those rows from the global system of equations, and in the remaining rows to transfer the terms containing these known values onto the right-hand-sides, before solving for the remaining unknown d.o.f.s. This is however very cumbersome to program, in general. the exception is with the frontal solution method, in which the degrees of freedom are dealt with in turn; it will be found that specified d.o.f.s are treated in this way in front and afront in ELADV.FOR and VPLAS.FOR.

The alternative method for imposing nodal fixities, used in ELAST.FOR, is sometimes called a `big spring' approach. To impose the condition that u_k \approx 0.0 for some k , we add a large number ( 10^{20} ) to the corresponding diagonal term in GSTIF; this is done in subroutine fixdof, called by solve before calling chodec. This will force the value of u_k arising from the solution process to be very small, and we later set such small values to zero.

Fixed values of u_k which are non-zero (e.g. specified displacements, or the specified values of the potential for nodes lying on Dirichlet boundaries) are dealt with similarly. If we want to ensure that u_k = g , we multiply K_{k,k} by 10^{12} , and set the corresponding entry in the r.h.s. vector as F_k = g*K_{k,k} . This process can be seen in fixdof in POISS.FOR.

When a displacement degree of freedom is fixed, there will be a reaction force in this direction at the node, to maintain equilibrium. This can be found after the main solution, by calculating the r.h.s. for the corresponding row. This is done in subroutine reactn in ELAST.FOR and FRAME.FOR; these reactions are written into the load vector ASLOD, which is written in the result file.


text:
Cabal-2.4.0.1: A framework for packaging Haskell software

Distribution.Compat.Directory

Synopsis

# Documentation

listDirectory dir returns a list of all entries in dir without the special entries (. and ..).

The operation may fail with:

• HardwareFault A physical I/O error has occurred. [EIO]
• InvalidArgument The operand is not a valid directory name. [ENAMETOOLONG, ELOOP]
• isDoesNotExistError / NoSuchThing The directory does not exist. [ENOENT, ENOTDIR]
• isPermissionError / PermissionDenied The process has insufficient privileges to perform the operation. [EACCES]
• ResourceExhausted Insufficient resources are available to perform the operation. [EMFILE, ENFILE]
• InappropriateType The path refers to an existing non-directory object. [ENOTDIR]

Since: directory-1.2.5.0

Convert a path into an absolute path. If the given path is relative, the current directory is prepended and then the combined result is normalised. If the path is already absolute, the path is simply normalised. The function preserves the presence or absence of the trailing path separator unless the path refers to the root directory /.

If the path is already absolute, the operation never fails. Otherwise, the operation may fail with the same exceptions as getCurrentDirectory.

Since: directory-1.2.2.0

format: plain


text:
# Substitute font's glyphs with glyphs from the same font

Is it possible to substitute certain font's glyphs with other glyphs from the same font? I am using Xelatex and Fontspec.

More specifically, I am using small caps with Linux Libertine O C and the greek letters' vowels are being typeset with accents. This is a problem as Capitals letters in greek never take accents. I could write the source file without accents, but I want a more general and clean solution (correct spellchecking etc).

What I want is to be able to write Κείμενο and have the same output as Κειμενο, this means that in the following mwe, the first line should typeset as the second one.

Is there, any other way to achieve this behavior instead of remapping the glyphs? I didn't find any relevant option in fontspec's documentation.

I want the first line to typeset as the second one.

\documentclass[]{article}
\usepackage{fontspec}

\begin{document}

\fontspec{Linux Libertine O C}
\textsc{Κείμενο}

\fontspec{Linux Libertine O C}
\textsc{Κειμενο}

\end{document}

In case it is relevant, other fonts like the GFS fonts have the desired (correct) behavior, although I am not sure if their small caps are real or just regular capital letters turned into small caps. Kerkis behaves the same way as Linux Liberine O C

## UPDATE

I contacted Kerkis author, Mr Tsolomitis, and he made (the obvious) suggestion, that the simpler solution is not to use accents in the first place. Either that or use the optional argument in constructs like \section[Κείμενο]{Κειμενο}

As an alternate way (with XeLaTeX) of changing the (default) behavior of Kerkis he suggested changing the *.pfb file with fontforge by erasing the accents of the vowels. Of course the name of font should be changed also.

I also filled a bug report for 'Linux Libertine' but so far no answer.

-
On/Off-Topic: I don't think that the statement "Capital letters in greek never take accents" is exactly right. Why not? Think, for example, poem titles, where accented capitals could make a nice difference. – Nikos Alexandris Dec 2 '11 at 16:31

Yes, it is possible. You can write a mapping and compile it with teckit-compile to get a .tec file which you can then use with the "Mapping" option of fontspec. (like the "tex-text" which maps -- to an en-dash) http://scripts.sil.org/cms/scripts/page.php?site_id=nrsi&id=TECkitDownloads. You can find example of map-files in your tex installation.

But I would suggest at first to report the problem to the font author. It is certainly better if they decide first if this problem can be solved e.g. by an open type feature.

Btw: You are not really consistent when describing your font. Your example uses Linux Libertine C (which doesn't exist), in the text you are using the names Linux Libertine O C (which you meant) and Linux Libertine O (which is another font).

-
Thank you. I corrected my example and I'll report it to the font authors. – pmav99 May 6 '11 at 8:51
After all, I think that Linux Liberine without O does exist. It is the ttf version. And that's why thw MWE compiled on my pc. – pmav99 May 12 '11 at 22:03

After doing some research, I think that the proper way to do it is to create a Glyph Substitution table (GSUB).

FontForge, is up to the task. I haven't tried it, and I don't know how difficult it can be (propably not much - especially for just a few glyphs). I believe that this is the "correct" way of substituting glyphs in OpenType fonts. More info can be found at the FontForge's documentation. Propably the new tables can be easily incorporated by the font authors (if they approve them).

Another more hackish way is to just erase the accents from the vowels or just copy paste the unaccented versions of the vowels using FontForge. This is very fast and it should work out of the box. This is neither clean nor a general solution for glyph substitution.

-
Please note that fontspec supports arbitrary .fea feature files and custom features, e.g. \setmainfont[FeatureFile=greeksmpcp.fea,RawFeature=+gksc;]{Linux Libertine O}. Also, I noticed that version 4.7 was the last release of Linux Libertine to include the “Capitals” variant Linux Libertine O C, which does support Greek small caps, but does not support bold or italic small caps. – Crissov Feb 12 '14 at 16:04

The solution proposed by @UlrikeFischer is xetex-specific; and the one by @pmav99 requires a change to the font, which is difficult and not generic at all.

Moreover, doing some research on the same issue myself (after encountering the same problem) I found that the situation is very inconsistent among the available fonts: as an extreme example, the new Microsoft fonts (eg. Constantia) changed behaviour between versions: older versions produced the "correct" (desired) output, whereas the version shipped with Windows 8 produces output like the one on the picture in the question. I even found a font that does not map ς (final sigma) to anything and retains it in the small-caps form.

So, my suggestion would be to create a macro much like \MakeUppercase, which does appropriate substitutions (Άνθρωπος -> Ανθρωποσ) and then use the small-caps form of the font.

-
That's a great idea! Could you provide a macro definition for something like that? It could even become a package in CTAN. – pmav99 Sep 1 '12 at 14:00
Unfortunately not myself; I do not know almost anything about writing packages. I will probably open a relevant question asking for guidelines... – nplatis Sep 4 '12 at 12:32


text:
Approximation of high-dimensional parametric PDEs

# Approximation of high-dimensional parametric PDEs

Albert Cohen and Ronald DeVore This research was supported by the ONR contracts N00014-11-1-0712 and N00014-12-1-0561, the NSF grant DMS 1222715, the Institut Universitaire de France and the ERC advanced grant BREAD.
###### Abstract

Parametrized families of PDEs arise in various contexts such as inverse problems, control and optimization, risk assessment, and uncertainty quantification. In most of these applications, the number of parameters is large or perhaps even infinite. Thus, the development of numerical methods for these parametric problems is faced with the possible curse of dimensionality. This article is directed at (i) identifying and understanding which properties of parametric equations allow one to avoid this curse and (ii) developing and analyzing effective numerical methodd which fully exploit these properties and, in turn, are immune to the growth in dimensionality.

The first part of this article studies the smoothness and approximability of the solution map, that is, the map where is the parameter value and is the corresponding solution to the PDE. It is shown that for many relevant parametric PDEs, the parametric smoothness of this map is typically holomorphic and also highly anisotropic in that the relevant parameters are of widely varying importance in describing the solution. These two properties are then exploited to establish convergence rates of -term approximations to the solution map for which each term is separable in the parametric and physical variables. These results reveal that, at least on a theoretical level, the solution map can be well approximated by discretizations of moderate complexity, thereby showing how the curse of dimensionality is broken. This theoretical analysis is carried out through concepts of approximation theory such as best -term approximation, sparsity, and -widths. These notions determine a priori the best possible performance of numerical methods and thus serve as a benchmark for concrete algorithms.

The second part of this article turns to the development of numerical algorithms based on the theoretically established sparse separable approximations. The numerical methods studied fall into two general categories. The first uses polynomial expansions in terms of the parameters to approximate the solution map. The second one searches for suitable low dimensional spaces for simultaneously approximating all members of the parametric family. The numerical implementation of these approaches is carried out through adaptive and greedy algorithms. An a priori analysis of the performance of these algorithms establishes how well they meet the theoretical benchmarks.

## 1 Overview

### 1.1 Parametric and stochastic PDEs

Partial differential equations (PDEs) are commonly used to model complex systems in a variety of physical contexts. When solving a given PDE, one typically fixes certain parameters: the shape of the physical domain, the diffusion or velocity field, the source term, the flux or reaction law, etc. We use the terminology parametric PDEs when some of these parameters are allowed to vary over a certain range of interest. When treating such parametric PDEs, one is interested in finding the solution for all parameters in the range of interest.

To describe such problems in their full generality, we adopt the formulation

 P(u,a)=0, (1.1)

where denotes the parameters, is the unknown of the problem, and

 P:V×X→W, (1.2)

is a linear or nonlinear partial differential operator, with a triplet of Banach spaces. We assume that the parameter ranges over a compact set , and that for any in this range there exists a unique solution to (1.1). This allows us to define the solution map

 u:a↦u(a), (1.3)

which acts from onto and is well defined over . We also define the solution manifold as the family

 M=u(A)={u(a):a∈A}, (1.4)

which gathers together all solutions as the parameter varies within its range.

One simple guiding example, which will be often used throughout this article, is the linear elliptic equation

 −div(a∇u) = f,on D, u = 0on ∂D, (1.5)

set on a Lipschitz domain . Here, we fix the right side as a real valued function and consider real valued diffusion coefficients as the parameter. The corresponding operator is therefore given by

 P(u,a)=f+div(a∇u). (1.6)

A possible choice for the triplet of spaces is then

 (X,V,W)=(L∞(D),H10(D),H−1(D)). (1.7)

Indeed, if , and , one then defines as an element of by requiring that

 ⟨P(u,a),v⟩=⟨f,v⟩+∫Da∇u⋅∇v,v∈V, (1.8)

where is the duality bracket between and . Lax-Milgram theory ensures the existence and uniqueness of a solution to (1.1) from , if for some , the diffusion satisfies the ellipticity condition

 a(x)≥r,x∈D. (1.9)

Therefore, a typical parameter range is a set , which in addition is assumed to be compact in .

Although elementary, the above example gathers important features that are present in other relevant examples of parametric PDEs. In particular, the solution map acts from an infinite dimensional space into another infinite dimensional space. Also note that, while the operator of (1.6) is linear both in and (up to the constant additive term ) the solution map is nonlinear. Because of the high dimensionality of the parameter space , such problems represent a significant challenge when trying to capture this map numerically. One objective of this article is to understand which properties of this map allow for a successful numerical treatment. Concepts such as holomorphy, sparsity, and adaptivity are at the heart of our development.

The solution map may also be viewed as a function

 (x,a)↦u(x,a) (1.10)

of both the physical variable and the parametric variable . The parametric variable has a particular status because its different instances are uncoupled: for any fixed instance , we may solve the PDE exactly or approximately and therefore compute for all values of , while ignoring all other values . This plays an important role for certain numerical methods which are based on solving the parametric PDE for different particular values of , since this task can be parallelized.

Parametric PDEs occur in a variety of modeling contexts. We draw the following major distinctions in their setting:

• Deterministic modeling: the parameters are deterministic design or control variables which may be tuned by the user so that the solution , or some quantity of interest , has prescribed properties. For instance, if the elliptic equation (1.5) is used to model the heat conduction in a component produced by an industrial process, one may want to design the material in order to minimize the heat flux on a certain part on the boundary , in which case the quantity of interest to be optimized is

 Q(u)=∫Γa∇u⋅n, (1.11)

where is the outer normal. This amounts to optimizing the function

 a↦F(a):=Q(u(a)). (1.12)

over .

• Stochastic modeling: the parameters are random variables with prescribed probability laws, which account for uncertainties in the model. Therefore has a certain probability distribution supported on . One is then typically interested in the resulting probabilistic properties of the solution , which is itself a random variable over with values in , or in the probabilistic properties of a quantity of interest . For instance, if the elliptic equation (1.5) is used to model oil or ground water diffusion, a common way to deal with the uncertainty of the underground porous media is to define as a random field with some prescribed law. Then one might want to estimate the mean solution

 ¯¯¯u:=E(u), (1.13)

or its variance

 V(u):=E(∥u−¯¯¯u∥2V), (1.14)

or the average flux through a certain interface , that is, with as in (1.11), or the probability that this flux exceeds a certain quantity, etc.

In both the deterministic and stochastic settings, the given application may require evaluating for a large number of instances of the parameter . This is the case, for instance, when using a descent method for optimizing a quantity of interest in the deterministic framework, or a Monte-Carlo method for evaluating an expectation in the stochastic framework. Since each individual instance is the solution of a PDE, its exact evaluation is typically out of reach. Instead, each query for is approximately evaluated by a numerical solver, which may itself be computationally intensive for high accuracy.

In order to significantly reduce the number of computations needed for attaining a prescribed accuracy, alternate strategies, commonly referred to as reduced modeling, have been developed. Understanding which of these strategies are effective, in the case where the parameter has large or infinite dimension, and why, is the subject of this article.

### 1.2 Affine representation of the parameters

So far our description of the set of parameters allows it to be any compact subset of . An important ingredient in both our theoretical and numerical developments is to identify any through a sequence of real numbers. We are especially interested in affine representations of . We say that a sequence of functions is an affine representer for , or representer for short, if we can write each as

 a=a(y)=¯¯¯a+∑j≥1yjψj,y:=(yj)j≥1,yj=yj(a), (1.15)

where the are real numbers, is a fixed function from , and the series converges in the norm of for each . We are making a slight abuse of notation here since we use to represent a general element of and also use to represent the map

 a:y↦a(y), (1.16)

from to . But the meaning will always be clear from the context.

It is easy to see that for any compact set in a Banach space affine representers exist. For example, if has a Schauder basis then any such basis will be a representer. Even if does not admit a Schauder basis, as is the case for our example , we can still find representers as follows. Choose any . Since is compact there exist finite dimensional spaces , with , such that

 dist(K,Xn)X:=supa∈Kminb∈Xn∥a−b∥X→0asn→∞. (1.17)

We can also take the spaces to be nested, that is

 Xn⊂Xn+1,n≥0. (1.18)

Let be any basis for and define . The sequence

 ψj:=ϕj−Nn,n,j=Nn,…,Nn+1, (1.19)

contains each of the bases for all . Given any , let be a best approximation in to from ,with . Then, we can write

 a=¯a+∞∑n=1(an−an−1). (1.20)

Each term is in and hence can be written as a linear combination of the . Therefore, is a representer for .

Affine representations (1.15) often occur in the natural formulation of the parametric problem. For instance, if the diffusion coefficient in (1.5) is piecewise constant over a fixed partition of the physical domain , then it is natural to set

 a(y)=¯¯¯a+d∑j=1yj\largeχDj, (1.21)

where is a constant and the are the characteristic functions of the subdomains . Similarly, if the parameter describes the shape of the boundary of the physical domain in a computer-aided design setting, a typical format is

 a(y)=¯¯¯a+d∑j=1yjBj, (1.22)

where represents a nominal shape and the are B-spline functions associated to control points. In these two examples, is finite, yet possibly very large.

In the statistical context, if is a second order random field over a domain , a frequently used choice in (1.15) is , the average field, and , the Karhunen-Loeve basis, that is, the eigenfunctions of the covariance operator

 v↦Rav:=∫DCa(⋅,x)v(x)dx,Ca(z,x):=E((a(z)−¯¯¯a(z))(a(x)−¯¯¯a(x)). (1.23)

Then, the resulting scalar variables are centered and uncorrelated, that is, and when .

Even if an affine representation of the form (1.15) is not given in the formulation of the problem, one can be derived by taking any representation system in the Banach space . For example, if admits a Schauder basis, then one can take any such basis for and arrive at such an expansion. In classical spaces , such as or Sobolev spaces, standard systems of approximation, such as Fourier series, splines, or wavelets can be used.

The advantage of the representation (1.15) is that can now be identified through the sequence . When considering all , we obtain a family of such sequences. Note that this family can be quite complicated. In order to simplify matters, we normalize the , so that for any ,

 supa∈A|yj(a)|=1. (1.24)

Such a renormalization is usually possible because is compact and depends continuously on . After this normalization, for each , the sequence belongs to the infinite dimensional cube

 U:=[−1,1]N. (1.25)

Notice that taking a general sequence from this cube, there may not be an with , . Also, if is not a basis , the representation (1.15) may not be unique. We define

 UA:={(yj)j≥1∈U:∑j≥1yjψj∈A}. (1.26)

We are mainly interested in representers , for which

 ¯a+∑j≥1yjψj (1.27)

converges in for each . We call such representers complete. In this case, we may define

 a(U):={a=a(y)=¯¯¯a+∑j≥1yjψj:(yj)j≥1∈U}, (1.28)

so that

 A⊂a(U). (1.29)

A typical case of a complete representer is when is a sequence in .

Once an affine representation has been chosen, the initial solution map becomes equivalent to the map which is defined on . With an abuse of notation, we write this new solution map as

 y↦u(y):=u(a(y)). (1.30)

This is a Banach space valued function of an infinite number of variables. Note that in the case where the affine representation has a finite number of terms, the range of is . However, the infinite dimensional case subsumes the finite dimensional case, since the latter may be viewed as a particular case with for .

In the case of a complete representer, is defined on all of . However, we do not know whether the solution map is defined on all of . To guarantee this, the following assumption will be used often.

Assumption A: The parameter set has a complete representer and the solution map is well defined on the whole set , or equivalently the solution map is well defined on the whole set .

This assumption naturally holds when the set is exactly defined as .

### 1.3 Smoothness of the solution map

One objective of this article is to develop efficient numerical approximations to the solution maps of (1.3) or (1.30). One of the main difficulties is that these maps are high or infinite dimensional, in the sense that the dimension of the variable or is high or infinite. In order to understand what might be good strategies for constructing such approximations, we need first to understand the inherent properties of these maps that might allow us to circumvent this difficulty.

We initiate such a program in §2, where we first analyze the smoothness of the solution map . In the case of the elliptic equation (1.5), it is easily seen that this map is not only infinitely differentiable, but also admits a holomorphic extension to certain subdomains of the complex valued . We propose two general approaches which allow us to establish similar holomorphy properties for other relevant instances of linear and nonlinear parametric PDEs. One first approach is based on the Ladyzenskaia-Babushka-Brezzi theory. It applies to a range of linear PDEs where the operator and the right hand side have holomorphic dependence in . These include parabolic and saddle-point problems, such as the heat equations or the Stokes problem, with parameter in the diffusion term, similar to (1.5). One second approach is based on the implicit function theorem in complex valued Banach spaces. In contrast to the first approach, it can be applied to certain nonlinear PDEs.

Using the affine representation (1.15) of , we then study the solution map under Assumption A, which means that it is defined on the whole of . In addition to holomorphy, an important property of can be extracted from the affine representation (1.15). The functions appearing in (1.15) have norms of varying size. Since the variable is scaled to be in , when is small, this variable has a reduced effect on the variations of . Thus the variables are not democratic, but rather they have varying importance. In other words, the map is highly anisotropic. More specifically, we derive holomorphic extension results for this map on certain multivariate complex domains of tensor product type. In particular, we consider polydiscs of the general form

 Uρ:=⊗{|zj|≤ρj}={z=(zj)j≥1∈CN:|zj|≤ρj} (1.31)

where is a positive sequence which serves to describe the anisotropy of the solution map. We also consider polyellipses which deviate less far from the real axis. These holomorphy domains play a key role in the derivation of approximation results.

###### Remark 1.1

While we are generally interested in real valued solutions to the parametric PDE (1.1), corresponding to real valued parameters or , our analysis of holomorphic smoothness leads us naturally to complex valued solutions, corresponding to complex valued parameters. For this reason, the spaces are always assumed to be complex valued Banach spaces throughout this paper.

### 1.4 Approximation of the solution map

Reduced modeling methods seek to take advantage of the properties of the solution maps or such as the holomorphy and anisotropy mentioned above. These properties suggest strategies for approximating these map by simple functions in which the physical variables and the parametric variable or are separated and hence take the form

 (x,a)↦un(x,a):=n∑i=1vi(x)ϕi(a), (1.32)

or

 (x,y)↦un(x,y):=n∑i=1vi(x)ϕi(y), (1.33)

where are functions of living in the solution space and are functions of or with values in or .

We may view as a rank approximation to , in analogy with low rank approximation of matrices. We adopt the notations

 a↦un(a)=un(⋅,a)andy↦un(y)=un(⋅,y), (1.34)

for the above approximations.

Let us discuss the potential accuracy of separable approximations of the form (1.32). If our objective is to capture for all with a prescribed accuracy , this means that we search for an error bound in the uniform sense, i.e., of the form

 ∥u−un∥L∞(A,V):=supa∈A∥u(a)−un(a)∥V≤ε(n). (1.35)

For certain applications, in particular in the stochastic framework, we may instead decide to measure the error on average, for instance by searching for an error bound in the mean square sense,

 E(∥u−un∥2V):=∥u−un∥2L2(A,V,μ)=∫A∥u(a)−un(a)∥2Vdμ(a)≤ε(n)2, (1.36)

where is the probability measure for the distribution of over . Since is a probability measure, one has for any ,

 ∥v∥L2(A,V,μ)≤∥v∥L∞(A,V). (1.37)

Therefore the uniform bound is stronger than the average bound, in the sense that (1.35) implies (1.36) with the same value of .

Likewise, for the approximation of the map , we may search for a uniform bound

 ∥u−un∥L∞(UA,V):=supy∈UA∥u(y)−un(y)∥V≤ε(n) (1.38)

or a mean square bound

 E(∥u−un∥2V):=∥u−un∥2L2(UA,V,μ)=∫UA∥u(y)−un(y)∥2Vdμ(y)≤ε(n)2, (1.39)

where is the probability measure for the distribution of over .

###### Remark 1.2

We do not indicate the measure in our notation or , since in all relevant examples considered in this article we always consider the exact supremum over in or over in , rather than the essential supremum.

For any , the approximation belongs to

 Vn:=span{v1,…,vn}, (1.40)

which is a fixed -dimensional subspace of . Ideal benchmarks for the performance of separable expansions of the form (1.32) may thus be defined by selecting optimal -dimensional spaces for the approximation of in either a uniform or an average sense.

For uniform error bounds, this benchmark is given by the concept of Kolmogorov’s -width, which is well known in approximation theory. If is a compact set in a Banach space , we define its Kolmogorov -width as

 dn(K)V:=infdim(Vn)≤nsupv∈ÊKminw∈Vn∥v−w∥V. (1.41)

This quantity, first introduced in [57], describes the best achievable accuracy, in the norm of , when approximating all possible elements of by elements from a linear -dimensional space . Obviously, the optimal choice of for approximation of in a uniform sense corresponds to the space, if it exists, that reaches the above infimum when is taken to be the solution manifold . The best achievable error in the uniform sense is thus given by

 ε(n):=dn(M)V. (1.42)

There exist many other notions of widths that are used to measure the size of compact sets. Here we only work with the above one and refer to [74] for a more general treatment.

For mean square bounds, in the case where is a Hilbert space, the corresponding benchmark is related to the concept of principal component analysis, which is of common use in statistics. By choosing an arbitrary orthonormal basis of , we may expand according to

 u(a)=∑i≥0ziei,zi=zi(a):=⟨u(a),ei⟩V. (1.43)

Approximation of from an -dimensional space of is then equivalent to the approximation of from an -dimensional space of . The optimal space is then obtained through the study of the correlation operator

 R=(Ri,j)i,j≥0,Ri,j:=E(zi(a)zj(a))=∫a∈Azi(a)zj(a)dμ(a). (1.44)

This operator is symmetric, positive, compact and of trace class. It therefore admits an orthonormal basis of eigenvectors associated to a positive, non-increasing and summable sequence of eigenvalues. The space

 Gn:=span{g1,…,gn}, (1.45)

minimizes over all -dimensional spaces the mean square error between and its orthogonal projection , with

 E(∥z−PGnz∥2ℓ2)=∑k>nλk. (1.46)

In turn, the optimal space is spanned by the functions

 vk:=∑i≥0gk,iei,k=1,…,n, (1.47)

where is the -th component of , that is, . The best achievable error in the mean square sense is thus given by

 ε(n)2:=∑k>nλk. (1.48)

Note that, in view of (1.37), one has the comparison

 ∑k>nλk≤dn(M)2V. (1.49)

The above described optimal spaces are usually out of reach, both from an analytic and computational point of view, and it is therefore interesting to consider sub-optimal approximations. In addition, when considering the solution map , the tensor product structure of allows us to consider approximations based on further separation between the parametric variables, that is, where each function in (1.33) is itself a product of univariate functions of the different . The simplest example of such approximations are multivariate polynomials, which have the general form

 un(x,y)=∑ν∈Λnvν(x)yν,yν:=∏j≥1yνjj, (1.50)

where each index is a finitely supported sequence of positive integers, or equivalently such that , and is a set of such sequences.

In §3, we obtain such polynomial approximations by taking finite portions of infinite polynomial expansions of . Here, we work under Assumption A, which means that is defined on the whole of . We consider two types of expansions:

• Power series the form

 u(y)=∑ν∈Ftνyν, (1.51)

where is the set of all finitely supported sequence of positive integers.

• Orthogonal series of the form

 u(y)=∑ν∈FwνPν(y),Pν(y):=∏j≥1Pνj(yj), (1.52)

where is the Legendre polynomial of degree defined on .

Using the holomorphy and anisotropy properties of the solution map established in §2 for specific classes of parametric PDEs, we derive a priori bounds on the -norms and of the coefficients which appear in these expansions. In this way, we are able to establish algebraic convergence rates for certain truncations of the above expansions, where is the cardinality of the truncation set .

One critical aspect of the truncation strategy is that we retain the largest coefficients, which is a form of nonlinear approximation, also known as sparse or best -term approximation. With such a choice for , the exponent in the convergence rate is related to the available summability for of the -norms of the coefficients in the considered infinite expansion. In particular, for uniform approximation estimates, that is, in , one has

 s=1p−1, (1.53)

once the -summability of these -norms has been proven. The main result from §3 shows that, under suitable assumptions, the summability of the sequence implies that the norms of the coefficients in the expansions (1.51) or (1.52) are also summable.

The fact that we obtain the algebraic convergence rate despite the infinite dimensional nature of the variable reveals that the curse of dimensionality can be avoided in the approximation of relevant parametric PDEs.

### 1.5 The n-widths of the solution manifold M

In both the uniform or mean square ways of measuring error, as described above, the success of reduced modeling can be proven if converges to zero sufficiently fast as . Thus, the study of these widths constitutes a major subject in the theoretical justification of reduced modeling.

A common way to measure the widths of compact sets in classical spaces is to embed into an appropriate smoothness space such as a Sobolev or Besov space. For example, in our model parametric elliptic equation (1.5), the space is , and this approach would lead us to examine the Sobolev smoothness of the individual functions for . Classical elliptic regularity theory says that for smooth domains the smoothness of can be inferred from the smoothness of the right side . However, for general domains, there are severe limits on this regularity due to the irregularity of the boundary. Therefore, bounding the decay of widths of through such regularity results will generally prove only slow decay for and therefore is not useful for obtaining the fast decay rates we seek. Indeed, recall, that classical smoothness spaces, such as Sobolev or Besov spaces of order in variables, have widths that decay, at best, like as when these widths are measured in an norm for some . For example, it is known that if is the unit ball of , then its -width in satisfies

 cdn−m/d≤dn(K)L∞≤Cdn−m/d,n≥1. (1.54)

Likewise, if is the unit ball of , then its -width in satisfies

 cdn−(m−1)/d≤dn(K)H1≤Cdn−(m−1)/d,n≥1. (1.55)

Also note that the regularity of , as measured by membership in Sobolev and Besov spaces, is closely related to the performance of piecewise polynomial approximation such as that used in finite element methods and is the reason why these algorithms have rather slow convergence. So the fast decay of to zero cannot be obtained by such an approach.

The widths go to zero fast not because the individual elements in are smooth in the physical variable, but rather because is the image of the solution map which is, as previously discussed, smooth and anisotropic in the parametric variable. However, let us remark that it is by no means trivial to deduce fast decay for from this fact alone, because of what is called the curse of dimensionality. Namely, approximation rates for a given target function are generally proved by showing that the target function, in our case the function , has sufficiently high regularity. But, in high dimensions, regularity by itself is usually not enough. Indeed, returning to the bounds (1.55), we see that the large dimension affects the approximation rate in two detrimental ways. The exponent in the smoothness rate is divided by and the constants are known to grow exponentially with increasing . In our case can even be infinite and this makes the derivation of approximation rates a subtle problem.

In §4, we discuss general principles for estimating by above the -widths of the solution manifold . One immediate consequence of the approximation results in §3 is that, for specific classes of parametric PDEs, when Assumption AÊ holds and , then these widths decay at least like with given by (1.53). We extend this analysis to the general case when Assumption AÊ does not necessarily hold. Since the manifold is not directly accessible, one would like to understand what properties of , which is assumed to be completely known to us, imply decay rates for . We show that the asymptotic decay of the -width of in is related to that of the -width of in . This follows from general results on widths of images of compact sets under holomorphic maps. Namely, if is holomorphic in a neighborhood of a general compact set , then we prove the following result on the -width of the image in :

 supn≥1nrdn(A)X<∞⇒supn≥1nsdn(M)V<∞,s

This result shows that from the view point of preserving -widths, holomorphic maps behave almost as good as linear maps, up a loss of in the convergence rate. One open problem is to understand if this loss is sharp or could be improved.

### 1.6 Numerical methods for reduced modeling

The above mentioned results on holomorphic extensions, sparse expansions for , and -widths of the manifold , can be thought of as theoretical justifications for the role of reduced modeling in solving parametric and stochastic equations. They provide evidence that reduced modeling numerical methods should yield significant computational savings over traditional methods such as finite element solvers for parametric problems or Monte Carlo methods for stochastic problems. However, they do not constitute actual numerical methods.

The second part of our article turns to the construction of numerical algorithms motivated by these theoretical results. Such algorithms compute specific separable approximations of the form (1.32) or (1.33) for a given value of at an affordable computational cost. Our objective, in this regard, is not to give an exhaustive description of numerical reduced modeling methods and their numerous variants. Rather, our main focus is to introduce some important representative examples of these methods for which an a priori analysis quantifies the gains in numerical performance of these methods.

One important distinction is between non-adaptive and adaptive methods. In the first ones, the choice of the functions or of used in (1.32) or (1.33), for a given value of is made in an a priori manner, if possible using the available information on the problem. In the second ones, the computations executed for lower values of are exploited in order to monitor the choice at stage . One desirable feature of an adaptive algorithm is that it should be incremental or greedy: only one new function or is added at each stage to the previously selected functions which are left unchanged. Adaptive algorithms are known to often perform better than their non-adaptive counterpart, but their convergence analysis is usually more delicate.

A second important distinction between the various numerical methods is whether they are non-intrusiveÊ or intrusive. A non-intrusive algorithm builds on an existing exact or approximate solver for the PDE which may be computationally expensive. It derives approximations of the form (1.32) or (1.33) by choosing instances or and using the values or computed by the solver. Non-intrusive algorithms may be implemented even when this solver is a black box, and therefore with a possibly limited knowledge on the exact PDE model. An intrusive algorithm, on the other hand, directly exploits the precise form of the PDE for computing the approximation (1.32) or (1.33), and therefore requires the full knowledge of the PDE model for its implementation.

It should be noted that instances as well as the functions used in (1.32) or (1.33) can only be computed with a certain level of spatial discretization, for example resulting from a finite element solver. In such a case they belong to a finite dimensional space . If the same finite element space is used to discretize all instances up to a precision that is satisfactory for the user, this means that we are actually trying to capture the approximate solution maps

 a↦uh(a)∈Vh, (1.57)

or

 y↦uh(y)∈Vh. (1.58)

The analysis of the performance of numerical reduced modeling methods needs to incorporate the additional error produced by this discretization.

### 1.7 Sparse polynomial approximation algorithms

One first class of methods that we analyze consists in finding a numerically computable polynomial approximation of the form (1.50). There are two major tasks in constructing such a numerical approximation: (i) find good truncation sets and (ii) numerically compute an approximation to the coefficients for each . By far, the most significant issue in numerical methods based on polynomial expansions is to find a good choice for the sets . If everything was known to us, we would simply take for the set of indices corresponding to the largest -norms of the coefficients in (1.51) or (1.52). However, finding such an optimal would require in principle that we compute the coefficients for all values which is obviously out of reach. In addition, the structure of the optimal set can be quite complicated. One saving factor is that our analysis in §3 gives a priori bounds on the size of these coefficients. These bounds can be used in order to make an a priori selection of the sets . This is a non-adaptive approach, which generally gives suboptimal performance due to the possible lack of sharpness in the a priori bounds. This leads one to try to enhance performance by combining the a priori bounds together with an adaptive selection of the sets .

The numerical methods are facilitated by imposing that the selected index sets are downward closed (or lower sets), i.e. satisfy the following property:

 ν∈Λnand~ν≤ν⇒~ν∈Λn, (1.59)

where means that for all .

In §6, we discuss algorithms which compute the polynomial approximation by an interpolation process. These algorithms are non-intrusive and apply to a broad scope of problems. One key issue is the choice of the interpolation points, which is facilitated by the following result: if is a sequence of pairwise distinct points in and if , then for any downward closed set , any polynomial of the form (1.50) is uniquely characterized by its value on the grid . This allows us to construct the interpolation in a hierarchical manner, by simultaneously incrementing the polynomial space and the interpolation grid. The sets can either be a priori chosen based on the bounds on the coefficients established in §3 or adaptively generated. These sets generally differ from the ideal sets corresponding to the largest coefficients (which may not fullfill the downward closedness property). We show that certain choices of the univariate sequence , known as Leja or -Leja points lead to stable interpolation processes (in the sense that the Lebesgue constant have moderate growth with the number of points) allowing us to retrieve by interpolation the same algebraic convergence rate which are proved for the best polynomial approximations.

In §7, we discuss another class of algorithms which recursively compute the exact coefficients in the Taylor series of the approximate solution map (1.58). These algorithms are intrusive and they only apply to problems where is linear both in and up to a constant term, such as the elliptic equation (1.5) which serves a guiding example. The recursive computation is facilitated by imposing that the index sets in the truncated Taylor expansion are downward closed. Similar to the interpolation algorithm from §6, the sets can either be a priori chosen based on the available bounds on the coefficients established in §3, or adaptively generated. One main result shows that adaptive algorithms based on a so-called bulk chasing procedure have the same convergence rate as the one which is established when using the index sets corresponding to the largest Taylor coefficients.

### 1.8 Reduced basis algorithms

The second class of methods that we analyze seeks to find, in an offline stage, a set of functions for which the resulting -dimensional space is close to an optimal linear -dimensional approximation space. For mean square estimates, one such approach, known as proper orthogonal decomposition, consists in building the functions based on an approximation of the exact covariance operator (1.44) computed from a sufficiently dense sampling of the random solution . Another approach, which targets uniform estimates, is the reduced basis method, which consists in generating by a selection of particular solution instances chosen from a very large set of potential candidates. In both cases, the offline stage is potentially very costly.

Once such a space is chosen, one builds an online solver, such that for any given , the approximate solution is an element from . There are several possibilities on how to build this online solver. The most prominent of these is to take the Galerkin projection of onto which consist of finding by solving the system of equations

 ⟨P(un(a),a),w⟩=0,w∈Vn (1.60)

for a suitable duality product . This online computation determines for each the values for which appears in (1.32). The advantage of using the Galerkin solver, is that, for certain problems such as elliptic ones, it is known to give the best error in approximating by elements of when error is measured in the norm of . Its disadvantage its computational cost in finding given the query . For this reason other projections onto are also studied, some of these based on interpolation.

The key issue when using reduced basis methods is how to find a good space , i.e. how to find good basis functions. In §8, we discuss an elementary greedy strategy for the offline selection of the instances , that consists in picking the -th instance which deviates the most from the space generated from the previously selected ones. The approximation error

 σn(M):=supv∈Minfw∈Vn∥v−w∥, (1.61)

produced by such spaces may be significantly larger than the ideal benchmark of the -width of the solution manifold for a given value of . However, a striking result is that both are comparable in terms of rate of decay: for any , there is a constant such that

 supn≥1nsσn(M)≤Cssupn≥1nsdn(M)V. (1.62)

Similar results are established for exponential convergence rates.

While both classes of numerical methods aim to construct separable approximations of the form (1.32) or (1.33), there is a significant distinction between them in the way they organize computation. For the first class of polynomial approximation methods, the offline stage fixes the polynomial functions through the selection of the set and computes the coefficients . Then the online stage is in some sense trivial since it simply computes through the linear combination (1.50). For the second class, the online stage still requires solving PDE approximately in the chosen reduced space . This offline/online splitting makes it difficult to draw a fair comparison between the different methods from the point of view of computational time vs accuracy.

Notation for constants: Numerous multiplicative constants appear throughout this paper, for example in convergence estimates. We use the generic notation , which may therefore change value between different formulas, and if necessary we indicate the parameters on which depends. We use a more specific notation if we want to express the dependence of the constant with respect to a cetain parameter (for example the dimension in (1.55)) or if we want to refer to this specific constant later in the paper.

### 1.9 Historical orientation

Numerical methods for parametric and stochastic PDEs using polynomials (or other approximation tools) in the parametric variable have been widely studied since the 1990’s. We refer in particular to [40, 46, 55, 56, 90] for general introductions to these approaches, and to [1, 41, 37, 53, 54, 83, 84] for related work prior to the results exposed in our paper.

The approximation results presented in §3 have been obtained by the authors and their co-authors in a series of paper [23, 24, 19], and the results in §4 on the evaluation of -widths are from [25]. These results establish for the first time convergence rates immune to the curse of dimensionality, in the sense that they hold with infinitely many variables, see also [44] for a survey dealing in particular with these issues. In a similar infinite dimensional framework, and not covered in our paper, let us mention the following related works: (i) similar holomorphy and approximation results are established in [47, 48, 58] for specific type of PDEs and control problems, (ii) approximation of integrals by quadratures is discussed in [60, 61], (iii) inverse problems are discussed in [78, 79, 82], following the Bayesian perspective from [87], and (iv) diffusion problems with lognormal coefficients are treated in [50, 43, 45].

The sparse interpolation method presented in §6 is introduced and studied in [18]. See also [2, 4, 71, 72] for related work on collocation methods. Other non-intrusive methods are based on least-square regression as discussed in [16, 34, 36, 69], or on pseudo-spectral approximation as discussed in [91, 26]. The Taylor approximation algorithm presented in §7 is introduced and studied in [17]. Other intrusive methods based on Galerkin projection are discussed in [4, 5, 23, 42].

Reduced basis methods have been studied since the 1980’s [73]. The greedy algorithms presented in §8 have been introduced and discussed in [89, 66, 67, 68, 88], and their convergence analysis was given in [6] and [32]. We refer to [52] for a general introduction on the related POD method, which is not discussed in our paper.

## 2 Holomorphic extensions

In this section, we discuss smoothness properties of the solution maps and which are central to the development of efficient numerical methods that are immune to the curse of dimensionality. We show that, under suitable assumptions on the parametric PDE, these maps admit holomorphic extensions to certain complex domains. Recall that a map from a complex Banach space to a second complex Banach space is said to be holomorphic on an open set if for each , has a Frechet derivative at . Here is a linear operator mapping to such that

 ∥F(x+h)−F(x)−dF(x)h∥Y=o(∥h∥X),h∈X. (2.1)

### 2.1 Extension of a↦u(a) for the model elliptic equation

In order to formulate results on holomorphy, we need to introduce existence-uniqueness theory for solutions to (1.1) in the case that is complex valued.

We begin with our guiding example of the elliptic equation (1.5). We now consider

 (X,V,W)=(L∞(D),H10(D),H−1(D)) (2.2)

as spaces of complex valued functions, and extend the standard variational formulation to such spaces in a straightforward manner: for a given , and with , find such that

 ∫Da∇u⋅∇v=⟨f,v⟩W,V,v∈V, (2.3)

where in the left integrand,

 ∇u⋅∇v:=m∑i=1∂xiu¯¯¯¯¯¯¯¯¯∂xiv, (2.4)

is the standard hilbertian inner product, and is the anti-duality pairing between and , which when is given by the hilbertian inner product

 ⟨f,v⟩=∫Df¯¯¯v. (2.5)

We recall that

 ∥v∥V:=∥∇v∥L2(D), (2.6)

and

 ∥v∥W:=sup{⟨v,w⟩:∥w∥V≤1}. (2.7)

This is therefore a particular case of a linear problem with the following general variational formulation. Let denote the set of all sesquilinear forms defined on and let be the set of all antilinear functionals defined on , i.e., is the antidual of . We define the following norm on :

 ∥B∥:=sup∥v∥V≤1, ∥w∥V≤1|B(v,w)|. (2.8)

Problem: Given and , find such that

 B(u,v)=L(v),∀v∈V. (2.9)

The existence-uniqueness theory for such problems can be proven from the complex version of Lax-Milgram theorem given in Theorem 2.1. This theorem is a particular case of Theorem 2.2 proved below. To formulate this theorem and for later use, we introduce the notation for the space of all linear operators mapping the Banach space into the Banach space with its usual norm

 ∥T∥L(X,Y):=supx∈X∥Tx∥Y∥x∥X. (2.10)

Given , one has that is an anti-linear functional and hence for any , there is a such that

 B(u,v)=⟨Bu,v⟩W,V,v∈V, (2.11)

where is the anti-duality pairing between and . Therefore, is a linear operator from into and its norm is the same as that of :

 ∥B∥L(V,W)=∥B∥. (2.12)

So the operator is bounded and hence continuous. The problem (2.9) is equivalent to the equation

 Bu=L, (2.13)

set in . With this notation and remarks in hand, we can now state the complex version of the Lax-Milgram theorem.

###### Theorem 2.1

Assume that, is a sesquilinear form on such that

 |B(u,u)|≥α∥u∥2V,u∈V, (2.14)

for some . Then, is is invertible and its inverse satisfies

 ∥B−1∥L(W,V)≤1α. (2.15)

Thus, for each , the problem (2.9) has a unique solution which satisfies the a priori estimate

 ∥uL∥V≤∥L∥Wα. (2.16)

For the particular and given by the left and right side of (2.3), the ellipticity condition (2.14) holds with under the assumption

 R(a(x))≥r,x∈D, (2.17)

since the latter implies, for all ,

 |B(v,v)|≥R(B(v,v))=∫DR(a)|∇v|2≥r∥∇v∥2L2=r∥v∥2V. (2.18)

Therefore, for any we may extend the solution map of the elliptic problem (1.5) to the complex domain

 Dr:={a∈X:R(a)≥r}, (2.19)

with the uniform bound

 ∥u∥L∞(Dr,V)=sup{∥u(a)∥V:a∈Dr}≤∥f∥Wr. (2.20)

This extension is therefore defined on the open set .

The fact that this map is holomorphic immediately follows by viewing it as a chain of holomorphic maps: introducing for any the operator acting from into , we can decompose into the chain of maps

 a↦B(a)↦B(a)−1↦B(a)−1f=u(a). (2.21)

The first and third maps are continuous linear and therefore holomorphic, from into and from onto respectively. The second map is the operator inversion which is holomorphic at any invertible .

For further purposes, it is interesting to compute the Frechet complex derivative for the elliptic problem (1.5). Fix an and let be such that . Then, the solution map is also defined at . Substracting the variational formulations (2.3) for and , we find that

 ∫Da∇(u(a+h)−u(a))⋅∇v=−∫Dh∇u(a+h)⋅∇v. (2.22)

We first use this identity to obtain a Lipschitz continuity bound: by taking and taking the real part of both sides, we find that

 r∥u(a+h)−u(a)∥2V≤∥h∥X∥u(a+h)∥V∥u(a+h)−u(a)∥V≤2∥f∥Wr∥h∥X∥u(a+h)−u(a)∥V,

and therefore

 ∥u(a+h)−u(a)∥V≤C∥h∥X,C:=2∥f∥Wr2. (2.23)

We next show that can be defined as the solution to the problem

 ∫Da∇w⋅∇v=−∫Dh∇u(a)⋅∇v,v∈V, (2.24)

which is well-posed in the sense of the above Lax-Milgram theorem. Indeed, on the one hand, the dependence of on is linear and it is continuous because taking we find that

 ∥w∥V≤C∥h∥X,C:=∥u(a)∥Vr. (2.25)

On the other hand, the remainder is the solution to

 ∫Da∇g⋅∇v=∫Dh(∇u(a)−∇u(a+h))⋅∇v,v∈V, (2.26)

which by taking and using (2.23) gives the quadratic bound

 ∥g∥V≤1r∥h∥X∥u(a)−u(a+h)∥V≤C∥h∥2X,C:=2∥f∥Wr3. (2.27)

This confirms that .

### 2.2 Extensions by the Ladyzhenskaya-Babushka-Brezzi theory

Our next goal is to treat more general linear parametric problems that are not necessarily elliptic. In particular, we have in mind parabolic problems such as the heat equation, or saddle points problems such as the Stokes equations. In order to formulate this general class of problems, we suppose that and are two complex Hilbert spaces with inner products and , respectively. So, for example, for every ,


text:
### Discussion 22 posts |  Only members can post. Join?

Title Author Replies Latest Post
Small wireless conitnuous light? Cayce Lee Photography 2 38 months ago
Which camera does long exposures without huge delay? glenneroo 3 38 months ago
Light Painting Experimental Doc earsaregood 0 39 months ago
New group - "What's happiness for you?" - Post your pics FLAVIA ADRIANE 5 42 months ago
Crap in Group pool! Admin do something Kushal Goyal 0 44 months ago
light a.m.i.snaps 0 47 months ago

LIGHT, where light is the primary subject. Since all photography involves light (there is one notable exception at the bottom of this page!) any photo or group might seem to fit. But this group is about photos where the main reason for taking the photo was to capture a special lighting effect.

For example, If I see a cool car I take a picture of it because I like that car. If I see a particular interplay of light and shadow on the hood then I take it for that reason, whether it's a particularly interesting car or not. Then if I get really lucky I can get the cool car and the interplay of light and shadow together in the same photo!

If your photo fits into any of the groups below (except DARK) then feel free to post it here!

Related groups in the [directory] system are =
BUG [directory] 虫
CAR [directory] 車
COLOR [directory] 色
CRITIQUE [directory] 批評
DIRECTORY [directory] ダイレクトリー
JAPAN [directory] 日本
PLANT [directory] 植物

!! . . . !! . . . !! . . . ATMOSPHERE . . . !! . . . !! . . . !!
Misty hazy smokey

!! . . . !! . . . !! . . . FIRE . . . !! . . . !! . . . !!
Ablazin'
Burning Sensations
CAMPERS & campfires
Campfires and bonfires
Fire
Fire and water
Fire Shots
Fire Spinning
Taken by Candlelight
Wildfire

!! . . . !! . . . !! . . . FIREWORKS . . . !! . . . !! . . . !!
Fireworks
Pyrotechnics
Sparklers

!! . . . !! . . . !! . . . LIGHT . . . !! . . . !! . . . !!
Beautiful Light
BLACKLIGHT Experience
bright light
Burnt ccd
Candles, Lamps & Moonlight
Chiaroscuro
EXPERIMENTS WITH LIGHT
the FANTASTIC
Illuminated Watercolors
The Light Fantastic
www.flickr.com/groups/light_camera_/
Light in the Face [FACES only]
lights & skin
Memories of shattered light
Playing with lights
Romantic Light
Urban Backlight / Contraluz Urbano

!! . . . !! . . . !! . . . LIGHT (artificial) . . . !! . . . !! . . . !!
Candles, Lamps & Moonlight
Christmas Lights
Colored Lights at Night
Fluorescent Lamps
Idea/One Lamp
Illumination
Lamps and Candles
Lava Lamp's We all love them.
Light Etching
Light Poles
Light pollution
Lighting Design
Lights and Lamps
Little Green Man - Pedestrian Traffic Lights Around the World
Little Lights
Pro Sound & Light
_public_digital_displays
Red Light/Green Light
Studio Lighting
Super Christmas House Lights
Theatrical Lighting Design

!! . . . !! . . . !! . . . LIGHT DRAWING & PAINTING . . . !! . . . !! . . . !!
painting with light
LIGHTDRAWING
LightPainting
Light Writings
smears of light
light junkies

!! . . . !! . . . !! . . . LIGHTNING . . . !! . . . !! . . . !!
Lightning
Lightning Images.

!! . . . !! . . . !! . . . LIGHT (low, available, long exposure) . . . !! . . . !! . . . !!
Available Light (a.k.a Existing Light)
2 seconds last long exposure
Candles, Lamps & Moonlight
Lightstream
LONGEXPOSURE
Long Exposure
Long Exposures
Long Exposure Times
Low light photography
No tripod - Now what?
A steady hand (1/4s+, no tripods allowed)

!! . . . !! . . . !! . . . LIGHT (night) . . . !! . . . !! . . . !!
Colored Lights at Night
Long exposures at night
Night Crawlers
Night Over-Exposure
night over water - night at the waterside -- nowat
Night Shot
Nightshots
No-Flash Night Shots
~~~The Warpedlight Zone~~~

!! . . . !! . . . !! . . . PLASMA . . . !! . . . !! . . . !!
plasma
plasmarites

!! . . . !! . . . !! . . . RAINBOW . . . !! . . . !! . . . !!
rainbow
Rainbow Warriors
Rainbows: created by light (not water)
Regenbogen Rainbow Colorfull

!! . . . !! . . . !! . . . SHADE & SHADOW . . . !! . . . !! . . . !!

!! . . . !! . . . !! . . . THE SKY . . . !! . . . !! . . . !!
sky & clouds
The Skies Of Texas
ThE sKy'S tHe LiMiT
storm light

!! . . . !! . . . !! . . . STARS, MOONS, PLANETS . . . !! . . . !! . . . !!
Astronomy
Astrophotography
Edmonton Observers Group Network
Star Trails

!! . . . !! . . . !! . . . THE SUN . . . !! . . . !! . . . !!
The Sun
solar eclipse
Solar Max

!! . . . !! . . . !! . . . SUNLIGHT . . . !! . . . !! . . . !!
The Blue Hour
Gilded
Heavenly Rays of Light
Lusco-Fusco
The Magic Hour
twilight and rust

!! . . . !! . . . !! . . . SUNRISE & SUNSET . . . !! . . . !! . . . !!
Fire in the Sky!
High Altitude Sunrises and Sunsets
Rojo atardecer - Red sunset
Scottish Sunsets and Sunrises
Sunrise and Sunset
Sunrise, Sunset -- Anything Sun!
Sun_Set
Sun sets
SUNSET
SUNSET&RISEPIX
Sunset Photo
Sunsets
Sunsets and sunrises around the world
Sunsets over the City
Top 20 Sunrise/Sunset Pix
Whitstable sunsets

!! . . . !! . . . !! . . . of Interest to Light Lovers . . . !! . . . !! . . . !!
images of ra

!! . . . !! . . . !! . . . The Opposite Group . . . !! . . . !! . . . !!
DARK　-　暗い

This is a public group.

• Members can post 50 things to the pool each month.
• Accepted media types:
• Photos
• Video
• Accepted content types:
• Photos / Videos
• Screenshots / Screencasts
• Illustration/Art / Animation/CGI
• Accepted safety levels:
• Safe
Feed – Subscribe to LIGHT [directory] 光 discussion threads


text:
# us76_u86_4-spectra#

The us76_u86_4-spectra data set is an absorption cross section data set for the the us76_u86_4 absorbing gas mixture and computed using the SPECTRA Information System.

## Description#

### Absorber#

The us76_u86_4 absorber is a gas mixture defined by the mixing ratios provided in the table below.

 Species N2 O2 CO2 CH4 Mixing ratio [%] 78.084 20.9476 0.0314 0.0002

Each gas species in this gas mixture is an absorber. The us76_u86_4 absorber is named in this way because it corresponds to the gas mixture defined by the us76 atmosphere model [NNU76] in the region of altitudes below 86 km, and restricted to the 4 main molecular species, namely N2, O2, CO2 and CH4.

### Data source#

The absorption cross sections data in this data set were computed using the SPECTRA Information System, with inputs available here. SPECTRA uses the HITRAN database (2016-edition) [GRH+17].

### Structure#

The us76_u86_4-spectra data set includes two data variables:

• absorption cross section (xs)

• mixing ratios (mr)

and three dimension coordinates:

• wavenumber (w)

• pressure (p) and

• molecule (m).

The absorption cross section data variable is tabulated with respect to the wavenumber and pressure coordinates. The mixing ration data variable is tabulated with respect to the pressure and molecule coordinates.

Note

The data set does not include a temperature dimension coordinate because the temperature coordinate is mapped onto the pressure coordinate according to the pressure-temperature relationship defined by the us76 thermophysical profile. This is why this data set is tied to the us76 thermophysical profile.

The graph below illustrates the pressure and temperature points used to compute the absorption cross section. There are 64 points (in red). The pressure values are geometrically spaced between $$10^{-6} \mathrm{atm}$$ ($$0.101325 \mathrm{Pa}$$) and $$1 \mathrm{atm}$$ ($$101325 \mathrm{Pa}$$). The temperature values are defined as the ordinate values when the abscissa values are set to these 64 pressure values, and the pressure-temperature relationship is given by the US76 model (blue curve).

#### Wavenumber coordinate#

The wavenumber coordinate is a linearly spaced mesh between $$4000 \, \mathrm{cm}^{-1}$$ ($$2500 \, \mathrm{nm}$$) and $$25711 \, \mathrm{cm}^{-1}$$ ($$390 \, \mathrm{nm}$$) with a constant wavenumber step of $$0.01 \, \mathrm{cm}^{-1}$$ (corresponding wavelength step varies between $$~ 0.00015 \, \mathrm{nm}$$ and $$~ 0.006 \, \mathrm{nm}$$).

#### Pressure coordinate#

The pressure coordinate is a geometrically spaced mesh between $$0.101325 \, \mathrm{Pa}$$ ($$10^{-6} \, \mathrm{atm}$$) and $$101325 \, \mathrm{Pa}$$ with 64 mesh points.

The minimum pressure value of $$10^{-6} \, \mathrm{atm}$$ corresponds to a maximum altitude of approximately 93 km, in the us76 thermophysical profile. This minimum value is a restriction of SPECTRA. For higher altitudes, that is to say for lower pressure values, the absorption cross section value is approximated to 0.

Note

Given that the maximum value of the absorption cross section at 93 km is:

$\max_{\nu} \sigma_{a} \approx 9 \, 10^{-23} \, \mathrm{cm}^2,$

which corresponds to a maximal absorption coefficient value of:

$\max_{\nu} k_{a} \approx 4 \, 10^{-4} \, \mathrm{km}^{-1},$

this approximation seems reasonable.

## Interpolation accuracy#

We assess the accuracy of the wavenumber and pressure interpolation of the us76_u86_4-spectra absorption data set.

### Method#

We down-sample the original data set on the corresponding axis, i.e. either wavenumber or pressure, by selecting data corresponding to even-indices coordinate values. Then, we interpolate the down-sampled data set on the odd-indices coordinate values. We assess the interpolation accuracy by computing the relative errors:

(1)#$\epsilon(p_i, q) = \frac{ \lvert \sigma_{\mathrm{a, interpolated}}(p_i, q) - \sigma_{\mathrm{a, original}}(p_i, q) \rvert}{ \sigma_{\mathrm{a, original}}(p_i, q)}$

where:

• $$\sigma_{\mathrm{a, original}}$$ is the original data set absorption cross section,

• $$\sigma_{\mathrm{a, interpolated}}$$ is the interpolated down-sampled data set absorption cross section,

• $$p_i$$ is the interpolation coordinate (i.e. wavenumber or pressure) where $$i$$ is odd, and

• $$q$$ is the other coordinate (i.e. pressure or wavenumber).

Note

Since the assessed interpolation accuracy here is that of the downsampled data set, we can expect that the interpolation accuracy of the original data set to be better.

We discard the cross section data where the following condition is met:

(2)#$k_{\mathrm{a}}(\nu_{i}, p) < 10^{-3} \, \mathrm{cm}^{-1},$

where $$k_{\mathrm{a}}(\nu_{i}, p)$$ is computed with:

$k_{\mathrm{a}}(\nu_{i}, p) = n(p) \, \sigma_{\mathrm{a}}(\nu_{i}, p)$

where $$n(p)$$ is the number density corresponding to the pressure $$p$$ in the us76 thermophysical profile.

We apply (2) because we estimate that these absorption coefficient values are too small to influence radiative transfer computations significantly.

Note

By interpolating the down-sampled absorption cross section data set at the odd-index coordinate values, we maximise the interpolation error with respect to the reference data set because each interpolation point is exactly at the middle of each data interval. This means that the relative errors are computed in the worst case scenarios, hence provide a conservative estimate of the interpolation accuracy.

### Results#

Below are illustrated some histograms of the values of $$\epsilon$$ computed with (1).

#### Wavenumber interpolation#

Below is given the histogram of the relative errors $$\epsilon (\nu_i, p)$$ for $$p=101325 \, \mathrm{Pa}$$ when interpolating the down-sampled data set on to the even-index wavenumber values of the original data set.

#### Pressure interpolation#

Below is illustrated the histogram of the relative errors $$\epsilon(p_i, \nu)$$ for $$p_i = 65349 \, \mathrm{Pa}$$ and for all wavenumbers values of the original data set.

### Analysis#

#### Wavenumber interpolation#

We make the following observations:

• At standard pressure (ground level altitude), the interpolation accuracy is fine – better than 5 % except for some outliers.

• As the pressure decreases (corresponding altitude increases), the interpolation accuracy gets poorer and poorer, although the errors counts decrease simultaneously.

• The counts number decreases with decreasing pressure because lower pressure means lower number density hence fewer absorption cross section data points satisfy (2).

#### Pressure interpolation#

We make the following observations:

• At ~4 km altitude ($$p = 65349 \, \mathrm{Pa}$$), the interpolation accuracy is relatively bad – up to 50 % – except for a few outliers.

• The interpolation accuracy generally improves with increasing altitudes (decreasing pressure).

• For altitudes ~ 21 km and higher ($$p < 4703 \, \mathrm{Pa}$$), the interpolation accuracy is better than 5 %, except for some outliers.

• The counts number decreases with decreasing pressure because lower pressure means lower number density hence fewer absorption cross section data points satisfy (2).

### Conclusion#

It is not easy to draw any conclusion regarding the interpolation accuracy of the us76_u86_4-spectra absorption data set from the above analysis. Since the interpolation accuracy of the down-sampled data set is not quite satisfactory, there is not guarantee that the interpolation accuracy of the actual data set is. On the other hand, the interpolation accuracy of the down-sampled data set was estimated in the worst case scenarios, so that the average accuracy is in fact better than the reported accuracy. Nevertheless, we plan to generate a new version of the data set with finer wavenumber and pressure meshes to improve the interpolation accuracy.


text:
# CommonMark testing started here!

We circled back on this, we may enable it by default in a few weeks, just living with it for a bit now

5 Likes

you would escape it like so: $\$a\$

I really do not think so, can you think of an example where this is needed?

No, i just seen something strange on math.stackexchange and asked question about reasons, because you have more experience in this area

Very marginally, I have only been doing MathJax for a few days

1 Like

Yeah, you have to use cookAsync from now on, will check in a deprecation message there.

https://github.com/discourse/discourse/blob/master/app/assets/javascripts/discourse/lib/text.js.es6#L26-L33

So your new code would be something like:

  let props = {
raw: new_raw,
edit_reason: 'checklist change',
};

cookAsync(new_raw).then(cooked => {
props.cooked = cooked.string;
viewPost.save(props);
});

3 Likes

Will image alignment be made possible with commonmark?

Not much progress in the CommonMark world there…

But … since we already only use HTML based image metadata, whitelisting a few classes is super easy in a plugin.

@codinghorror I wonder if we should just whitelist a few here out of the box?

2 Likes
4 Likes

Cool - thanks for thinking about image alignment, which may seem unimportant in the grand scheme of things but would be a big help to us.

Right now, my colleague is trying to display some profile pictures with biographical details attractively over here on my discourse.

Interestingly, I was able to get something close to attractive on the commonmark.js demo site:

Here’s the code I used to create that:

[<img src="https://community.namati.org/uploads/default/original/2X/4/44c8fd629b36d97fe0c6e618fc99bb84497768bd.jpg" width="170" height="125" align="left">](https://namati.org/network/member/LucyClaridge) Lucy Claridge (@LucyClaridge)  is Legal Director at Minority Rights Group International, where she manages the legal department.  Her work uses strategic litigation, advocacy and capacity building to improve access to justice for minority and indigenous communities worldwide, with a particular focus on land rights, political participation and wider anti-discrimination issues.  Casework includes the recent Ogiek land and conservation case against Kenya before the African Court of Human and Peoples' Rights.

Nah we are fine as is. I don’t want to whitelist a ton of random stuff that we don’t need, every one is a security risk in the long term.

@sam, so that works great for when markdown-it is enabled, however, when it is disabled, I now get

ember:19818 Error: Could not find module pretty-text/engines/markdown-it/helpers imported from discourse/plugins/discourse-plugin-checklist/lib/discourse-markdown/checklist
at requireFrom (discourse-loader:128:13)
at reify (discourse-loader:106:22)
at mod.state (discourse-loader:163:17)
at tryFinally (discourse-loader:30:14)
at require (discourse-loader:162:5)
at eval (pretty-text/engines/discourse-markdown:381:26)
at Array.forEach ()
at DialectHelper.setup (pretty-text/engines/discourse-markdown:379:40)
at setup (pretty-text/engines/discourse-markdown:456:12)
at buildOptions (pretty-text/pretty-text:59:36)

It seems to not like my importing of the helper when markdown-it is not enabled. Is there a way to avoid that error? I tried using

  if (helper.markdownIt)
{
const prettyTextHelper = require('pretty-text/engines/markdown-it/helpers');
helper.registerPlugin(md => {
const ruler = md.inline.ruler;
ruler.push('checklist-strikethrough', prettyTextHelper.inlineRegexRule(md, {
start: '--',
matcher: /--(.*)--/i,
emitter: applyStrikethrough
}));

ruler.push('checklist-empty-checkbox', prettyTextHelper.inlineRegexRule(md, {
start: '[',
matcher: /$([\s_\-x\*]{0,1})$/i,
emitter: applyCheckbox
}));
});
}

versus an import at the top, but that didn’t work. Otherwise, I may just have to leave this in a branch until markdown-it becomes the default (which I’m not opposed to)

Yeah… you are going to need this PR

https://github.com/cpradio/discourse-plugin-checklist/pull/9

Note you can do the strikethrough like you had with the inline regex, but checkboxes you can not. You need to apply at the end of the chain otherwise there be

Just import the helper using:

// import like this for backwards compat, add a comment
require('pretty-text/engines/markdown-it/helpers').inlineRegexRule

I think I will clean the API up a bit so we just hang the helpers of the md object then you don’t need to import anything.

Hmm, I’m struggling getting that PR to work. I’ve pulled it in locally, kicked up the docker dev environment, but I can’t get the checkboxes to show with or without markdown-it enabled.

I’ve even completely destroyed the docker dev environment and re-initialized.

Edit:
Figured it out, this errant return; statement
https://github.com/cpradio/discourse-plugin-checklist/pull/9/commits/1ff70f501da887b1486a2d97f265e2d6a7dbf298#diff-58d619630012ba3d44eea8bdb4f2440bR114

However, there is already : with the current implementation, given the following, it breaks

[] Item One
[*] Item Two (should be checked)
[ ] Item Three
[*] Item Four (should also be checked)

Actual output:, Items Two through Four are italized.

I’ve got that fixed, however, now with markdown-it turned off, the fix for the italize is causing checking a checkbox to not show properly (preview pane shows it properly though when you edit the post afterwards).

It is only breaking with markdown-it disabled. But for the life of me, I haven’t a clue why, as when I trace it down, it is replacing [\*] with the appropriate span tag and classes, but when it saves the post, the [\*] reappears from seemingly nowhere…

Granted, this can happen with [_] as well on new markdown-it, but that value can’t be generated by clicking the checkbox, so the user should see the issue and escape it before it is submitted. Only [*] is capable of being written based on the user clicking a box.

Updated PR to look at (which has your PR and the fixes I described above – and the issue with markdown-it disabled)
https://github.com/cpradio/discourse-plugin-checklist/pull/9

2 Likes

I see that headings have changed from ##heading to ## heading (with a space).
What happens to the old topics that use ##heading instead of ## heading after the upgrade?

We have to do a remap as describe here Replace a string in all posts ?

If you don’t rebake old posts they will keep the current generated HTML.

3 Likes

And become a time bomb, waiting to break in the future.

7 Likes

I’ve already seen a topic of mine here on Meta break. Sub-bullets now need two spaces, not one, so when someone edited a wiki post of mine my bulleted list broke.

I wonder @sam if any sort of migration script is planned?

1 Like

Very super undecided on a migration script. The concern I have is that it is possible a migration script can mess stuff up more than it will actually help.

I am semi open to making a script that fixes up bad multiline quotes (one’s that start mid line) and maybe the title thing.

Will think about it.

4 Likes

If you store bake time and switch to new engine time, you can punt on a script to deal with later, or to just highlight warnings when appropriate.

“You are editing a post that predates the current parser. Please double check that the syntax still works as expected.”

2 Likes

No, this is not needed at all. I would not bother. People can fix up the rare affected posts themselves.

2 Likes


text:
Cancel
Showing results for
Did you mean:

# Femap Model Bounding Box

Community Manager

# Femap Tips and Tricks: Model Box

What is the Femap Model Bounding Box, and why is it important?

The Femap Model Bounding Box is a hexahedral cuboid constructed from planes that bound the model
in the x, y, and z axes. The longest diagonal of the model bounding box is used to determine certain
default values in Femap.

For example, the default tolerance used when merging entities is calculated to be bounding box
diagonal divided by by 10,000.

Merge tolerance = Diagonal / 10,000

Also the default maximum deformation is calculated to be the bounding box diagonal multiplied by the
applied scale factor and root 3, to average the edge length of the model bounding box.

Maximum deformation = Scale factor * Diagonal * √3

(view in My Videos)

• ,
• ### Femap Tip and Trick

Contributors
Latest articles
Subject Date Posted

by

by

by