source: ThirdParty/vmg/src/solver/solver_regular.cpp@ 33a694

Action_Thermostats Add_AtomRandomPerturbation Add_RotateAroundBondAction Add_SelectAtomByNameAction Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_StructOpt_integration_tests AutomationFragmentation_failures Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator Docu_Python_wait EmpiricalPotential_contain_HomologyGraph_documentation Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph Fix_ChronosMutex Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion GeometryObjects Gui_displays_atomic_force_velocity IndependentFragmentGrids_IntegrationTest JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks RotateToPrincipalAxisSystem_UndoRedo StoppableMakroAction TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps Ubuntu_1604_changes
Last change on this file since 33a694 was 7faa5c, checked in by Frederik Heber <heber@…>, 8 years ago

Merge commit 'de061d9d851257a04e924d4472df4523d33bb08b' as 'ThirdParty/vmg'

  • Property mode set to 100644
File size: 3.9 KB
Line 
1/*
2 * vmg - a versatile multigrid solver
3 * Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
4 *
5 * vmg is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * vmg is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19/**
20 * @file solver_regular.cpp
21 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
22 * @date Mon Apr 18 13:11:32 2011
23 *
24 * @brief VMG::SolverRegular
25 *
26 */
27
28#ifdef HAVE_CONFIG_H
29#include <libvmg_config.h>
30#endif
31
32#include <iostream>
33
34#include "base/discretization.hpp"
35#include "base/stencil.hpp"
36#include "grid/multigrid.hpp"
37#include "solver/solver_regular.hpp"
38#include "mg.hpp"
39
40using namespace VMG;
41
42// TODO: Implement global communication here
43// TODO: Implement this more efficiently
44// TODO: This has to be reviewed for parallelization
45
46void SolverRegular::AssembleMatrix(const Grid& rhs)
47{
48 Grid::iterator grid_iter;
49 Stencil::iterator stencil_iter;
50 int mat_index, mat_index2;
51 vmg_float prefactor_inv = 1.0 / MG::GetDiscretization()->OperatorPrefactor(rhs);
52 const Stencil& A = MG::GetDiscretization()->GetStencil();
53
54#ifdef DEBUG_MATRIX_CHECKS
55 rhs.IsConsistent();
56#endif
57
58 this->Realloc(rhs.Global().GlobalSize().Product());
59
60 for (grid_iter = rhs.Iterators().Local().Begin(); grid_iter != rhs.Iterators().Local().End(); ++grid_iter) {
61
62 mat_index = rhs.GlobalLinearIndex(*grid_iter + rhs.Global().LocalBegin());
63
64 assert(mat_index >= 0 && mat_index<this->Size());
65
66 this->Sol(mat_index) = 0.0;
67 this->Rhs(mat_index) = prefactor_inv * rhs.GetVal(*grid_iter);
68
69 for (int l=0; l<this->Size(); l++)
70 this->Mat(mat_index, l) = 0.0;
71
72 this->Mat(mat_index, mat_index) = A.GetDiag();
73
74 for (stencil_iter = A.begin(); stencil_iter != A.end(); ++stencil_iter) {
75
76 mat_index2 = rhs.GlobalLinearIndex(*grid_iter + rhs.Global().LocalBegin() + stencil_iter->Disp());
77
78 assert(mat_index2 >= 0 && mat_index2<this->Size());
79
80 this->Mat(mat_index, mat_index2) += stencil_iter->Val();
81
82 }
83 }
84
85 for (int i=0; i<3; ++i) {
86
87 for (grid_iter = rhs.Iterators().Boundary1()[i].Begin(); grid_iter != rhs.Iterators().Boundary1()[i].End(); ++grid_iter) {
88
89 mat_index = rhs.GlobalLinearIndex(*grid_iter + rhs.Global().LocalBegin());
90
91 assert(mat_index >= 0 && mat_index<this->Size());
92
93 this->Sol(mat_index) = this->Rhs(mat_index) = rhs.GetVal(*grid_iter);
94
95 for (int l=0; l<this->Size(); l++)
96 this->Mat(mat_index, l) = 0.0;
97
98 this->Mat(mat_index, mat_index) = 1.0;
99
100 }
101
102 for (grid_iter = rhs.Iterators().Boundary2()[i].Begin(); grid_iter != rhs.Iterators().Boundary2()[i].End(); ++grid_iter) {
103
104 mat_index = rhs.GlobalLinearIndex(*grid_iter + rhs.Global().LocalBegin());
105
106 assert(mat_index >= 0 && mat_index<this->Size());
107
108 this->Sol(mat_index) = this->Rhs(mat_index) = rhs.GetVal(*grid_iter);
109
110 for (int l=0; l<this->Size(); l++)
111 this->Mat(mat_index, l) = 0.0;
112
113 this->Mat(mat_index, mat_index) = 1.0;
114
115 }
116
117 }
118
119}
120
121void SolverRegular::ExportSol(Grid& sol, Grid& rhs)
122{
123 int index;
124 Index offset;
125
126 for (int i=0; i<3; ++i)
127 offset[i] = (sol.Local().HaloEnd1()[i] > 0 ? 1 : 0);
128
129 for (Grid::iterator iter = sol.Iterators().CompleteGrid().Begin(); iter != sol.Iterators().CompleteGrid().End(); ++iter) {
130 index = sol.GlobalLinearIndex(sol.Global().LocalBegin() + *iter - offset);
131 sol(*iter) = this->Sol(index);
132 }
133
134#ifdef DEBUG_MATRIX_CHECKS
135 sol.IsConsistent();
136#endif
137}
Note: See TracBrowser for help on using the repository browser.