source: ThirdParty/vmg/src/smoother/gsrb_poisson_4.cpp@ be848d

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

Merge commit 'de061d9d851257a04e924d4472df4523d33bb08b' as 'ThirdParty/vmg'
  • Property mode set to 100644
File size: 5.5 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 gsrb_poisson_4.cpp
21 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
22 * @date Fri May 11 18:30:20 2012
23 *
24 * @brief Gauss-Seidel Red Black method, specialized to
25 * the Poisson equation. Performance improved by
26 * explicit loop unrolling.
27 *
28 */
29
30#ifdef HAVE_CONFIG_H
31#include <libvmg_config.h>
32#endif
33
34#ifdef HAVE_MPI
35#include <mpi.h>
36#endif
37
38#include "base/helper.hpp"
39#include "comm/comm.hpp"
40#include "grid/grid.hpp"
41#include "smoother/gsrb_poisson_4.hpp"
42#include "mg.hpp"
43
44using namespace VMG;
45
46static inline void ComputePartial(Grid& sol, Grid& rhs,
47 const Index& begin, const Index& end,
48 const vmg_float& prefactor, const int& off)
49{
50 const vmg_float fac_1 = 1.0 / 12.0;
51 const vmg_float fac_2 = 1.0 / 24.0;
52
53 for (int i=begin.X(); i<end.X(); ++i)
54 for (int j=begin.Y(); j<end.Y(); ++j)
55 for (int k=begin.Z() + (i + j + begin.Z() + off) % 2; k<end.Z(); k+=2)
56 sol(i,j,k) = prefactor * rhs.GetVal(i,j,k) + fac_1 * (sol.GetVal(i-1,j ,k ) +
57 sol.GetVal(i+1,j ,k ) +
58 sol.GetVal(i ,j-1,k ) +
59 sol.GetVal(i ,j+1,k ) +
60 sol.GetVal(i ,j ,k-1) +
61 sol.GetVal(i ,j ,k+1))
62 + fac_2 * (sol.GetVal(i-1,j-1,k ) +
63 sol.GetVal(i-1,j+1,k ) +
64 sol.GetVal(i+1,j-1,k ) +
65 sol.GetVal(i+1,j+1,k ) +
66 sol.GetVal(i-1,j ,k-1) +
67 sol.GetVal(i-1,j ,k+1) +
68 sol.GetVal(i+1,j ,k-1) +
69 sol.GetVal(i+1,j ,k+1) +
70 sol.GetVal(i ,j-1,k-1) +
71 sol.GetVal(i ,j-1,k+1) +
72 sol.GetVal(i ,j+1,k-1) +
73 sol.GetVal(i ,j+1,k+1));
74}
75
76void GaussSeidelRBPoisson4::Compute(Grid& sol, Grid& rhs)
77{
78 const vmg_float prefactor_inv = Helper::pow_2(sol.Extent().MeshWidth().Max()) / 4.0;
79 const int off = rhs.Global().LocalBegin().Sum() - rhs.Global().GlobalBegin().Sum() + rhs.Local().HaloSize1().Sum();
80 const LocalIndices& local = rhs.Local();
81 Comm& comm = *MG::GetComm();
82
83 /*
84 * Compute first halfstep
85 */
86
87 // Start asynchronous communication
88 comm.CommToGhostsAsyncStart(sol);
89
90 // Smooth part not depending on ghost cells
91 ComputePartial(sol, rhs,
92 local.Begin()+1, local.End()-1,
93 prefactor_inv, off+1);
94
95 // Finish asynchronous communication
96 comm.CommToGhostsAsyncFinish(sol);
97
98 /*
99 * Smooth near boundary cells
100 */
101
102 ComputePartial(sol, rhs,
103 local.Begin(),
104 Index(local.Begin().X()+1, local.End().Y(), local.End().Z()),
105 prefactor_inv, off+1);
106
107 ComputePartial(sol, rhs,
108 Index(local.End().X()-1, local.Begin().Y(), local.Begin().Z()),
109 local.End(),
110 prefactor_inv, off+1);
111
112 ComputePartial(sol, rhs,
113 Index(local.Begin().X()+1, local.Begin().Y(), local.Begin().Z()),
114 Index(local.End().X()-1, local.Begin().Y()+1, local.End().Z()),
115 prefactor_inv, off+1);
116
117 ComputePartial(sol, rhs,
118 Index(local.Begin().X()+1, local.End().Y()-1, local.Begin().Z()),
119 Index(local.End().X()-1, local.End().Y(), local.End().Z()),
120 prefactor_inv, off+1);
121
122 ComputePartial(sol, rhs,
123 Index(local.Begin().X()+1, local.Begin().Y()+1, local.Begin().Z()),
124 Index(local.End().X()-1, local.End().Y()-1, local.Begin().Z()+1),
125 prefactor_inv, off+1);
126
127 ComputePartial(sol, rhs,
128 Index(local.Begin().X()+1, local.Begin().Y()+1, local.End().Z()-1),
129 Index(local.End().X()-1, local.End().Y()-1, local.End().Z()),
130 prefactor_inv, off+1);
131
132 /*
133 * Compute second halfstep
134 */
135
136 // Start asynchronous communication
137 comm.CommToGhostsAsyncStart(sol);
138
139 // Smooth part not depending on ghost cells
140 ComputePartial(sol, rhs,
141 local.Begin()+1, local.End()-1,
142 prefactor_inv, off);
143
144 // Finish asynchronous communication
145 comm.CommToGhostsAsyncFinish(sol);
146
147 /*
148 * Smooth near boundary cells
149 */
150
151 ComputePartial(sol, rhs,
152 local.Begin(),
153 Index(local.Begin().X()+1, local.End().Y(), local.End().Z()),
154 prefactor_inv, off);
155
156 ComputePartial(sol, rhs,
157 Index(local.End().X()-1, local.Begin().Y(), local.Begin().Z()),
158 local.End(),
159 prefactor_inv, off);
160
161 ComputePartial(sol, rhs,
162 Index(local.Begin().X()+1, local.Begin().Y(), local.Begin().Z()),
163 Index(local.End().X()-1, local.Begin().Y()+1, local.End().Z()),
164 prefactor_inv, off);
165
166 ComputePartial(sol, rhs,
167 Index(local.Begin().X()+1, local.End().Y()-1, local.Begin().Z()),
168 Index(local.End().X()-1, local.End().Y(), local.End().Z()),
169 prefactor_inv, off);
170
171 ComputePartial(sol, rhs,
172 Index(local.Begin().X()+1, local.Begin().Y()+1, local.Begin().Z()),
173 Index(local.End().X()-1, local.End().Y()-1, local.Begin().Z()+1),
174 prefactor_inv, off);
175
176 ComputePartial(sol, rhs,
177 Index(local.Begin().X()+1, local.Begin().Y()+1, local.End().Z()-1),
178 Index(local.End().X()-1, local.End().Y()-1, local.End().Z()),
179 prefactor_inv, off);
180}
Note: See TracBrowser for help on using the repository browser.