GlobalLinSysIterativeFull.cpp 10.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
///////////////////////////////////////////////////////////////////////////////
//
// File GlobalLinSys.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: GlobalLinSys definition
//
///////////////////////////////////////////////////////////////////////////////

36
37
#include <map>
#include <MultiRegions/GlobalLinSysIterativeFull.h>
38
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
39
40
41
42
43
44
45
46

namespace Nektar
{
    namespace MultiRegions
    {
        /**
         * @class GlobalLinSysIterativeCG
         *
47
48
49
         * This class implements a conjugate gradient matrix solver.
         * Preconditioning is implemented using a Jacobi (diagonal)
         * preconditioner.
50
51
52
53
54
         */

        /**
         * Registers the class with the Factory.
         */
55
        string GlobalLinSysIterativeFull::className
56
                = GetGlobalLinSysFactory().RegisterCreatorFunction(
57
58
59
                    "IterativeFull",
                    GlobalLinSysIterativeFull::create,
                    "Iterative solver for full matrix system.");
60
61


62
63
64
65
66
67
68
        /**
         * Constructor for full direct matrix solve.
         * @param   pKey        Key specifying matrix to solve.
         * @param   pExp        Shared pointer to expansion list for applying
         *                      matrix evaluations.
         * @param   pLocToGloMap Local to global mapping.
         */
69
        GlobalLinSysIterativeFull::GlobalLinSysIterativeFull(
70
                    const GlobalLinSysKey &pKey,
71
                    const boost::weak_ptr<ExpList> &pExp,
72
                    const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
73
                : GlobalLinSysIterative(pKey, pExp, pLocToGloMap)
74
        {
75
            ASSERTL1(m_linSysKey.GetGlobalSysSolnType()==eIterativeFull,
76
77
                     "This routine should only be used when using an Iterative "
                     "conjugate gradient matrix solve.");
78
79
        }

80
81
82
83

        /**
         *
         */
84
        GlobalLinSysIterativeFull::~GlobalLinSysIterativeFull()
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
        {

        }


        /**
         * Solve a global linear system with Dirichlet forcing using a
         * conjugate gradient method. This routine performs handling of the
         * Dirichlet forcing terms and wraps the underlying iterative solver
         * used for the remaining degrees of freedom.
         *
         * Consider solving for \f$x\f$, the matrix system \f$Ax=b\f$, where
         * \f$b\f$ is known. To enforce the Dirichlet terms we instead solve
         * \f[A(x-x_0) = b - Ax_0 \f]
         * where \f$x_0\f$ is the Dirichlet forcing.
         *
         * @param           pInput      RHS of linear system, \f$b\f$.
         * @param           pOutput     On input, values of dirichlet degrees
103
104
         *                              of freedom with initial guess on other values.
         *                              On output, the solution \f$x\f$.
105
106
107
         * @param           pLocToGloMap    Local to global mapping.
         * @param           pDirForcing Precalculated Dirichlet forcing.
         */
108
        void GlobalLinSysIterativeFull::v_Solve(
109
110
                    const Array<OneD, const NekDouble>  &pInput,
                          Array<OneD,       NekDouble>  &pOutput,
111
                    const AssemblyMapSharedPtr &pLocToGloMap,
112
                    const Array<OneD, const NekDouble>  &pDirForcing)
113
        {
114
            boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
115
            bool vCG;
Dave Moxey's avatar
Dave Moxey committed
116
117
            if ((m_locToGloMap = boost::dynamic_pointer_cast<AssemblyMapCG>(
                     pLocToGloMap)))
118
119
120
            {
                vCG = true;
            }
Dave Moxey's avatar
Dave Moxey committed
121
122
            else if ((m_locToGloMap = boost::dynamic_pointer_cast<
                          AssemblyMapDG>(pLocToGloMap)))
123
124
125
            {
                vCG = false;
            }
126
127

            bool dirForcCalculated = (bool) pDirForcing.num_elements();
128
129
130
            int nDirDofs  = pLocToGloMap->GetNumGlobalDirBndCoeffs();
            int nGlobDofs = pLocToGloMap->GetNumGlobalCoeffs();
            int nDirTotal = nDirDofs;
131
            
132
            expList->GetComm()->AllReduce(nDirTotal, LibUtilities::ReduceSum);
133
            
Dave Moxey's avatar
Dave Moxey committed
134
            Array<OneD, NekDouble> tmp(nGlobDofs), tmp2;
135
136

            if(nDirTotal)
137
138
139
140
141
142
            {
                // calculate the Dirichlet forcing
                if(dirForcCalculated)
                {
                    Vmath::Vsub(nGlobDofs, pInput.get(), 1,
                                pDirForcing.get(), 1,
143
                                tmp.get(), 1);
144
145
146
147
148
                }
                else
                {
                    // Calculate the dirichlet forcing B_b (== X_b) and
                    // substract it from the rhs
149
                    expList->GeneralMatrixOp(
150
                        m_linSysKey, pOutput, tmp, eGlobal);
151

152
153
154
                    Vmath::Vsub(nGlobDofs, pInput.get(), 1,
                                           tmp.get(),    1,
                                           tmp.get(),    1);
155
                }
156
157
                if (vCG)
                {
158
                    Array<OneD, NekDouble> out(nGlobDofs,0.0);
Dave Moxey's avatar
Dave Moxey committed
159

160
                    // solve for perturbation from intiial guess in pOutput
Dave Moxey's avatar
Dave Moxey committed
161
162
163
164
                    SolveLinearSystem(
                        nGlobDofs, tmp, out, pLocToGloMap, nDirDofs);
                    Vmath::Vadd(nGlobDofs-nDirDofs,    &out    [nDirDofs], 1,
                                &pOutput[nDirDofs], 1, &pOutput[nDirDofs], 1);
165
166
167
168
169
                }
                else
                {
                    ASSERTL0(false, "Need DG solve if using Dir BCs");
                }
170
171
172
            }
            else
            {
Chris Cantwell's avatar
Chris Cantwell committed
173
                Vmath::Vcopy(nGlobDofs, pInput, 1, tmp, 1);
174
                SolveLinearSystem(nGlobDofs, tmp, pOutput, pLocToGloMap);
175
176
177
            }
        }

178

179
180
181
182
183
184
185
        /**
         *
         */
        void GlobalLinSysIterativeFull::v_DoMatrixMultiply(
                const Array<OneD, NekDouble>& pInput,
                      Array<OneD, NekDouble>& pOutput)
        {
186
            boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
187
            // Perform matrix-vector operation A*d_i
188
            expList->GeneralMatrixOp(m_linSysKey,
189
                                     pInput, pOutput, eGlobal);
190
191
192
193

            // retrieve robin boundary condition information and apply robin
            // boundary conditions to the solution.
            const std::map<int, RobinBCInfoSharedPtr> vRobinBCInfo
194
                                                = expList->GetRobinBCInfo();
195
196
197
198
199
200
201
202
203
204
            if(vRobinBCInfo.size() > 0)
            {
                ASSERTL0(false,
                        "Robin boundaries not set up in IterativeFull solver.");
                int nGlobal = m_locToGloMap->GetNumGlobalCoeffs();
                int nLocal  = m_locToGloMap->GetNumLocalCoeffs();
                int nDir    = m_locToGloMap->GetNumGlobalDirBndCoeffs();
                int nNonDir = nGlobal - nDir;
                Array<OneD, NekDouble> robin_A(nGlobal, 0.0);
                Array<OneD, NekDouble> robin_l(nLocal,  0.0);
205
206
207
                Array<OneD, NekDouble> tmp;
                NekVector<NekDouble> robin(nNonDir,
                                           tmp = robin_A + nDir, eWrapper);
208
209
210
211
212
213
214

                // Operation: p_A = A * d_A
                // First map d_A to local solution
                m_locToGloMap->GlobalToLocal(pInput, robin_l);

                // Iterate over all the elements computing Robin BCs where
                // necessary
215
                for (int n = 0; n < expList->GetNumElmts(); ++n)
216
                {
217
218
219
                    int nel = expList->GetOffset_Elmt_Id(n);
                    int offset = expList->GetCoeff_Offset(n);
                    int ncoeffs = expList->GetExp(nel)->GetNcoeffs();
220
221
222
223
224

                    if(vRobinBCInfo.count(nel) != 0) // add robin mass matrix
                    {
                        RobinBCInfoSharedPtr rBC;
                        Array<OneD, NekDouble> tmp;
225
                        StdRegions::StdExpansionSharedPtr vExp = expList->GetExp(nel);
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246

                        // add local matrix contribution
                        for(rBC = vRobinBCInfo.find(nel)->second;rBC; rBC = rBC->next)
                        {
                            vExp->AddRobinEdgeContribution(rBC->m_robinID,rBC->m_robinPrimitiveCoeffs, tmp = robin_l + offset);
                        }
                    }
                    else
                    {
                        Vmath::Zero(ncoeffs, &robin_l[offset], 1);
                    }
                }

                // Map local Robin contribution back to global coefficients
                m_locToGloMap->LocalToGlobal(robin_l, robin_A);
                // Add them to the output of the GeneralMatrixOp
                Vmath::Vadd(nGlobal, pOutput, 1, robin_A, 1, pOutput, 1);
            }

        }

247
248
249
250
251
252
253
254
        /**
         *
         */
        void GlobalLinSysIterativeFull::v_UniqueMap()
        {
            m_map = m_locToGloMap->GetGlobalToUniversalMapUnique();
        }

255
256
    }
}