- Thank you received: 29
DG and Newtons method linearization
6 years 10 months ago #330
by cwinters
Replied by cwinters on topic DG and Newtons method linearization
Hi Bart,
HDG was not fully implemented for 1D. An updated version as available on Github and you can pull it from there. If you use an installer, you can update/download ngsolve_nightly tomorrow.
If you want to use HDG, you have to use an discontinuous finite element space for your element variable (e.g. L2). For a 1d mesh, you just have one degree of freedom on the facets, thus you have a 0th order facet finite element space.
A full example solving a laplace problem with HDG is attached.
Best,
Christoph
HDG was not fully implemented for 1D. An updated version as available on Github and you can pull it from there. If you use an installer, you can update/download ngsolve_nightly tomorrow.
If you want to use HDG, you have to use an discontinuous finite element space for your element variable (e.g. L2). For a 1d mesh, you just have one degree of freedom on the facets, thus you have a 0th order facet finite element space.
Code:
order = 4
V1 = L2(mesh, order=order)
V2 = FacetFESpace(mesh, order=0, dirichlet=[1])
V = FESpace([V1,V2])
A full example solving a laplace problem with HDG is attached.
Best,
Christoph
Attachments:
The following user(s) said Thank You: Matejczyk
6 years 9 months ago - 6 years 9 months ago #337
by Matejczyk
Replied by Matejczyk on topic DG and Newtons method linearization
Hi,
Thanks for the previous answer. I defined my spaces as you advised and now they look like
.
As always there is a but. Now I struggle with the definitions of the BC and the initial guess for the Newton loop.
What I have is
which (as I believe) provides the right boundary conditions for the facet space and should provide the BC for the Newton loop. Unfortunately this way the code gets stuck on the line
I fixed this issue and it is not stuck any more. I needed to increse the heapsize by seting
Now I end up with a singular matrix so I believe the BC/initialization problem stays the same.
Is there a better way to implement the BC and the initial guess than the one I'm using?
Best wishes,
Bart
Thanks for the previous answer. I defined my spaces as you advised and now they look like
Code:
Y = L2(mesh, order=ord) # space can be fixed for the whole routine
Z = FacetFESpace(mesh, order=0,dirichlet = [1,2])
X = FESpace([Y,Z,Y,Z,Y,Z])
As always there is a but. Now I struggle with the definitions of the BC and the initial guess for the Newton loop.
What I have is
Code:
u.components[1].vec[nel] = Vapl
u.components[3].vec[nel]=z1*beta*Vapl + log(bc1)
u.components[5].vec[nel]=z2*beta*Vapl + log(bc1)
Code:
a.Apply(u.vec, res)
Code:
a.Apply(u.vec, res, heapsize=100000000)
Is there a better way to implement the BC and the initial guess than the one I'm using?
Best wishes,
Bart
Last edit: 6 years 9 months ago by Matejczyk.
6 years 9 months ago #338
by cwinters
Replied by cwinters on topic DG and Newtons method linearization
Hi,
for the 1d case, your way of setting the boundary conditions should work. (The Set() function which you used before does not support points by now.)
Which matrix is singular? The result of "a.Apply(...)" is a vector.
If you do not set anything else than those three boundary values, all other values in "u.vec" are 0. Is that the vector you want to start with?
Best,
Christoph
for the 1d case, your way of setting the boundary conditions should work. (The Set() function which you used before does not support points by now.)
Which matrix is singular? The result of "a.Apply(...)" is a vector.
If you do not set anything else than those three boundary values, all other values in "u.vec" are 0. Is that the vector you want to start with?
Best,
Christoph
6 years 9 months ago #339
by Matejczyk
Replied by Matejczyk on topic DG and Newtons method linearization
Hi,
I thought that I'm doing something wrong with the Set function. Is good to know that it is not only me.
For the BC and the initial Guess.
I set 3 boundary conditions manually by
so that the facet space vector is 0 everywhere but the boundary. I don't set anything for the values of the functions in the L2 space.
Then my system enters the Newtons loop and what I get is
[...]
linear form done
Iteration 0
<A u 0 , A u 0 >_{-1}^0.5 = 190.3901437327294
Iteration 1
<A u 1 , A u 1 >_{-1}^0.5 = 2.0534197662160665
UMFPACK V5.7.4 (Feb 1, 2016): WARNING: matrix is singular
[...]
so I get a singular matrix in my Newton loop which looks like
. The problematic line is the inverse now. Any ideas why the HDG formulation might lead to a singular matrix? I was solving the same system just in H1 space and I didn't have this problem.
Best wishes,
B
I thought that I'm doing something wrong with the Set function. Is good to know that it is not only me.
For the BC and the initial Guess.
I set 3 boundary conditions manually by
Code:
u.components[1].vec[nel]=
Then my system enters the Newtons loop and what I get is
[...]
linear form done
Iteration 0
<A u 0 , A u 0 >_{-1}^0.5 = 190.3901437327294
Iteration 1
<A u 1 , A u 1 >_{-1}^0.5 = 2.0534197662160665
UMFPACK V5.7.4 (Feb 1, 2016): WARNING: matrix is singular
[...]
so I get a singular matrix in my Newton loop which looks like
Code:
res = u.vec.CreateVector()
du = u.vec.CreateVector()
SetHeapSize(100000000)
################################# newton loop
for it in range(250):
a.Apply(u.vec, res)
a.AssembleLinearization(u.vec)
du.data = a.mat.Inverse(X.FreeDofs()) * res
u.vec.data -= du
#stopping criteria
stopcritval = sqrt(abs(InnerProduct(du,res)))
print ("Iteration",it)
print ("<A u",it,", A u",it,">_{-1}^0.5 = ", stopcritval)
Best wishes,
B
6 years 9 months ago #340
by cwinters
Replied by cwinters on topic DG and Newtons method linearization
Hi,
could you please attach your updated file with the HDG formulation.
Best,
Christoph
could you please attach your updated file with the HDG formulation.
Best,
Christoph
6 years 9 months ago #341
by Matejczyk
Replied by Matejczyk on topic DG and Newtons method linearization
Hi,
The problem with the singular matrices arises because in my bilinear form I have some functions that are x dependent A(x) and Sigma(x) (you can see the code attached). Before I was just assigning values to the vectors manually using Sigma.vec = g(i) . As now Sigma is a GridFunction(L2) this way is a no longer valid method of assigning values. I changed it then and in the code, I attach I just set
to have something. Now the Newton loop converges but my ability to assign values is somehow limited. We are somehow back to the Set function we discussed before. Is there any way to assign the values of this functions (A and sigma) such that I can use for example indicator functions over x?
Best wishes,
Bart
The problem with the singular matrices arises because in my bilinear form I have some functions that are x dependent A(x) and Sigma(x) (you can see the code attached). Before I was just assigning values to the vectors manually using Sigma.vec = g(i) . As now Sigma is a GridFunction(L2) this way is a no longer valid method of assigning values. I changed it then and in the code, I attach I just set
Code:
Sigma.Set(x)
Best wishes,
Bart
Attachments:
Time to create page: 0.135 seconds