Skip to content

Commit b317cd0

Browse files
niccolotubininiccolotubini
committed
Update
1 parent 316e6ec commit b317cd0

1 file changed

Lines changed: 42 additions & 39 deletions

File tree

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

Lines changed: 42 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -55,7 +55,7 @@ public class NestedNewtonThomas {
5555
private double[] bb;
5656
private double[] cc;
5757
private double[] dis;
58-
private int[] rheologyID;
58+
private int[] equationStateID;
5959
private int[] parameterID;
6060
private double[] x_outer;
6161

@@ -89,7 +89,7 @@ public NestedNewtonThomas(int nestedNewton, double newtonTolerance, int MAXITER_
8989

9090

9191

92-
public void set(double[] x, double[] y, double[] mainDiagonal, double[] upperDiagonal, double[] lowerDiagonal, double[] rhss, int KMAX, int[] parameterID, int[] rheologyID){
92+
public void set(double[] x, double[] y, double[] mainDiagonal, double[] upperDiagonal, double[] lowerDiagonal, double[] rhss, int KMAX, int[] parameterID, int[] equationStateID){
9393

9494
this.x = x;
9595
this.y = y;
@@ -101,7 +101,7 @@ public void set(double[] x, double[] y, double[] mainDiagonal, double[] upperDia
101101
this.KMAX = KMAX;
102102

103103
this.parameterID = parameterID;
104-
this.rheologyID = rheologyID;
104+
this.equationStateID = equationStateID;
105105

106106
}
107107

@@ -117,7 +117,7 @@ public double[] solver(){
117117
} else {
118118
for(int element = 0; element < KMAX; element++) {
119119
//x[element] = Math.min(x[element], xStar[element]-1 );
120-
x[element] = equationState.get(rheologyID[element]).initialGuess(x[element],parameterID[element],element);
120+
x[element] = equationState.get(equationStateID[element]).initialGuess(x[element],parameterID[element],element);
121121
}
122122
}
123123

@@ -127,25 +127,27 @@ public double[] solver(){
127127
for(int i = 0; i < MAXITER_NEWT; i++) {
128128
// I have to assign 0 to outerResidual otherwise I will take into account of the previous error
129129
outerResidual = 0.0;
130-
for(int element = 0; element < KMAX; element++) {
131-
if(element==0) {
132-
fs[element] = equationState.get(rheologyID[element]).equationState(x[element],y[element],parameterID[element],element) - rhss[element] + mainDiagonal[element]*x[element] + upperDiagonal[element]*x[element+1];
133-
dis[element] = equationState.get(rheologyID[element]).dEquationState(x[element],y[element],parameterID[element],element);
134-
// System.out.println(element+" "+fs[element]);
135-
} else if(element==KMAX-1) {
136-
fs[element] = equationState.get(rheologyID[element]).equationState(x[element],y[element],parameterID[element],element) - rhss[element] + lowerDiagonal[element]*x[element-1] + mainDiagonal[element]*x[element];
137-
dis[element] = equationState.get(rheologyID[element]).dEquationState(x[element],y[element],parameterID[element],element);
138-
// System.out.println(element+" "+fs[element]);
139-
} else {
140-
fs[element] = equationState.get(rheologyID[element]).equationState(x[element],y[element],parameterID[element],element) - rhss[element] + lowerDiagonal[element]*x[element-1] + mainDiagonal[element]*x[element] + upperDiagonal[element]*x[element+1];
141-
dis[element] = equationState.get(rheologyID[element]).dEquationState(x[element],y[element],parameterID[element],element);
142-
// System.out.println(element+" "+fs[element]);
143-
}
144-
130+
for(int element = 1; element < KMAX-1; element++) {
131+
fs[element] = equationState.get(equationStateID[element]).equationState(x[element],y[element],parameterID[element],element) - rhss[element] + lowerDiagonal[element]*x[element-1] + mainDiagonal[element]*x[element] + upperDiagonal[element]*x[element+1];
132+
dis[element] = equationState.get(equationStateID[element]).dEquationState(x[element],y[element],parameterID[element],element);
133+
// System.out.println(element+" "+fs[element]);
145134
outerResidual += fs[element]*fs[element];
146135
}
136+
// element==0
137+
fs[0] = equationState.get(equationStateID[0]).equationState(x[0],y[0],parameterID[0],0) - rhss[0] + mainDiagonal[0]*x[0] + upperDiagonal[0]*x[1];
138+
dis[0] = equationState.get(equationStateID[0]).dEquationState(x[0],y[0],parameterID[0],0);
139+
// System.out.println("0 "+fs[0]);
140+
outerResidual += fs[0]*fs[0];
141+
142+
// element==KMAX-1
143+
fs[KMAX-1] = equationState.get(equationStateID[KMAX-1]).equationState(x[KMAX-1],y[KMAX-1],parameterID[KMAX-1],KMAX-1) - rhss[KMAX-1] + lowerDiagonal[KMAX-1]*x[KMAX-2] + mainDiagonal[KMAX-1]*x[KMAX-1];
144+
dis[KMAX-1] = equationState.get(equationStateID[KMAX-1]).dEquationState(x[KMAX-1],y[KMAX-1],parameterID[KMAX-1],KMAX-1);
145+
// System.out.println("KMAX-1 "+fs[element]);
146+
outerResidual += fs[KMAX-1]*fs[KMAX-1];
147+
148+
147149
outerResidual = Math.pow(outerResidual,0.5);
148-
// System.out.println("\tOuter iteration " + i + " with residual " + outerResidual);
150+
System.out.println("\tOuter iteration " + i + " with residual " + outerResidual);
149151
if(outerResidual < newtonTolerance) {
150152
break;
151153
}
@@ -176,29 +178,30 @@ public double[] solver(){
176178
for(int j = 0; j < MAXITER_NEWT; j++) {
177179
// I have to assign 0 to innerResidual otherwise I will take into account of the previous error
178180
innerResidual = 0.0;
179-
for(int element=0; element < KMAX; element++) {
180-
if(element==0) {
181-
fks[element] = equationState.get(rheologyID[element]).pIntegral(x[element],y[element],parameterID[element],element) - ( equationState.get(rheologyID[element]).qIntegral(x_outer[element],y[element],parameterID[element],element) + equationState.get(rheologyID[element]).q(x_outer[element],y[element],parameterID[element],element)*(x[element]-x_outer[element]) ) - this.rhss[element] + mainDiagonal[element]*x[element] + upperDiagonal[element]*x[element+1];
182-
dis[element] = ( equationState.get(rheologyID[element]).p(x[element],y[element],parameterID[element],element) - equationState.get(rheologyID[element]).q(x_outer[element],y[element],parameterID[element],element) );
183-
// System.out.println(element+" "+fks[element]);
184-
// System.out.println(element+" "+dis[element]);
185-
} else if(element==KMAX-1) {
186-
fks[element] = equationState.get(rheologyID[element]).pIntegral(x[element],y[element],parameterID[element],element) - ( equationState.get(rheologyID[element]).qIntegral(x_outer[element],y[element],parameterID[element],element) + equationState.get(rheologyID[element]).q(x_outer[element],y[element],parameterID[element],element)*(x[element]-x_outer[element]) ) - this.rhss[element] + lowerDiagonal[element]*x[element-1] + mainDiagonal[element]*x[element];
187-
dis[element] = ( equationState.get(rheologyID[element]).p(x[element],y[element],parameterID[element],element) - equationState.get(rheologyID[element]).q(x_outer[element],y[element],parameterID[element],element) );
188-
// System.out.println(element+" "+fks[element]);
189-
// System.out.println(element+" "+dis[element]);
190-
} else {
191-
fks[element] = equationState.get(rheologyID[element]).pIntegral(x[element],y[element],parameterID[element],element) - ( equationState.get(rheologyID[element]).qIntegral(x_outer[element],y[element],parameterID[element],element) + equationState.get(rheologyID[element]).q(x_outer[element],y[element],parameterID[element],element)*(x[element]-x_outer[element]) ) - this.rhss[element] + lowerDiagonal[element]*x[element-1] + mainDiagonal[element]*x[element] + upperDiagonal[element]*x[element+1];
192-
dis[element] = ( equationState.get(rheologyID[element]).p(x[element],y[element],parameterID[element],element) - equationState.get(rheologyID[element]).q(x_outer[element],y[element],parameterID[element],element) );
193-
// System.out.println(element+" "+fks[element]);
194-
// System.out.println(element+" "+dis[element]);
195-
}
196-
181+
for(int element=1; element < KMAX-1; element++) {
182+
fks[element] = equationState.get(equationStateID[element]).pIntegral(x[element],y[element],parameterID[element],element) - ( equationState.get(equationStateID[element]).qIntegral(x_outer[element],y[element],parameterID[element],element) + equationState.get(equationStateID[element]).q(x_outer[element],y[element],parameterID[element],element)*(x[element]-x_outer[element]) ) - this.rhss[element] + lowerDiagonal[element]*x[element-1] + mainDiagonal[element]*x[element] + upperDiagonal[element]*x[element+1];
183+
dis[element] = ( equationState.get(equationStateID[element]).p(x[element],y[element],parameterID[element],element) - equationState.get(equationStateID[element]).q(x_outer[element],y[element],parameterID[element],element) );
184+
// System.out.println(element+" "+fks[element]);
185+
// System.out.println(element+" "+dis[element]);
197186
innerResidual += fks[element]*fks[element];
198187
}
188+
// element==0
189+
fks[0] = equationState.get(equationStateID[0]).pIntegral(x[0],y[0],parameterID[0],0) - ( equationState.get(equationStateID[0]).qIntegral(x_outer[0],y[0],parameterID[0],0) + equationState.get(equationStateID[0]).q(x_outer[0],y[0],parameterID[0],0)*(x[0]-x_outer[0]) ) - this.rhss[0] + mainDiagonal[0]*x[0] + upperDiagonal[0]*x[1];
190+
dis[0] = ( equationState.get(equationStateID[0]).p(x[0],y[0],parameterID[0],0) - equationState.get(equationStateID[0]).q(x_outer[0],y[0],parameterID[0],0) );
191+
// System.out.println("0 "+fks[element]);
192+
// System.out.println("0 "+dis[element]);
193+
innerResidual += fks[0]*fks[0];
194+
195+
// element==KMAX-1
196+
fks[KMAX-1] = equationState.get(equationStateID[KMAX-1]).pIntegral(x[KMAX-1],y[KMAX-1],parameterID[KMAX-1],KMAX-1) - ( equationState.get(equationStateID[KMAX-1]).qIntegral(x_outer[KMAX-1],y[KMAX-1],parameterID[KMAX-1],KMAX-1) + equationState.get(equationStateID[KMAX-1]).q(x_outer[KMAX-1],y[KMAX-1],parameterID[KMAX-1],KMAX-1)*(x[KMAX-1]-x_outer[KMAX-1]) ) - this.rhss[KMAX-1] + lowerDiagonal[KMAX-1]*x[KMAX-2] + mainDiagonal[KMAX-1]*x[KMAX-1];
197+
dis[KMAX-1] = ( equationState.get(equationStateID[KMAX-1]).p(x[KMAX-1],y[KMAX-1],parameterID[KMAX-1],KMAX-1) - equationState.get(equationStateID[KMAX-1]).q(x_outer[KMAX-1],y[KMAX-1],parameterID[KMAX-1],KMAX-1) );
198+
// System.out.println("KMAX-1 "+fks[KMAX-1]);
199+
// System.out.println("KMAX-1 "+dis[KMAX-1]);
200+
innerResidual += fks[KMAX-1]*fks[KMAX-1];
201+
199202
innerResidual = Math.pow(innerResidual,0.5);
200203

201-
// System.out.println("\t\t-Inner iteration " + j + " with residual " + innerResidual);
204+
System.out.println("\t\t-Inner iteration " + j + " with residual " + innerResidual);
202205

203206
if(innerResidual < newtonTolerance) {
204207
break;

0 commit comments

Comments
 (0)