Skip to content

Commit a82a38e

Browse files
niccolotubininiccolotubini
committed
Add the NCZ algorithm with the Conjugate Gradient method and the Matop interface
1 parent 5aaedcd commit a82a38e

File tree

2 files changed

+162
-157
lines changed

2 files changed

+162
-157
lines changed

‎src/main/java/it/geoframe/blogspot/numerical/matop/Matop.java‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
/*
22
* GNU GPL v3 License
33
*
4-
* Copyright 2019 Niccolo` Tubini
4+
* Copyright 2021 Niccolo` Tubini
55
*
66
* This program is free software: you can redistribute it and/or modify
77
* it under the terms of the GNU General Public License as published by
@@ -19,14 +19,16 @@
1919

2020
package it.geoframe.blogspot.numerical.matop;
2121

22-
import java.util.Map;
22+
import java.util.List;
2323

2424
/**
2525
*
2626
* @author Niccolo' Tubini, Riccardo Rigon, Michael Dumbser
2727
*/
2828
public abstract class Matop {
2929

30-
public abstract Map<Integer, Double> solve(Map<Integer, Double> dis, Map<Integer, Double> variable);
30+
31+
public abstract List<Double> solve(List<Double> dis, List<Double> variable);
32+
3133

3234
}

‎src/main/java/it/geoframe/blogspot/numerical/newtonalgorithm/NestedNewtonCG.java‎

Lines changed: 157 additions & 154 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,29 @@
1+
/*
2+
* GNU GPL v3 License
3+
*
4+
* Copyright 2021 Niccolo` Tubini
5+
*
6+
* This program is free software: you can redistribute it and/or modify
7+
* it under the terms of the GNU General Public License as published by
8+
* the Free Software Foundation, either version 3 of the License, or
9+
* (at your option) any later version.
10+
*
11+
* This program is distributed in the hope that it will be useful,
12+
* but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
* GNU General Public License for more details.
15+
*
16+
* You should have received a copy of the GNU General Public License
17+
* along with this program. If not, see <http://www.gnu.org/licenses/>.
18+
*/
19+
120
package it.geoframe.blogspot.numerical.newtonalgorithm;
221

3-
import java.util.HashMap;
4-
import java.util.Map;
22+
import java.util.ArrayList;
23+
import java.util.Arrays;
24+
import java.util.List;
525

6-
import bidimensionalDomain.Geometry;
7-
import bidimensionalDomain.Topology;
26+
import it.geoframe.blogspot.closureequation.equationstate.EquationState;
827
import it.geoframe.blogspot.numerical.linearsystemsolver.*;
928
import it.geoframe.blogspot.numerical.matop.*;
1029

