PLaSK library

Loading...
Searching...
No Matches
bisection.cpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program 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, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#include "bisection.hpp"
15#include <plask/parallel.hpp>
16
17 namespace plask { namespace optical { namespace effective {
18
19#define CALL_FUN(result, re, im) \
20 if (error) continue; \
21 try { \
22 result = fun(dcomplex(re, im)); \
23 } catch (...) { \
24 PLASK_PRAGMA(omp critical) \
25 error = std::current_exception();\
26 }
27
28 Contour::Contour(const Solver* solver, const std::function<dcomplex(dcomplex)>& fun, dcomplex corner0, dcomplex corner1, size_t ren, size_t imn):
29 solver(solver), fun(fun), re0(real(corner0)), im0(imag(corner0)), re1(real(corner1)), im1(imag(corner1))
30{
31 bottom.reset(ren+1);
32 right.reset(imn+1);
33 top.reset(ren+1);
34 left.reset(imn+1);
35
36 double dr = (re1 - re0) / double(ren);
37 double di = (im1 - im0) / double(imn);
38
39 std::exception_ptr error;
40 PLASK_OMP_PARALLEL
41 {
42 #pragma omp for nowait
43 for (int i = 0; i < int(ren); ++i) {
44 CALL_FUN(bottom[i], re0+i*dr, im0)
45 }
46 #pragma omp for nowait
47 for (int i = 0; i < int(imn); ++i) {
48 CALL_FUN(right[i], re1, im0+i*di)
49 }
50 #pragma omp for nowait
51 for (int i = 1; i <= int(ren); ++i) {
52 CALL_FUN(top[i], re0+i*dr, im1)
53 }
54 #pragma omp for
55 for (int i = 1; i <= int(imn); ++i) {
56 CALL_FUN(left[i], re0, im0+i*di)
57 }
58 }
59 if (error) std::rethrow_exception(error);
60 // Wrap values
61 bottom[ren] = right[0];
62 right[imn] = top[ren];
63 top[0] = left[imn];
64 left[0] = bottom[0];
65}
66
67 namespace detail {
68 // Inform that there is zero very close to the boundary
69 static void contourLogZero(size_t i, size_t n, const Solver* solver, double r0, double i0, double r1, double i1) {
70 const double f = double(2*std::ptrdiff_t(i)-1) / double(2*std::ptrdiff_t(n)-2);
71 const double re = r0 + f*(r1-r0), im = i0 + f*(i1-i0);
72 solver->writelog(LOG_WARNING, "Zero at contour in {0} (possibly not counted)", str(dcomplex(re,im)));
73 }
74}
75
76 int Contour::crossings(const DataVector<dcomplex>& line, double r0, double i0, double r1, double i1) const
77{
78 int wind = 0;
79 for (size_t i = 1; i < line.size(); ++i) {
80 if (real(line[i-1]) < 0. && real(line[i]) < 0.) {
81 if (imag(line[i-1]) >= 0. && imag(line[i]) < 0.) {
82 if (real(line[i-1]) >= 0. || real(line[i]) >= 0.) detail::contourLogZero(i, line.size(), solver, r0, i0, r1, i1);
83 ++wind;
84 } else if (imag(line[i-1]) < 0. && imag(line[i]) >= 0.) {
85 if (real(line[i-1]) >= 0. || real(line[i]) >= 0.) detail::contourLogZero(i, line.size(), solver, r0, i0, r1, i1);
86 --wind;
87 }
88 }
89 }
90 return wind;
91}
92
93#ifdef NDEBUG
94std::pair<Contour,Contour> Contour::divide(double /*reps*/, double ieps) const
95#else
96std::pair<Contour,Contour> Contour::divide(double reps, double ieps) const
97#endif
98{
99 Contour contoura(solver, fun), contourb(solver, fun);
100
101 assert(re1-re0 > reps || im1-im0 > ieps);
102
103 if (bottom.size() > right.size() || im1-im0 <= ieps) { // divide real axis
104 double re = 0.5 * (re0 + re1);
105 contoura.re0 = re0; contoura.im0 = im0; contoura.re1 = re; contoura.im1 = im1;
106 contourb.re0 = re; contourb.im0 = im0; contourb.re1 = re1; contourb.im1 = im1;
107
108 size_t n = (bottom.size()-1) / 2; // because of the > in the condition, bottom.size() is always >= 3
109
110 size_t imn = right.size() - 1;
111 DataVector<dcomplex> middle(imn+1);
112 middle[0] = bottom[n]; middle[imn] = top[n];
113 double di = (im1 - im0) / double(imn);
114
115 std::exception_ptr error;
116 PLASK_OMP_PARALLEL_FOR
117 for (plask::openmp_size_t i = 1; i < imn; ++i) {
118 if (error) continue;
119 try {
120 middle[i] = fun(dcomplex(re, im0+double(i)*di));
121 } catch (...) {
122 error = std::current_exception();
123 }
124 }
125 if (error) std::rethrow_exception(error);
126
127 contoura.left = left;
128 contoura.right = middle;
129 contoura.bottom = DataVector<dcomplex>(const_cast<dcomplex*>(&bottom[0]), n+1);
130 contoura.top = DataVector<dcomplex>(const_cast<dcomplex*>(&top[0]), n+1);
131
132 contourb.left = middle;
133 contourb.right = right;
134 contourb.bottom = DataVector<dcomplex>(const_cast<dcomplex*>(&bottom[n]), n+1);
135 contourb.top = DataVector<dcomplex>(const_cast<dcomplex*>(&top[n]), n+1);
136
137 return std::make_pair(std::move(contoura), std::move(contourb));
138
139 } else { // divide imaginary axis
140 double im = 0.5 * (im0 + im1);
141 contoura.re0 = re0; contoura.im0 = im0; contoura.re1 = re1; contoura.im1 = im;
142 contourb.re0 = re0; contourb.im0 = im; contourb.re1 = re1; contourb.im1 = im1;
143
144 size_t ren = bottom.size() - 1;
145 DataVector<dcomplex> middle(ren+1);
146 if (right.size() <= 2) {
147 assert(ren == 1); // real axis also has only one segment
148 middle[0] = fun(dcomplex(re0, im));
149 middle[1] = fun(dcomplex(re1, im));
150 contoura.left = DataVector<dcomplex>({left[0], middle[0]});
151 contoura.right = DataVector<dcomplex>({right[0], middle[1]});
152 contourb.left = DataVector<dcomplex>({middle[0], left[1]});
153 contourb.right = DataVector<dcomplex>({middle[1], right[1]});
154 } else {
155 size_t n = (right.size()-1) / 2; // no less than 1
156 middle[0] = left[n]; middle[ren] = right[n];
157 double dr = (re1 - re0) / double(ren);
158
159 std::exception_ptr error;
160 PLASK_OMP_PARALLEL_FOR
161 for (plask::openmp_size_t i = 1; i < ren; ++i) {
162 if (error) continue;
163 try {
164 middle[i] = fun(dcomplex(re0+double(i)*dr, im));
165 } catch (...) {
166 error = std::current_exception();
167 }
168 }
169
170 contoura.left = DataVector<dcomplex>(const_cast<dcomplex*>(&left[0]), n+1);
171 contoura.right = DataVector<dcomplex>(const_cast<dcomplex*>(&right[0]), n+1);
172 contourb.left = DataVector<dcomplex>(const_cast<dcomplex*>(&left[n]), n+1);
173 contourb.right = DataVector<dcomplex>(const_cast<dcomplex*>(&right[n]), n+1);
174 }
175
176 contoura.bottom = bottom;
177 contoura.top = middle;
178
179 contourb.bottom = middle;
180 contourb.top = top;
181
182 return std::make_pair(std::move(contoura), std::move(contourb));
183
184 }
185}
186
187 namespace detail {
188
189 struct ContourBisect {
190 double reps, ieps;
191 std::vector<std::pair<dcomplex,dcomplex>>& results;
192 ContourBisect(double reps, double ieps, std::vector<std::pair<dcomplex,dcomplex>>& results): reps(reps), ieps(ieps), results(results) {}
193
194 int operator()(const Contour& contour) {
195 int wind = contour.winding();
196 if (wind == 0)
197 return 0;
198 if (contour.re1-contour.re0 <= reps && contour.im1-contour.im0 <= ieps) {
199 for(int i = 0; i != abs(wind); ++i)
200 results.push_back(std::make_pair(dcomplex(contour.re0, contour.im0), dcomplex(contour.re1, contour.im1)));
201 return wind;
202 }
203 auto contours = contour.divide(reps, ieps);
204 std::size_t w1 = (*this)(contours.first);
205 std::size_t w2 = (*this)(contours.second);
206 if (int(w1 + w2) < wind)
207 contour.solver->writelog(LOG_WARNING, "Lost zero between {0} and {1}", str(dcomplex(contour.re0, contour.im0)), str(dcomplex(contour.re1, contour.im1)));
208 else if (int(w1 + w2) > wind)
209 contour.solver->writelog(LOG_WARNING, "New zero between {0} and {1}", str(dcomplex(contour.re0, contour.im0)), str(dcomplex(contour.re1, contour.im1)));
210 return wind;
211 }
212 };
213}
214
215std::vector<std::pair<dcomplex,dcomplex>> findZeros(const Solver* solver, const std::function<dcomplex(dcomplex)>& fun,
216 dcomplex corner0, dcomplex corner1, size_t resteps, size_t imsteps, dcomplex eps)
217{
218 // Find first power of 2 not smaller than range/precision
219 size_t Nr = 1, Ni = 1;
220 for(; resteps > Nr; Nr <<= 1);
221 for(; imsteps > Ni; Ni <<= 1);
222
223 double reps = real(eps), ieps = imag(eps);
224
225 std::vector<std::pair<dcomplex,dcomplex>> results;
226 detail::ContourBisect bisection(reps, ieps, results);
227 Contour contour(solver, fun, corner0, corner1, Nr, Ni);
228 int zeros = abs(contour.winding());
229 solver->writelog(LOG_DETAIL, "Looking for {4} zero{5} between {0} and {1} with {2}/{3} real/imaginary intervals",
230 str(corner0), str(corner1), Nr, Ni, zeros, (zeros!=1)?"s":"");
231 bisection(contour);
232 return results;
233}
234
235
236}}} // namespace plask::optical::effective