@@ -13,34 +32,29 @@ public class NestedNewtonCG {
1332
private double outerResidual;
1433
private double innerResidual;
1534

16-
int nestedNewton;
17-
int MAXITER_NEWT;
18-
19-
double newtonTolerance;
20-
double tmp;
21-
double timeDelta;
22-
23-
24-
Map<Integer, Double> rhss;
25-
Map<Integer, Double> mainDiagonal;
26-
Map<Integer, Double> variable;
27-
28-
private Map<Integer, Double> Apsi;
29-
private Map<Integer, Double> fs;
30-
private Map<Integer, Double> fks;
31-
private Map<Integer, Double> dis;
32-
private Map<Integer, Double> dpsis;
33-
private Map<Integer, Double> psis_outer;
34-
private Map<Integer, Double> psism;
35-
36-
Topology topology;
37-
Geometry geometry;
38-
// SoilWaterRetentionCurve swrc;
39-
// //TotalDepth totalDepth;
40-
// //Thomas thomasAlg = new Thomas();
41-
//
42-
Matop matop;
43-
ConjugateGradientMethod cg;
35+
private int MAXITER_NEWT;
36+
37+
private double newtonTolerance;
38+
private double tmp;
39+
40+
41+
private List<Double> rhss;
42+
private List<Double> mainDiagonal;
43+
private List<Double> x;
44+
private List<Double> y;
45+
private List<Integer> elementEquationStateID;
46+
private List<Integer> elementParameterID;
47+
private List<EquationState> equationState;
48+
private List<Double> Apsi;
49+
private List<Double> fs;
50+
private List<Double> fks;
51+
private List<Double> dis;
52+
private List<Double> dx;
53+
private List<Double> x_outer;
54+
55+
56+
private Matop matop;
57+
private ConjugateGradientMethod cg;
4458

4559

4660

@@ -56,24 +70,25 @@ public class NestedNewtonCG {
5670
* @param thetaR vector containing the adimensional residual water contentfor each control volume, it is a vector of length NUM_CONTROL_VOLUMES-1
5771
* @param thetaS vector containing the adimensional water content at saturation for each control volume, it is a vector of length NUM_CONTROL_VOLUMES-1
5872
*/
59-
public NestedNewtonCG(int nestedNewton, double newtonTolerance, int MAXITER_NEWT, SoilWaterRetentionCurve swrc, Matop matop, double cgTolerance){
73+
public NestedNewtonCG(double newtonTolerance, int MAXITER_NEWT, List<EquationState> equationState, Matop matop, double cgTolerance,
74+
List<Integer> elementParameterID, List<Integer> elementEquationStateID) {
6075

61-
this.nestedNewton = nestedNewton;
6276
this.newtonTolerance = newtonTolerance;
6377
this.MAXITER_NEWT = MAXITER_NEWT;
64-
this.topology = Topology.getInstance();
65-
this.geometry = Geometry.getInstance();
66-
// this.swrc = swrc;
67-
// this.matop = matop;
68-
this.cg = new ConjugateGradientMethod(matop, cgTolerance);
69-
70-
71-
fs = new HashMap<Integer, Double>();
72-
fks = new HashMap<Integer, Double>();
73-
dis = new HashMap<Integer, Double>();
74-
dpsis = new HashMap<Integer, Double>();
75-
psis_outer = new HashMap<Integer, Double>();
76-
psism = new HashMap<Integer, Double>();
78+
this.equationState = equationState;
79+
this.elementParameterID = elementParameterID;
80+
this.elementEquationStateID = elementEquationStateID;
81+
this.cg = new ConjugateGradientMethod(matop, cgTolerance, elementParameterID);
82+
83+
this.matop = matop;
84+
85+
86+
fs = new ArrayList<Double>(Arrays.asList(new Double[elementParameterID.size()]));
87+
fks = new ArrayList<Double>(Arrays.asList(new Double[elementParameterID.size()]));
88+
dis = new ArrayList<Double>(Arrays.asList(new Double[elementParameterID.size()]));
89+
dx = new ArrayList<Double>(Arrays.asList(new Double[elementParameterID.size()]));
90+
x_outer = new ArrayList<Double>(Arrays.asList(new Double[elementParameterID.size()]));
91+
7792
}
7893

7994

@@ -85,143 +100,131 @@ public NestedNewtonCG(int nestedNewton, double newtonTolerance, int MAXITER_NEWT
85100
* @param lowerDiagonal lower diagonal of the coefficient matrix A of the linear system, it is a vector of length NUM_CONTROL_VOLUMES
86101
* @param rhss right hand side term of the linear system, it is a vector of length NUM_CONTROL_VOLUMES
87102
*/
88-
public void set(Map<Integer, Double> variable, Map<Integer, Double> rhss, Map<Integer, Double> mainDiagonal){
103+
public void set(List<Double> x, List<Double> y, List<Double> rhss, List<Double> mainDiagonal) {
89104

90-
this.variable = variable;
91-
this.rhss = rhss;
92-
this.mainDiagonal = mainDiagonal;
105+
this.x = new ArrayList<Double>(x);
106+
this.y = new ArrayList<Double>(y);
107+
this.rhss = new ArrayList<Double>(rhss);
108+
this.mainDiagonal = new ArrayList<Double>(mainDiagonal);
93109

94110
}
95111

96112

97113

98-
public Map<Integer, Double> solver(){
114+
public List<Double> solver() {
99115

100116

101117

102118
// Initial guess for the outer iteration
103-
for(Integer element : topology.s_i.keySet()) {
104-
tmp = Math.min(variable.get(element), SoilParameters.psiStar1[SoilParameters.elementsLabel.get(element)]);
105-
// Variables.waterSuctions.put(element, tmp);
119+
for(int element=1; element<elementParameterID.size(); element++) {
120+
tmp = equationState.get(elementEquationStateID.get(element)).initialGuess(x.get(element), elementParameterID.get(element), element);
121+
x.set(element, tmp);
106122
}
107123

108124
//// OUTER CYCLE ////
109125
for(int i = 0; i < MAXITER_NEWT; i++) {
110126
// I have to assign 0 to outerResidual otherwise I will take into account of the previous error
111127
outerResidual = 0.0;
112-
for(Integer element : topology.s_i.keySet()) {
113-
dis.put(element, 0.0);
128+
for(int element=1; element<elementParameterID.size(); element++) {
129+
130+
dis.set(element, 0.0);
131+
114132
}
115-
Apsi = matop.solve(dis, variable);
116-
// System.out.println("Apsi outer:");
117-
// for(Integer element : Topology.s_i.keySet()) {
118-
// System.out.println("\t"+ element + "\t" + Apsi.get(element));
119-
// }
120-
121-
// System.out.println("fs :");
122-
for(Integer element : topology.s_i.keySet()) {
123-
tmp = swrc.dWaterContent(variable.get(element),element)*geometry.elementsArea.get(element);
124-
dis.put(element, tmp );
125-
tmp = swrc.waterContent(variable.get(element),element)*geometry.elementsArea.get(element) - rhss.get(element) + Apsi.get(element);
126-
fs.put(element, tmp);
127-
// System.out.println(element +": " +swrc.waterContent(Variables.waterSuctions.get(element),element) +" "+Geometry.elementsArea.get(element) +" "+rhss.get(element)+" "+ Apsi.get(element)+" " +fs.get(element) + " " + dis.get(element));
128-
// System.out.println("\t"+ element + "\t" + fs.get(element));
129-
133+
134+
Apsi = matop.solve(dis, x);
135+
// System.out.println("Apsi outer:");
136+
// for(int element=1; element<elementParameterID.size(); element++) {
137+
// System.out.println("\t"+ element + "\t" + Apsi.get(element));
138+
// }
139+
//
140+
// System.out.println("fs :");
141+
for(int element=1; element<elementParameterID.size(); element++) {
142+
143+
tmp = equationState.get(elementEquationStateID.get(element)).dEquationState(x.get(element), y.get(element), elementParameterID.get(element), element);
144+
dis.set(element, tmp);
145+
146+
tmp = equationState.get(elementEquationStateID.get(element)).equationState(x.get(element), y.get(element), elementParameterID.get(element), element) - rhss.get(element) + Apsi.get(element);
147+
fs.set(element, tmp);
148+
// System.out.println("\t"+ element + "\t" + fs.get(element));
149+
130150
outerResidual += tmp*tmp;
131151
}
152+
132153
outerResidual = Math.pow(outerResidual,0.5);
133-
System.out.println("\t\t-Outer iteration " + i + " with residual " + outerResidual);
154+
// System.out.println("\t\t-Outer iteration " + i + " with residual " + outerResidual);
134155
if(outerResidual < newtonTolerance) {
156+
135157
break;
158+
136159
}
137-
if(nestedNewton == 0){
138160

139-
}else{
140161

141-
// Initial guess for the inner iteration (optional)
142-
for(Integer element : topology.s_i.keySet()) {
143-
// psis_outer.put(element, Variables.waterSuctions.get(element));
144-
// Variables.waterSuctions.put(element, Math.max(Variables.waterSuctions.get(element), SoilParameters.psiStar1[SoilParameters.elementsLabel.get(element)]));
145-
}
162+
for(int element=1; element<elementParameterID.size(); element++) {
146163

147-
//// INNER CYCLE ////
148-
for(int j = 0; j < MAXITER_NEWT; j++) {
149-
// I have to assign 0 to innerResidual otherwise I will take into account of the previous error
150-
innerResidual = 0.0;
151-
for(Integer element : topology.s_i.keySet()) {
152-
dis.put(element, 0.0);
153-
}
154-
Apsi = matop.solve(dis, variable);
155-
// System.out.println("Apsi inner:");
156-
// for(Integer element : Topology.s_i.keySet()) {
157-
// System.out.println("\t"+ element + "\t" + Apsi.get(element));
158-
// }
159-
// System.out.println("inner psi-psi_outer :");
160-
for(Integer element : topology.s_i.keySet()) {
161-
162-
tmp = (swrc.p(variable.get(element),element) - swrc.q(psis_outer.get(element),element))*geometry.elementsArea.get(element);
163-
dis.put(element, tmp);
164-
tmp = swrc.pIntegral(variable.get(element),element)*geometry.elementsArea.get(element)
165-
- ( swrc.qIntegral(psis_outer.get(element),element) + swrc.q(psis_outer.get(element),element)*(variable.get(element) - psis_outer.get(element)) )*Geometry.elementsArea.get(element)
166-
- rhss.get(element) + Apsi.get(element);
167-
fks.put(element, tmp);
168-
169-
// System.out.println(element +" "+Variables.waterSuctions.get(element) + " " + psis_outer.get(element));
170-
// System.out.println("\t"+ element + "\t" + fks.get(element));
171-
172-
innerResidual += tmp*tmp;
173-
}
174-
175-
innerResidual = Math.pow(innerResidual,0.5);
176-
System.out.println("\t\t\t-Inner iteration " + j + " with residual " + innerResidual);
177-
if(innerResidual < newtonTolerance) {
178-
// System.out.println("Psi:");
179-
// for(Integer element : Topology.s_i.keySet()) {
180-
// System.out.println("\t"+ element + "\t" + Variables.waterSuctions.get(element));
181-
//
182-
// }
183-
//
184-
// System.out.println("\n\n");
185-
// System.out.println("dpsi:");
186-
// for(Integer element : Topology.s_i.keySet()) {
187-
// System.out.println("\t"+ element + "\t" + dpsis.get(element));
188-
//
189-
// }
190-
//
191-
break;
192-
}
193-
//// CONJUGATE GRADIENT METHOD////
194-
dpsis = cg.solve(dis, fks, mainDiagonal);
195-
// dpsis = cg.solve(dis, fks);
196-
197-
for(Integer element : topology.s_i.keySet()) {
198-
// psism.put(element, Variables.waterSuctions.get(element));
199-
tmp = variable.get(element)-dpsis.get(element);
200-
variable.put(element, tmp );
201-
}
202-
203-
// if(j>1) {
204-
// for(Integer element : Topology.s_i.keySet()) {
205-
// Variables.waterSuctions.put(element, Math.min( Variables.waterSuctions.get(element), psism.get(element)) );
206-
// if(i>1) {
207-
// Variables.waterSuctions.put(element, Math.max( Variables.waterSuctions.get(element), psis_outer.get(element)) );
208-
// }
209-
// }
210-
// }
211-
} //// INNER CYCLE END ////
164+
x_outer.set(element, x.get(element));
212165

213166
}
214-
// return variable;
215-
// if(i>1) {
216-
// for(Integer element : Topology.s_i.keySet()) {
217-
// Variables.waterSuctions.put(element, Math.max( Variables.waterSuctions.get(element), psis_outer.get(element)) );
218-
// }
219-
// }
220-
return variable;
221-
} //// OUTER CYCLE END ////
222167

223-
}
168+
//// INNER CYCLE ////
169+
for(int j = 0; j < MAXITER_NEWT; j++) {
170+
// I have to assign 0 to innerResidual otherwise I will take into account of the previous error
171+
innerResidual = 0.0;
172+
for(int element=1; element<elementParameterID.size(); element++) {
173+
dis.set(element, 0.0);
174+
}
175+
Apsi = matop.solve(dis, x);
176+
// System.out.println("Apsi inner:");
177+
// for(int element=1; element<elementParameterID.size(); element++) {
178+
// System.out.println("\t"+ element + "\t" + Apsi.get(element));
179+
// }
180+
//
181+
// System.out.println("fks :");
182+
for(int element=1; element<elementParameterID.size(); element++) {
183+
184+
tmp = equationState.get(elementEquationStateID.get(element)).p(x.get(element), y.get(element), elementParameterID.get(element), element) - equationState.get(elementEquationStateID.get(element)).q(x_outer.get(element), y.get(element), elementParameterID.get(element), element);
185+
// System.out.println("\t\telement "+element+"\t"+x.get(element)+"\t"+equationState.get(elementEquationStateID.get(element)).p(x.get(element), y.get(element), elementParameterID.get(element), element)+"\t"+equationState.get(elementEquationStateID.get(element)).q(x_outer.get(element), y.get(element), elementParameterID.get(element), element));
186+
dis.set(element, tmp);
187+
188+
tmp = equationState.get(elementEquationStateID.get(element)).pIntegral(x.get(element), y.get(element), elementParameterID.get(element), element)
189+
- ( equationState.get(elementEquationStateID.get(element)).qIntegral(x_outer.get(element), y.get(element), elementParameterID.get(element), element) + equationState.get(elementEquationStateID.get(element)).q(x_outer.get(element), y.get(element), elementParameterID.get(element), element)*(x.get(element) - x_outer.get(element)) )
190+
- rhss.get(element) + Apsi.get(element);
191+
fks.set(element, tmp);
192+
// System.out.println("\t"+ element + "\t" + fks.get(element));
193+
194+
innerResidual += tmp*tmp;
195+
}
196+
197+
innerResidual = Math.pow(innerResidual,0.5);
198+
// System.out.println("\t\t\t-Inner iteration " + j + " with residual " + innerResidual);
199+
if(innerResidual < newtonTolerance) {
200+
201+
break;
224202

203+
}
204+
205+
206+
/*
207+
* CONJUGATE GRADIENT METHOD
208+
*/
209+
dx = cg.solve(dis, fks, mainDiagonal);
210+
// System.out.println("CG");
211+
// System.out.println("dx psi");
212+
for(int element=1; element<elementParameterID.size(); element++) {
213+
214+
tmp = x.get(element)-dx.get(element);
215+
x.set(element, tmp);
216+
// System.out.println(dx.get(element) +" "+ tmp);
217+
218+
}
219+
220+
} //// OUTER CYCLE END ////
221+
222+
}
223+
224+
return x;
225+
226+
227+
}
225228

226229
}
227230

0 commit comments

Comments
 (0)