Skip to content
Snippets Groups Projects
Commit a23667f9 authored by Nam Nguyen's avatar Nam Nguyen
Browse files

Add FE implementations and tests for HMHF2D and RSHMHF

parent 079781c7
No related branches found
No related tags found
No related merge requests found
Showing
with 1529 additions and 2 deletions
from netgen.geom2d import SplineGeometry
import ngsolve as ngs
from .parameters_HMHF2D import *
def Create2DMesh(space_params):
domain = space_params.attributes.get("domain")
use_dirichlet = space_params.attributes.get("Dirichlet")
dirichlet_bnd_name = ""
if use_dirichlet:
dirichlet_bnd_name = space_params.attributes.get("Dirichlet_bnd_name")
if domain != "unit_disk":
raise RuntimeError(" Domain type ", domain, " is not supported! For now only unit_disk!")
geo = SplineGeometry()
geo.AddCircle((0,0),1,bc=dirichlet_bnd_name)
mesh_size = space_params.attributes.get("mesh_size")
fes_order = space_params.attributes.get("fes_order")
mesh = ngs.Mesh(geo.GenerateMesh(maxh=mesh_size))
#print(dirichlet_bnd_name)
## use isoparametric elements for higher order
if fes_order > 1:
#mesh.Curve(fes_order)
mesh.ngmesh.SecondOrder()
return mesh
import math
import sys
import fileinput
class ParametersHMHF2D:
def __init__(self):
self.attributes = dict([])
def ReadParametersFromFile(self,filename=""):
if len(sys.argv) < 2 and filename == "":
raise RuntimeError("No input file!")
if filename == "":
filename = sys.argv[1]
for line in fileinput.input(files=filename):
#print("line ", fileinput.lineno())
if line[0]=="#":
continue
arr_string = line.split()
if len(arr_string) == 0:
#print("Nothing in here")
continue
if len(arr_string) < 2:
raise RuntimeError('In line {} of the file "{}" an argument is missing!'.format(fileinput.lineno(),fileinput.filename()))
if len(arr_string)==2:
if arr_string[0] in self.attributes:
self.attributes[arr_string[0]]= arr_string[1]
def PrintParameters(self):
attributes_list = list(self.attributes)
for iterator in attributes_list:
print(iterator, " = ", self.attributes.get(iterator))
class TimeInitialParametersHMHF2D(ParametersHMHF2D):
def __init__(self):
ParametersHMHF2D.__init__(self)
self.attributes.update({"tstart":0.0, "tend":1.0,\
"tstep":0.1,"BDF":1, "solver_type":"PPFEM"})
def ReadParametersFromFile(self,filename=""):
ParametersHMHF2D.ReadParametersFromFile(self,filename)
self.attributes["tstart"] = float(self.attributes.get("tstart"))
self.attributes["tend"] = float(self.attributes.get("tend"))
self.attributes["tstep"] = float(self.attributes.get("tstep"))
self.attributes["BDF"] = int(self.attributes.get("BDF"))
self.attributes["solver_type"] = self.attributes.get("solver_type")
def PrintParameters(self):
ParametersHMHF2D.PrintParameters(self)
class SpaceInitialParametersHMHF2D(ParametersHMHF2D):
def __init__(self):
ParametersHMHF2D.__init__(self)
self.attributes.update({"domain":"unit_disk",\
"mesh_size":0.1,\
"Dirichlet":1, "Dirichlet_bnd_name":"circle",\
"fes_order":1, "CPFEM_max_iter":100, "CPFEM_eps":1e-10})
def ReadParametersFromFile(self,filename=""):
ParametersHMHF2D.ReadParametersFromFile(self,filename)
self.attributes["domain"] = self.attributes.get("domain")
self.attributes["mesh_size"] = float(self.attributes.get("mesh_size"))
self.attributes["Dirichlet"] = int(self.attributes.get("Dirichlet"))
self.attributes["Dirichlet_bnd_name"] = self.attributes.get("Dirichlet_bnd_name")
self.attributes["fes_order"] = int(self.attributes.get("fes_order"))
### parameter for fix point iteration of CPFEM
self.attributes["CPFEM_max_iter"] = int(self.attributes.get("CPFEM_max_iter"))
self.attributes["CPFEM_eps"] = float(self.attributes.get("CPFEM_eps"))
def PrintParameters(self):
ParametersHMHF2D.PrintParameters(self)
### ToDo: Add more output options
class OutputParametersHMHF2D(ParametersHMHF2D):
def __init__(self):
ParametersHMHF2D.__init__(self)
self.attributes.update({"debug": 1, "save_mesh":"mesh.dat", \
"save_gfu":"gfu.dat", "save_energy":"energy.dat", \
"save_time":"time.dat", "save_tstep":"tstep.dat", \
"save_ndof":"ndof.dat", \
"save_temporal_error_indicator":"eta_tau.dat", \
"save_spacial_error_indicator":"eta_h.dat", \
"save":"0","vtk":0, "vtkFilename":"", "vtkInterval":100
})
def ReadParametersFromFile(self,filename=""):
ParametersHMHF2D.ReadParametersFromFile(self,filename)
self.attributes["debug"] = int(self.attributes.get("debug"))
self.attributes["save"] = int(self.attributes.get("save"))
self.attributes["vtk"] = int(self.attributes.get("vtk"))
self.attributes["vtkFilename"] = self.attributes.get("vtkFilename")
self.attributes["vtkInterval"] = int(self.attributes.get("vtkInterval"))
def PrintParameters(self):
ParametersHMHF2D.PrintParameters(self)
class AdaptivityParametersHMHF2D(ParametersHMHF2D):
def __init__(self):
ParametersHMHF2D.__init__(self)
self.attributes.update({"use_adaptivity":0,"tstep_min":1e-8, \
"tstep_max":1e-1, "alpha_minus":0.1, \
"alpha_plus":2.0, "eps_err":5e-1, \
"eps_mesh_size":1e-12, "theta":0.5, \
"space_adaptive_option":"", \
"time_adaptive_option":"", \
"marking_strategy":"Doerfler", \
"blowup_rate":0
})
def ReadParametersFromFile(self,filename=""):
ParametersHMHF2D.ReadParametersFromFile(self,filename)
self.attributes["use_adaptivity"] = int(self.attributes.get("use_adaptivity"))
self.attributes["tstep_min"] = float(self.attributes.get("tstep_min"))
self.attributes["tstep_max"] = float(self.attributes.get("tstep_max"))
self.attributes["alpha_minus"] = float(self.attributes.get("alpha_minus"))
self.attributes["alpha_plus"] = float(self.attributes.get("alpha_plus"))
self.attributes["eps_err"] = float(self.attributes.get("eps_err"))
self.attributes["eps_mesh_size"] = float(self.attributes.get("eps_mesh_size"))
self.attributes["theta"] = float(self.attributes.get("theta"))
def PrintParameters(self):
ParametersHMHF2D.PrintParameters(self)
def ReadAndPrintAllParametersFromFile(filename="",debug = 1):
tparams = TimeInitialParametersHMHF2D()
tparams.ReadParametersFromFile(filename)
sparams = SpaceInitialParametersHMHF2D()
sparams.ReadParametersFromFile(filename)
oparams = OutputParametersHMHF2D()
oparams.ReadParametersFromFile(filename)
aparams = AdaptivityParametersHMHF2D()
aparams.ReadParametersFromFile(filename)
if debug:
tparams.PrintParameters()
sparams.PrintParameters()
oparams.PrintParameters()
aparams.PrintParameters()
params = {"time":tparams,"space":sparams,"output":oparams,"adaptivity":aparams}
return params
from ngsolve.solvers import *
from ngsolve.comp import IntegrationRuleSpace
from .parameters_HMHF2D import *
from .meshing_2D import *
import sys
import numpy as np
class SolvingHMHF2D:
def __init__(self,initial_cond, params):
self.time_init_params = params.get("time")
self.space_init_params = params.get("space")
self.output_params = params.get("output")
self.adaptivity_params = params.get("adaptivity")
self.mesh = Create2DMesh(self.space_init_params)
fes_order = self.space_init_params.attributes.get("fes_order")
dirichlet_bnd_name = self.space_init_params.attributes.get("Dirichlet_bnd_name")
self.solver_type_name = self.time_init_params.attributes.get("solver_type")
#print("Using ", self.solver_type_name, " to solve HMHF2D!")
if self.adaptivity_params.attributes.get("use_adaptivity"):
raise RuntimeError("No adaptivity implemented yet!")
self.fes_component = ngs.H1(self.mesh, order=fes_order,dirichlet=dirichlet_bnd_name,autoupdate=True)
else:
self.fes_component = ngs.H1(self.mesh, order=fes_order,dirichlet=dirichlet_bnd_name)
if self.solver_type_name == "TFEM":
self.fes = ngs.FESpace([self.fes_component,self.fes_component,self.fes_component,self.fes_component])
else:
#self.fes = self.fes_component**3
self.fes = ngs.H1(self.mesh, order=fes_order,dirichlet=dirichlet_bnd_name, dim=3)
#print("ndof = ", self.fes.ndof)
self.initial_cond = initial_cond
self.gfu = [] # store numerical solution at current timestep j [0], previous j-1 [1], j-2 [2] and so on depending on how many steps needed in bdf scheme
if self.solver_type_name != "TFEM":
for i in range(self.time_init_params.attributes.get("BDF")+1):
self.gfu.append(ngs.GridFunction(self.fes))
self.gfu[i].Set(self.initial_cond)
else:
for i in range(self.time_init_params.attributes.get("BDF")+1):
self.gfu.append(ngs.GridFunction(self.fes))
self.gfu[i].components[0].Set(self.initial_cond[0])
self.gfu[i].components[1].Set(self.initial_cond[1])
self.gfu[i].components[2].Set(self.initial_cond[2])
## discrete time derivative
self.gfu_dot = ngs.GridFunction(self.fes)
#print("number of gfus = ",len(self.gfu), " for bdf ",self.time_init_params.attributes.get("BDF"))
self.gfu_extrap = ngs.GridFunction(self.fes)
if self.solver_type_name == "TFEM":
self.gfu_extrap.components[0].Set(self.initial_cond[0])
self.gfu_extrap.components[1].Set(self.initial_cond[1])
self.gfu_extrap.components[2].Set(self.initial_cond[2])
tstep =self.time_init_params.attributes.get("tstep")
self.dt = ngs.Parameter(tstep)
self.bf = ngs.BilinearForm(self.fes)
self.lf = ngs.LinearForm(self.fes)
### catch some cases
if not ( self.solver_type_name == "PPFEM" or self.solver_type_name == "TFEM" or self.solver_type_name == "CPFEM"):
raise RuntimeError("Solver type ",self.solver_type_name, " not supported! Only PPFEM, TFEM and CPFEM" )
if self.time_init_params.attributes.get("BDF") > 2 or self.time_init_params.attributes.get("BDF") < 1 :
raise RuntimeError("Only BDF 1 and 2 are supported")
if self.space_init_params.attributes.get("fes_order") > 2 or self.space_init_params.attributes.get("fes_order") < 1:
raise RuntimeError("Only P1 and P2 are supported")
if self.solver_type_name == "CPFEM":
if self.time_init_params.attributes.get("BDF") != 1:
raise RuntimeError("CPFEM is not supported by BDF with order higher than 1")
if fes_order != 1:
raise RuntimeError("CPFEM is not supported by finite element space with order higher than 1")
### Initialize VTK output
if self.solver_type_name != "TFEM":
self.vtk = ngs.VTKOutput(ma = self.mesh,coefs = [self.gfu[0]],names=["gfu"],filename=self.output_params.attributes.get("vtkFilename"),subdivision = 3)
else:
self.vtk= ngs.VTKOutput(ma = self.mesh,coefs = [self.gfu[0].components[0],self.gfu[0].components[1],self.gfu[0].components[2]],names=["gfu_x","gfu_y","gfu_z"],filename=self.output_params.attributes.get("vtkFilename"),subdivision = 3)
def ConstructForms(self,use_bdf=1):
if self.solver_type_name == "PPFEM":
if use_bdf == 1:
self.coeff_bdf = [1,0] # [1, -1] are the actual coeff, but we homogenize the dirichlet boundary condition by solving for dgfu := gfu[0] - gfu[1]
self.coeff_extrap = [1]
elif use_bdf == 2:
self.coeff_bdf =[1.5,-0.5,0.5] # [1.5, -2,0.5] same as bdf1
self.coeff_extrap = [2,-1]
else:
raise RuntimeError("Only BDF 1 and 2 supported!")
u = self.fes.TrialFunction()
v = self.fes.TestFunction()
## incremental form: "unknown = gfu - gfu_prev"
self.bf = ngs.BilinearForm(self.fes)
self.bf += self.coeff_bdf[0] * u * v * ngs.dx
self.bf += self.dt* ngs.InnerProduct(ngs.grad(u),ngs.grad(v))* ngs.dx
self.bf += -self.dt * ngs.InnerProduct(ngs.grad(self.gfu_extrap),ngs.grad(self.gfu_extrap)) * u * v * ngs.dx
self.lf = ngs.LinearForm(self.fes)
self.lf += -self.dt * ngs.InnerProduct(ngs.grad(self.gfu[1]),ngs.grad(v)) * ngs.dx
self.lf += self.dt*ngs.InnerProduct(ngs.grad(self.gfu_extrap),ngs.grad(self.gfu_extrap)) * self.gfu[1] * v *ngs.dx
for i in range(1,len(self.coeff_bdf)):
self.lf += -self.coeff_bdf[i]*self.gfu[i]*v*ngs.dx
elif self.solver_type_name == "TFEM":
if use_bdf == 1:
self.coeff_bdf = [1,-1]
self.coeff_extrap = [1]
elif use_bdf == 2:
self.coeff_bdf =[1.5,-2,0.5]
self.coeff_extrap = [2,-1]
else:
raise RuntimeError("Only BDF 1 and 2 supported!")
u_dot_x, u_dot_y, u_dot_z,lamb = self.fes.TrialFunction()
v_x,v_y,v_z,w = self.fes.TestFunction()
self.bf = ngs.BilinearForm(self.fes)
self.bf += (u_dot_x*v_x+u_dot_y*v_y+u_dot_z*v_z)*ngs.dx
self.bf += 1/self.coeff_bdf[0] * self.dt * (ngs.grad(u_dot_x)*ngs.grad(v_x)+ngs.grad(u_dot_y)*ngs.grad(v_y)+ngs.grad(u_dot_z)*ngs.grad(v_z))*ngs.dx
self.bf += (self.gfu_extrap.components[0]*v_x+self.gfu_extrap.components[1]*v_y+self.gfu_extrap.components[2]*v_z)*lamb*ngs.dx
self.bf += (self.gfu_extrap.components[0]*u_dot_x+self.gfu_extrap.components[1]*u_dot_y+self.gfu_extrap.components[2]*u_dot_z)*w*ngs.dx
self.lf = ngs.LinearForm(self.fes)
for i in range(1,len(self.coeff_bdf)):
self.lf += self.coeff_bdf[i]/self.coeff_bdf[0]*(ngs.grad(self.gfu[i].components[0])*ngs.grad(v_x)+ngs.grad(self.gfu[i].components[1])*ngs.grad(v_y)+ngs.grad(self.gfu[i].components[2])*ngs.grad(v_z))*ngs.dx
def Solve(self):
if self.solver_type_name != "CPFEM":
self.Solve2(self.time_init_params.attributes.get("tstart"),self.time_init_params.attributes.get("tend"),self.time_init_params.attributes.get("BDF"))
else:
self.Solve_CrossFEM(self.time_init_params.attributes.get("tstart"),self.time_init_params.attributes.get("tend"),self.space_init_params.attributes.get("CPFEM_max_iter"),self.space_init_params.attributes.get("CPFEM_eps"))
def Solve2(self,tstart,tend,use_bdf,one_step=False):
t = tstart
## first steps needed to do multi-step bdf
bdf_counter = use_bdf
counter = 0
while bdf_counter>1 :
bdf_counter = bdf_counter -1
## only compute one step
self.Solve2(t,tend,bdf_counter,True)
#print(" computed first time step with bdf ",bdf_counter)
t = t+self.dt.Get()
counter = counter+1
self.ConstructForms(use_bdf)
debug_mode = self.output_params.attributes.get("debug")
#if self.output_params.attributes.get("vtk") ==1 and use_bdf==1:
# self.vtk.Do(time=0)
while abs(tend-t)>1e-12 and tend>t:
t = t+self.dt.Get()
if debug_mode ==1:
print ("\rt=", t, end="")
#print ("\rt=", t)
counter = counter + 1
### update solution of previous steps
for i in range(len(self.gfu)-1,0,-1):
if self.solver_type_name != "TFEM":
self.gfu[i].vec.data = self.gfu[i-1].vec
else:
self.gfu[i].components[0].vec.data = self.gfu[i-1].components[0].vec
self.gfu[i].components[1].vec.data = self.gfu[i-1].components[1].vec
self.gfu[i].components[2].vec.data = self.gfu[i-1].components[2].vec
### update extrapolation
self.gfu_extrap.vec[:]=0
for i in range(len(self.coeff_extrap)):
if self.solver_type_name != "TFEM":
self.gfu_extrap.vec.data +=self.coeff_extrap[i]*self.gfu[i+1].vec
else:
self.gfu_extrap.components[0].vec.data +=self.coeff_extrap[i]*self.gfu[i+1].components[0].vec
self.gfu_extrap.components[1].vec.data +=self.coeff_extrap[i]*self.gfu[i+1].components[1].vec
self.gfu_extrap.components[2].vec.data +=self.coeff_extrap[i]*self.gfu[i+1].components[2].vec
if self.solver_type_name == "TFEM":
normalized_gfu = ngs.Normalize(ngs.CoefficientFunction((self.gfu_extrap.components[0],self.gfu_extrap.components[1],self.gfu_extrap.components[2] )))
self.gfu_extrap.components[0].Interpolate(normalized_gfu[0])
self.gfu_extrap.components[1].Interpolate(normalized_gfu[1])
self.gfu_extrap.components[2].Interpolate(normalized_gfu[2])
### Solve it
with ngs.TaskManager():
self.bf.Assemble()
self.lf.Assemble()
### TODO: put this into some postprocessing
if self.solver_type_name == "PPFEM":
self.gfu[0].vec.data += self.bf.mat.Inverse(freedofs=self.fes.FreeDofs())*self.lf.vec
normalized_gfu = ngs.Normalize(self.gfu[0])
self.gfu[0].Interpolate(normalized_gfu)
elif self.solver_type_name == "TFEM":
self.gfu_dot.vec.data = self.bf.mat.Inverse(freedofs=self.fes.FreeDofs())*self.lf.vec
self.gfu[0].components[0].vec.data = 1/self.coeff_bdf[0]*self.dt.Get() * self.gfu_dot.components[0].vec
self.gfu[0].components[1].vec.data = 1/self.coeff_bdf[0]*self.dt.Get() * self.gfu_dot.components[1].vec
self.gfu[0].components[2].vec.data = 1/self.coeff_bdf[0]*self.dt.Get() * self.gfu_dot.components[2].vec
for i in range(1,len(self.coeff_bdf)):
self.gfu[0].components[0].vec.data += -self.coeff_bdf[i]/self.coeff_bdf[0]*self.gfu[i].components[0].vec
self.gfu[0].components[1].vec.data += -self.coeff_bdf[i]/self.coeff_bdf[0]*self.gfu[i].components[1].vec
self.gfu[0].components[2].vec.data += -self.coeff_bdf[i]/self.coeff_bdf[0]*self.gfu[i].components[2].vec
### output
if self.output_params.attributes.get("vtk") ==1 and counter%self.output_params.attributes.get("vtkInterval") ==0:
self.vtk.Do(time=t)
#print("Energy = ",ngs.Integrate(ngs.InnerProduct(ngs.grad(self.gfu[0]),ngs.grad(self.gfu[0])),self.mesh))
if one_step:
#print("only did one step")
break
def ChangeTimeStep(self,tstep):
self.dt.Set(tstep)
def Solve_CrossFEM(self,tstart,tend,max_fp_iter,epsi):
debug_mode = self.output_params.attributes.get("debug")
v = self.fes.TestFunction()
w = self.fes.TrialFunction()
gfu_w = ngs.GridFunction(self.fes)
gfu_u = ngs.GridFunction(self.fes)
gfu_e = ngs.GridFunction(self.fes)
gfu_w_disc_laplacian = ngs.GridFunction(self.fes)
gfu_e_disc_laplacian = ngs.GridFunction(self.fes)
gfu_w.Set(self.initial_cond)
gfu_u.Set(self.initial_cond)
## Solution in prev step
gfu_w_prev = ngs.GridFunction(self.fes)
gfu_u_prev = ngs.GridFunction(self.fes)
## For reduced L2 product
massLumping = ngs.IntegrationRule( points = [(0,0), (1,0), (0,1)],
weights = [1/6, 1/6, 1/6] )
## Fixed point method formulation in incremental form:
a_2D = ngs.BilinearForm(self.fes)
a_2D+= 0.5 * self.dt * ngs.Cross(w,ngs.Cross(gfu_w_prev,gfu_w_disc_laplacian)) * v * ngs.dx(intrules={ngs.TRIG:massLumping})
a_2D+= w * v * ngs.dx(intrules={ngs.TRIG:massLumping})
f_2D = ngs.LinearForm(self.fes)
f_2D+= -0.5 * self.dt * ngs.Cross(gfu_w_prev, ngs.Cross(gfu_w_prev,gfu_w_disc_laplacian)) * v * ngs.dx(intrules={ngs.TRIG:massLumping})
f_2D+= -gfu_w_prev * v * ngs.dx(intrules={ngs.TRIG:massLumping})
f_2D+= gfu_u_prev * v * ngs.dx(intrules={ngs.TRIG:massLumping})
## Forms for discrete Laplacian
a_discrete_laplacian = ngs.BilinearForm(self.fes)
a_discrete_laplacian += - w * v * ngs.dx(intrules={ngs.TRIG:massLumping})
a_discrete_laplacian.Assemble()
A_discrete_laplacian_inv = a_discrete_laplacian.mat.Inverse(freedofs=self.fes.FreeDofs())
f_discrete_laplacian = ngs.LinearForm(self.fes)
f_discrete_laplacian += ngs.InnerProduct(ngs.grad(gfu_w_prev),ngs.grad(v)) * ngs.dx
f_e_discrete_laplacian = ngs.LinearForm(self.fes)
f_e_discrete_laplacian += ngs.InnerProduct(ngs.grad(gfu_e),ngs.grad(v)) * ngs.dx
t = tstart
count_fp_iter = 0
while abs(tend-t) > 1e-12:
t = t+self.dt.Get()
if debug_mode==1:
print ("t=", t)
gfu_u_prev.vec.data = gfu_u.vec
gfu_w.vec.data = gfu_u.vec
while count_fp_iter < max_fp_iter:
with ngs.TaskManager():
gfu_w_prev.vec.data = gfu_w.vec
f_discrete_laplacian.Assemble()
gfu_w_disc_laplacian.vec.data = A_discrete_laplacian_inv * f_discrete_laplacian.vec
a_2D.Assemble()
f_2D.Assemble()
inv_a = a_2D.mat.Inverse(freedofs=self.fes.FreeDofs())
gfu_e.vec.data = inv_a * f_2D.vec
gfu_w.vec.data += gfu_e.vec
f_e_discrete_laplacian.Assemble()
gfu_e_disc_laplacian.vec.data = A_discrete_laplacian_inv * f_e_discrete_laplacian.vec
r = ngs.sqrt(ngs.Integrate(ngs.InnerProduct( ngs.Cross(gfu_w,gfu_e_disc_laplacian)+ngs.Cross(gfu_e,gfu_w_disc_laplacian),ngs.Cross(gfu_w,gfu_e_disc_laplacian)+ngs.Cross(gfu_e,gfu_w_disc_laplacian) ),self.mesh) )
count_fp_iter = count_fp_iter +1
if debug_mode==1:
print(" iteration step: ", count_fp_iter, "with r=",r)
if r < epsi :
gfu_u.vec.data = 2 * gfu_w.vec - gfu_u_prev.vec
if debug_mode==1:
print(" number of iterations needed: ",count_fp_iter)
count_fp_iter = 0
break
if (count_fp_iter == max_fp_iter ):
raise RuntimeError("not enough fix point iteration steps")
### set to class variable
self.gfu[1].vec.data = gfu_u_prev.vec
self.gfu[0].vec.data = gfu_u.vec
def ComputeH1Error2D(gfu_1,gfu_2,mesh_ref):
with ngs.TaskManager():
l2error = ngs.sqrt(ngs.Integrate((gfu_1-gfu_2)**2*ngs.dx,mesh_ref))
graderror = ngs.sqrt(ngs.Integrate(ngs.InnerProduct(ngs.grad(gfu_1)-ngs.grad(gfu_2),ngs.grad(gfu_1)-ngs.grad(gfu_2))*ngs.dx,mesh_ref))
h1error =ngs.sqrt(graderror**2+l2error**2)
#print("\nL2 error = ",l2error)
#print("\nH1 error = ",h1error)
return l2error, h1error
def ComputeH1Error2DForCF(cf_1,cf_2,grad_cf_1,grad_cf_2,mesh):
with ngs.TaskManager():
l2error = ngs.sqrt(ngs.Integrate((cf_1-cf_2)**2*ngs.dx,mesh))
graderror = ngs.sqrt(ngs.Integrate(ngs.InnerProduct(grad_cf_1-grad_cf_2,grad_cf_1-grad_cf_2)*ngs.dx,mesh))
h1error =ngs.sqrt(graderror**2+l2error**2)
#print("\nL2 error = ",l2error)
#print("\nH1 error = ",h1error)
return l2error, h1error
from netgen.geom2d import SplineGeometry
import ngsolve as ngs
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve
import numpy as np
## transforms RSHMHF to HMHF on unit disk in 2D
## h_ref : mesh size for the mesh_2D of gfu_2D (comparable to mesh size of RSHMHF reference solution gfu_1D
## initial_cond: 3-vectorial 2-dim initial condition of HMHF2D
## FirstOrder = fes first order = nodal basis
def SphericalTransform_FirstOrder(gfu_1D,mesh_1D,h_ref,initial_cond):
geo = SplineGeometry()
geo.AddCircle((0,0),1,bc="unit_circle")
mesh_2D = ngs.Mesh(geo.GenerateMesh(maxh=h_ref))
fes_2D = ngs.H1(mesh_2D,order=1,dirichlet=geo.GetBCName(1))
gfu_2D_1 = ngs.GridFunction(fes_2D)
gfu_2D_1.Set(initial_cond[0])
gfu_2D_2 = ngs.GridFunction(fes_2D)
gfu_2D_2.Set(initial_cond[1])
gfu_2D_3 = ngs.GridFunction(fes_2D)
gfu_2D_3.Set(initial_cond[2])
j=0
for p in mesh_2D.ngmesh.Points():
r_p = ngs.sqrt(p[0]**2 + p[1]**2);
if abs(r_p -1) < 1e-12: # dont change anything at the boundary
j = j+1
continue
gfu_2D_1.vec[j] = p[0]/r_p * ngs.sin(gfu_1D(mesh_1D(r_p,0.0)))
gfu_2D_2.vec[j] = p[1]/r_p * ngs.sin(gfu_1D(mesh_1D(r_p,0.0)))
gfu_2D_3.vec[j] = ngs.cos(gfu_1D(mesh_1D(r_p,0.0)))
j = j+1
gfu_2D = ngs.CoefficientFunction((gfu_2D_1,gfu_2D_2,gfu_2D_3))
## Jacobian of gfu_2D, compareable to ngs.grad(gfu_2D)
gradgfu_2D = ngs.CoefficientFunction((ngs.grad(gfu_2D_1)[0],ngs.grad(gfu_2D_2)[0],ngs.grad(gfu_2D_3)[0],ngs.grad(gfu_2D_1)[1],ngs.grad(gfu_2D_2)[1],ngs.grad(gfu_2D_3)[1]), dims=(2,3))
## Gradient of gfu_2D, compareable to ngs.Grad(gfu_2D)
#gradgfu_2D = ngs.CoefficientFunction((ngs.grad(gfu_2D_1)[0],ngs.grad(gfu_2D_1)[1],ngs.grad(gfu_2D_2)[0],ngs.grad(gfu_2D_2)[1],ngs.grad(gfu_2D_3)[0],ngs.grad(gfu_2D_3)[1]), dims=(3,2))
return mesh_2D, fes_2D, gfu_2D, gradgfu_2D
## transforms RSHMHF (short "1D") to HMHF on unit disk in 2D
## h_ref : mesh size for the mesh_2D of gfu_2D (comparable to mesh size of RSHMHF reference solution gfu_1D
## initial_cond: to set for the boundary
## fes_order: should be same order as fes of gfu_1D
def SphericalTransform(gfu_1D,mesh_1D,h_ref,initial_cond, fes_order):
if fes_order == 1:
mesh_2D, fes_2D, gfu_2D, gradgfu_2D = SphericalTransform_FirstOrder(gfu_1D,mesh_1D,h_ref,initial_cond)
return mesh_2D, fes_2D, gfu_2D, gradgfu_2D
geo = SplineGeometry()
geo.AddCircle((0,0),1,bc="unit_circle")
mesh_2D = ngs.Mesh(geo.GenerateMesh(maxh=h_ref))
#mesh_2D.Curve(2)
mesh_2D.ngmesh.SecondOrder()
#fes_2D = ngs.H1(mesh_2D,order=fes_order,dirichlet=geo.GetBCName(1))
fes_2D = ngs.H1(mesh_2D,order=fes_order)
gfu_2D_1 = ngs.GridFunction(fes_2D)
gfu_2D_1.Set(initial_cond[0])
gfu_2D_2 = ngs.GridFunction(fes_2D)
gfu_2D_2.Set(initial_cond[1])
gfu_2D_3 = ngs.GridFunction(fes_2D)
gfu_2D_3.Set(initial_cond[2])
gfu_fes = ngs.GridFunction(fes_2D)
fes_ndof = fes_2D.ndof
A_1 = lil_matrix((fes_ndof,fes_ndof),dtype = np.float32)
f_1 = np.zeros([fes_ndof,1])
f_2 = np.zeros([fes_ndof,1])
f_3 = np.zeros([fes_ndof,1])
x_ng_pnt = np.ndarray(shape=fes_2D.ndof)
y_ng_pnt = np.ndarray(shape=fes_2D.ndof)
j=0
for p in mesh_2D.ngmesh.Points():
x_ng_pnt[j] = p[0]
y_ng_pnt[j] = p[1]
r_p =ngs.sqrt(p[0]**2+p[1]**2)
f_1[j,0] = p[0]/r_p * ngs.sin(gfu_1D(mesh_1D(r_p,0.0,0.0)))
f_2[j,0] = p[1]/r_p * ngs.sin(gfu_1D(mesh_1D(r_p,0.0,0.0)))
f_3[j,0] = ngs.cos(gfu_1D(mesh_1D(r_p,0.0,0.0)))
j+=1
pnt = mesh_2D(x_ng_pnt,y_ng_pnt)
for i in range(fes_ndof):
gfu_fes.vec[:] = 0
gfu_fes.vec[i] = 1
A_1[i] = lil_matrix(gfu_fes(pnt)).transpose()
A = A_1.transpose()
eps = 1e-12
for i in range(A.shape[0]):
row_data = A.data[i]
row_cols = A.rows[i]
A.data[i] = [val for val in row_data if abs(val) >= eps]
A.rows[i] = [col for val, col in zip(row_data, row_cols) if abs(val) >= eps]
A = A.tocsr()
gfu_2D_1.vec.data= spsolve(A,f_1)
gfu_2D_2.vec.data= spsolve(A,f_2)
gfu_2D_3.vec.data= spsolve(A,f_3)
#print("solving system successfully")
gfu_2D = ngs.CoefficientFunction((gfu_2D_1,gfu_2D_2,gfu_2D_3))
gradgfu_2D = ngs.CoefficientFunction((ngs.grad(gfu_2D_1)[0],ngs.grad(gfu_2D_2)[0],ngs.grad(gfu_2D_3)[0],ngs.grad(gfu_2D_1)[1],ngs.grad(gfu_2D_2)[1],ngs.grad(gfu_2D_3)[1]), dims=(2,3))
return mesh_2D, fes_2D, gfu_2D, gradgfu_2D
......@@ -2,9 +2,28 @@
## Installation
[matlab](https://www.mathworks.com/products/matlab.html) are needed to run the code.
[matlab](https://www.mathworks.com/products/matlab.html) (at least R2022b) and [ngsolve](https://www.ngsolve.org) (latest version) are needed to run the code.
## Usage
- In `RSHMHF_FD`: Finite difference discretizations with [matlab](https://www.mathworks.com/products/matlab.html) for radially symmetric HMHF (RSHMHF) used in [arXiv:2503.03525](https://arxiv.org/abs/2503.03525)
- `RSHMHF_FD`: BDF1/BDF2 + Finite difference discretizations with [matlab](https://www.mathworks.com/products/matlab.html) for radially symmetric HMHF (RSHMHF)
- "Discretization error analysis for a radially symmetric harmonic map heat flow problem" - Nguyen, Reusken (2025) [arXiv](https://arxiv.org/abs/2503.03525)
- `RSHMHF_FE`: BDF1/BDF2 + P1/P2 Finite element discretizations with [ngsolve](https://www.ngsolve.org) for RSHMHF
- `convex_lin`: "Error analysis for a finite element discretization of a radially symmetric harmonic map heat flow problem" - Nguyen, Reusken (2025)
- `max_lin`: "A comparative study of finite element methods for a class of harmonic map heat flow problems" - Nguyen, Reusken (2025)
- `HMHF_FE`: BDF1/BDF2 + P1/P2 Finite element discretizations wtih [ngsolve](https://www.ngsolve.org) for HMHF
- `Dim2d`:
- `PPFEM`: p.139 of "Computational Micromagnetism" - Andreas Prohl (2001) [DOI](https://doi.org/10.1007/978-3-663-09498-2)
- `TFEM`: "Higher-order linearly implicit full discretization of the Landau-Lifshitz-Gilbert equation" - Akrivis, Feischl, Kovács, Lubich (2021) [DOI](https://dx.doi.org/10.1090/mcom/3597)
- `CPFEM`: "Constraint preserving implicit finite element discretization of harmonic map flow into spheres" - Bartels, Prohl (2007) [DOI](https://dx.doi.org/10.1090/S0025-5718-07-02026-1)
## Tests
Run the following in command line to test mesh size and time step convergence for P1 and BDF1 of RSHMHF FE
```console
foo@bar:~$ python3 -m Tests.conv_RSHMHF_FE
Run the following in command line to test mesh size and time step convergence for P1 and BDF1 of PPFEM and TFEM of HMHF FE on unit disk in dimension 2
```console
foo@bar:~$ python3 -m Tests.conv_HMHF2D_FE
```
import ngsolve as ngs
# p0_fes := order 0 L2-space of mesh
# gfu is current solution and lives in solution space H1 order 1 on mesh
# return: piecewise gradient of gfu on mesh in p0_fes_space
def ComputeP0Gradient(mesh,p0_fes,gfu):
assert p0_fes.ndof+1 == len(gfu.vec)
p0_grad = ngs.GridFunction(p0_fes)
mesh_size_element_wise = ngs.Integrate(ngs.CoefficientFunction(1)*ngs.dx,mesh,element_wise=True)
for i in range(0, len(p0_grad.vec)):
p0_grad.vec[i] = (gfu.vec[i+1] - gfu.vec[i])/mesh_size_element_wise[i]
return p0_grad
# p1_fes := order 1 H1 space on mesh
# p0_grad : piecewise constant gradient in order 0 L2 on mesh
# return : ZZ post processed gradient in p1_fes on mesh
def ComputeZZGradient(mesh,p1_fes,p0_grad):
assert p1_fes.ndof == len(p0_grad.vec)+1
zz_grad_gfu = ngs.GridFunction(p1_fes)
#print("\n",p1_fes.ndof)
#print(len(zz_grad_gfu.vec))
for i in range(0,p1_fes.ndof):
if i==0:
zz_grad_gfu.vec[0] = p0_grad.vec[0]
continue
if i==p1_fes.ndof-1:
zz_grad_gfu.vec[p1_fes.ndof-1] = p0_grad.vec[p1_fes.ndof-2]
continue
zz_grad_gfu.vec[i] = 0.5*(p0_grad.vec[i] + p0_grad.vec[i-1])
return zz_grad_gfu
def ComputeZZIndicator(mesh,p1_fes,gfu):
W = ngs.L2(mesh,order = 0)
p0_grad_gfu = ComputeP0Gradient(mesh,W,gfu)
zz_grad_gfu = ComputeZZGradient(mesh,p1_fes,p0_grad_gfu)
err_elwise = ngs.Integrate((zz_grad_gfu-p0_grad_gfu)**2*ngs.dx,mesh,element_wise = True)
err_total = ngs.Integrate((zz_grad_gfu-p0_grad_gfu)**2*ngs.dx,mesh)
return err_elwise, err_total
from .parameters_RSHMHF import *
from .meshing_1D import *
from .ZZ_indicator_1D import *
from .marking_strategy import *
# gfu_list := [gfu at current time j, gfu at j-1, gfu at j-2,... ] in p1_fes on mesh
# p1_fes := order 1 H1 space on mesh
# tstep := current time step
# tcurrent := current time
# tend := end time
# number_elements := number of elements in mesh
# time_adaptivity_type: "TimeDerivativeGradient"
# space_adaptive_type: "ZZRecovery", "Energy"
# marking_strategy: "Doerfler", "equidistributed"
def SpaceTimeAdaptivity(mesh,gfu_list,gfu_extrap,p1_fes,tstep,tcurrent,tend,number_elements,adapt_params,counter_tmin,stop_space_adapt = False,stop_time_adapt=False,debug = 0):
use_adaptivity = adapt_params.attributes.get("use_adaptivity")
tstep_min = adapt_params.attributes.get("tstep_min")
tstep_max = adapt_params.attributes.get("tstep_max")
alpha_plus = adapt_params.attributes.get("alpha_plus")
alpha_minus = adapt_params.attributes.get("alpha_minus")
eps_err = adapt_params.attributes.get("eps_err")
eps_mesh_size = adapt_params.attributes.get("eps_mesh_size")
theta = adapt_params.attributes.get("theta")
space_adaptive_option = adapt_params.attributes.get("space_adaptive_option")
time_adaptive_option = adapt_params.attributes.get("time_adaptive_option")
marking_strategy = adapt_params.attributes.get("marking_strategy")
new_tstep = tstep
new_tcurrent = tcurrent
new_mesh = mesh
new_gfu_list = gfu_list
new_number_elements = number_elements
new_gfu_extrap = gfu_extrap
dirichlet_bnd = mesh.GetBoundaries()[0]+ "|" + mesh.GetBoundaries()[1]
#print(dirichlet_bnd)
new_p1_fes = ngs.H1(mesh,order=1,dirichlet = dirichlet_bnd)
refined_something = False
if not use_adaptivity:
return new_mesh,new_p1_fes,new_gfu_list, new_gfu_extrap, new_number_elements, new_tstep, new_tcurrent,refined_something
time_indicator = 0.0
#time_indicator_elwise = []
space_indicator = 0.0
space_indicator_elwise = []
W = ngs.L2(mesh,order = 0)
if time_adaptive_option == "TimeDerivativeGradient" and not stop_time_adapt:
p0_grad_gfu = ComputeP0Gradient(mesh,W,gfu_list[0])
p0_grad_gfu_prev = ComputeP0Gradient(mesh,W,gfu_list[1])
time_indicator = ngs.Integrate((p0_grad_gfu-p0_grad_gfu_prev)**2 * ngs.dx, mesh)
#time_indicator_elwise = ngs.Integrate((p0_grad_gfu-p0_grad_gfu_prev)**2 * ngs.dx, mesh, element_wise=True)
elif time_adaptive_option == "TimeDerivativeGradient" and stop_time_adapt:
time_indicator = 0.0
elif time_adaptive_option != "TimeDerivativeGradient":
raise RuntimeError("Time indicator of type ",time_adaptive_option, " not supported")
if space_adaptive_option == "ZZRecovery" and not stop_space_adapt:
p0_grad_gfu = ComputeP0Gradient(mesh,W,gfu_list[0])
zz_grad_gfu = ComputeZZGradient(mesh,p1_fes,p0_grad_gfu)
space_indicator_elwise = ngs.Integrate((zz_grad_gfu-p0_grad_gfu)**2*ngs.dx,mesh,element_wise = True)
space_indicator = ngs.Integrate((zz_grad_gfu-p0_grad_gfu)**2*ngs.dx,mesh)
#space_indicator_elwise,space_indicator = ComputeZZIndicator(mesh,p1_fes,gfu_list[0])
elif space_adaptive_option == "ZZRecovery" and stop_space_adapt:
space_indicator = 0
elif space_adaptive_option != "ZZRecovery":
raise RuntimeError("Space indicator of type ",space_adaptive_option, " not supported")
#print("\n eps_err = ",eps_err)
if space_indicator + time_indicator < eps_err:
refined_something = False
# incrase time step if possible
if not stop_time_adapt:
new_tstep = min([tstep_max,alpha_plus*tstep,tend-tcurrent])
counter_tmin = 0
else:
# reject current solution
gfu_list[0].vec.data = gfu_list[1].vec
refined_something = True
new_tcurrent = tcurrent-tstep
ref_flags = [False]*number_elements
# if minimal time step is used too many times without decreasing error, then refine instead
if counter_tmin > 3:
time_indicator = 0
counter_tmin = 0
# check which error dominates
if time_indicator > space_indicator:
# reduce time step if possible
new_tstep = max([tstep_min,alpha_minus*tstep])
if abs(new_tstep -tstep_min)<1e-12:
counter_tmin = counter_tmin +1
else:
index_to_be_refined = SortErrListAndReturnElementIndex(number_elements,space_indicator_elwise,theta,space_indicator)
for index in index_to_be_refined:
ref_flags[index] = True
new_mesh, new_number_elements = Refine1DMesh(mesh,ref_flags)
new_p1_fes = ngs.H1(new_mesh,order=1,dirichlet=dirichlet_bnd)
new_gfu_list = []
for i in range(0,len(gfu_list)):
new_gfu_list.append(ngs.GridFunction(new_p1_fes))
new_gfu_list[i].Interpolate(gfu_list[i])
new_gfu_extrap = ngs.GridFunction(new_p1_fes)
new_gfu_extrap.Interpolate(gfu_extrap)
counter_tmin = 0
if debug:
print("\n time indicator = ",time_indicator)
print("\n space indicator = ",space_indicator)
return new_mesh,new_p1_fes,new_gfu_list, new_gfu_extrap, new_number_elements, new_tstep, new_tcurrent,refined_something, counter_tmin
import ngsolve as ngs
import numpy as np
# Doerfler strategy
def SortErrListAndReturnElementIndex(num_elem,err_elem,theta,err_total):
assert len(err_elem) == num_elem
elerr_index = [i for i in range(0,num_elem)]
elerr_sorted = list( sorted(err_elem,reverse=True))
elerr_index_sorted =list( np.argsort(err_elem)[::-1][:num_elem])
index_refine =[]
elerr_refine = []
#print("total error",err_total)
#print("err list")
#for i,j in zip(elerr,elerr_index):
# print ("%.12f"%i,j)
#print("")
#print("sorted err list")
#for i,j in zip(elerr_sorted,elerr_index_sorted):
# print ("%.12f"%i,j)
#print("")
while sum(elerr_refine) < theta * err_total and len(elerr_index_sorted) > 0:
elerr_refine.append(elerr_sorted.pop(0))
index_refine.append(elerr_index_sorted.pop(0))
#for i,j in zip(elerr_refine,index_refine):
# print ("%.12f"%i,j)
#print("sum of elerr_refine so far ", "%.12f"%sum(elerr_refine), "; is it less than ", theta*err_total)
return index_refine
import netgen.meshing as ng
import ngsolve as ngs
def Print1DMesh(mesh):
current_discr_points =[p[0] for p in mesh.ngmesh.Points()]
print(current_discr_points)
def Create1DMesh(number_elements,interval_start,interval_end,region_name_left = "left",region_name_right="right",region_name_inner = "material"):
m = ng.Mesh(dim=1)
a = interval_start
b = interval_end
interval_length = b-a
assert(interval_length > 0)
assert(number_elements > 0)
pnums = []
for i in range(0, number_elements+1):
pnums.append (m.Add (ng.MeshPoint (ng.Pnt(interval_start + interval_length * i/number_elements, 0, 0))))
idx = m.AddRegion(region_name_inner, dim=1)
for i in range(0,number_elements):
m.Add (ng.Element1D ([pnums[i],pnums[i+1]], index=idx))
idx_left = m.AddRegion(region_name_left, dim=0)
idx_right = m.AddRegion(region_name_right, dim=0)
m.Add (ng.Element0D (pnums[0], index=idx_left))
m.Add (ng.Element0D (pnums[number_elements], index=idx_right))
mesh = ngs.Mesh(m)
return mesh
def Create1DMeshFromPoints(discr_points,region_name_left = "left",region_name_right="right",region_name_inner = "material"):
pnums = []
m = ng.Mesh(dim=1)
for i in range(0, len(discr_points)):
pnums.append (m.Add (ng.MeshPoint (ng.Pnt(discr_points[i], 0, 0))))
idx = m.AddRegion(region_name_inner, dim=1)
for i in range(0,len(discr_points)-1):
m.Add (ng.Element1D ([pnums[i],pnums[i+1]], index=idx))
idx_left = m.AddRegion(region_name_left, dim=0)
idx_right = m.AddRegion(region_name_right, dim=0)
m.Add (ng.Element0D (pnums[0], index=idx_left))
m.Add (ng.Element0D (pnums[len(discr_points)-1], index=idx_right))
mesh = ngs.Mesh(m)
return mesh
# returns refined mesh and number of elements of the refined mesh
# number_elements + 1 = number_discretization in 1D
def Refine1DMesh(mesh,refinement_flags):
current_discr_points =[p[0] for p in mesh.ngmesh.Points()]
number_discretization_pts = len(current_discr_points)
assert number_discretization_pts == len(current_discr_points)
assert number_discretization_pts == len(refinement_flags)+1
assert len(mesh.GetMaterials()) == 1
assert len(mesh.GetBoundaries()) == 2
ref_discr_points = current_discr_points
nmb_added_new_points = 0
for i in range(0,len(refinement_flags)):
if refinement_flags[i]:
ref_discr_points.insert(i+1+nmb_added_new_points, 0.5*(current_discr_points[i+nmb_added_new_points]+current_discr_points[i+1+nmb_added_new_points]))
nmb_added_new_points = nmb_added_new_points +1
#print(ref_discr_points)
number_elements_refined_mesh = len(ref_discr_points)-1
mesh_new = Create1DMeshFromPoints(ref_discr_points, mesh.GetBoundaries()[0],mesh.GetBoundaries()[1],mesh.GetMaterials()[0])
return mesh_new, number_discretization_pts + nmb_added_new_points-1
import math
import sys
import fileinput
class ParametersRSHMHF:
def __init__(self):
self.attributes = dict([])
def ReadParametersFromFile(self,filename=""):
if len(sys.argv) < 2 and filename == "":
raise RuntimeError("No input file!")
if filename == "":
filename = sys.argv[1]
for line in fileinput.input(files=filename):
#print("line ", fileinput.lineno())
if line[0]=="#":
continue
arr_string = line.split()
if len(arr_string) == 0:
#print("Nothing in here")
continue
if len(arr_string) < 2:
raise RuntimeError('In line {} of the file "{}" an argument is missing!'.format(fileinput.lineno(),fileinput.filename()))
if len(arr_string)==2:
if arr_string[0] in self.attributes:
self.attributes[arr_string[0]]= arr_string[1]
def PrintParameters(self):
attributes_list = list(self.attributes)
for iterator in attributes_list:
print(iterator, " = ", self.attributes.get(iterator))
class TimeInitialParametersRSHMHF(ParametersRSHMHF):
def __init__(self):
ParametersRSHMHF.__init__(self)
self.attributes.update({"tstart":0.0, "tend":1.0,\
"tstep":0.1,"BDF":1,"solver_type":"convex_lin"})
def ReadParametersFromFile(self,filename=""):
ParametersRSHMHF.ReadParametersFromFile(self,filename)
self.attributes["tstart"] = float(self.attributes.get("tstart"))
self.attributes["tend"] = float(self.attributes.get("tend"))
self.attributes["tstep"] = float(self.attributes.get("tstep"))
self.attributes["BDF"] = int(self.attributes.get("BDF"))
self.attributes["solver_type"] = self.attributes.get("solver_type")
def PrintParameters(self):
ParametersRSHMHF.PrintParameters(self)
class SpaceInitialParametersRSHMHF(ParametersRSHMHF):
def __init__(self):
ParametersRSHMHF.__init__(self)
self.attributes.update({"interval_start":0.0, "interval_end":1.0, \
"number_elements":10,"refinement_level":0,\
"refinement_until_this_point":1.0, \
"Dirichlet":1, "fes_order":1})
def ReadParametersFromFile(self,filename=""):
ParametersRSHMHF.ReadParametersFromFile(self,filename)
self.attributes["interval_start"] = float(self.attributes.get("interval_start"))
self.attributes["interval_end"] = float(self.attributes.get("interval_end"))
self.attributes["number_elements"] = int(self.attributes.get("number_elements"))
self.attributes["refinement_level"] = int(self.attributes.get("refinement_level"))
self.attributes["refinement_until_this_point"] \
= float(self.attributes.get("refinement_until_this_point"))
self.attributes["Dirichlet"] = int(self.attributes.get("Dirichlet"))
self.attributes["fes_order"] = int(self.attributes.get("fes_order"))
def PrintParameters(self):
ParametersRSHMHF.PrintParameters(self)
### ToDo: Add output options
class OutputParametersRSHMHF(ParametersRSHMHF):
def __init__(self):
ParametersRSHMHF.__init__(self)
self.attributes.update({"debug": 1, "save_mesh":"mesh.dat", \
"save_gfu":"gfu.dat", "save_energy":"energy.dat", \
"save_time":"time.dat", "save_tstep":"tstep.dat", \
"save_ndof":"ndof.dat", \
"save_temporal_error_indicator":"eta_tau.dat", \
"save_spacial_error_indicator":"eta_h.dat", \
"save":"0"
})
def ReadParametersFromFile(self,filename=""):
ParametersRSHMHF.ReadParametersFromFile(self,filename)
self.attributes["debug"] = int(self.attributes.get("debug"))
self.attributes["save"] = int(self.attributes.get("save"))
def PrintParameters(self):
ParametersRSHMHF.PrintParameters(self)
class AdaptivityParametersRSHMHF(ParametersRSHMHF):
def __init__(self):
ParametersRSHMHF.__init__(self)
self.attributes.update({"use_adaptivity":0,"tstep_min":1e-8, \
"tstep_max":1e-1, "alpha_minus":0.1, \
"alpha_plus":2.0, "eps_err":5e-1, \
"eps_mesh_size":1e-12, "theta":0.5, \
"space_adaptive_option":"", \
"time_adaptive_option":"", \
"marking_strategy":"Doerfler", \
"blowup_rate":0
})
def ReadParametersFromFile(self,filename=""):
ParametersRSHMHF.ReadParametersFromFile(self,filename)
self.attributes["use_adaptivity"] = int(self.attributes.get("use_adaptivity"))
self.attributes["tstep_min"] = float(self.attributes.get("tstep_min"))
self.attributes["tstep_max"] = float(self.attributes.get("tstep_max"))
self.attributes["alpha_minus"] = float(self.attributes.get("alpha_minus"))
self.attributes["alpha_plus"] = float(self.attributes.get("alpha_plus"))
self.attributes["eps_err"] = float(self.attributes.get("eps_err"))
self.attributes["eps_mesh_size"] = float(self.attributes.get("eps_mesh_size"))
self.attributes["theta"] = float(self.attributes.get("theta"))
self.attributes["blowup_rate"] = float(self.attributes.get("blowup_rate"))
def PrintParameters(self):
ParametersRSHMHF.PrintParameters(self)
def ReadAndPrintAllParametersFromFile(filename="",debug = 1):
tparams = TimeInitialParametersRSHMHF()
tparams.ReadParametersFromFile(filename)
sparams = SpaceInitialParametersRSHMHF()
sparams.ReadParametersFromFile(filename)
oparams = OutputParametersRSHMHF()
oparams.ReadParametersFromFile(filename)
aparams = AdaptivityParametersRSHMHF()
aparams.ReadParametersFromFile(filename)
if debug:
tparams.PrintParameters()
sparams.PrintParameters()
oparams.PrintParameters()
aparams.PrintParameters()
return tparams,sparams,oparams,aparams
from ngsolve.solvers import *
from .adaptivity_1D import *
import sys
import numpy as np
class Blowuprate():
def __init__(self):
self.dgfudx = []
self.time = []
self.value = 0.0
class SolvingRSHMHF:
def __init__(self, mesh, fes, initial_cond, params):
self.mesh = mesh
self.fes = fes
self.initial_cond = initial_cond
self.time_init_params = params.get("time")
### ToDo: add some output options
self.output_params = params.get("output")
self.adaptivity_params = params.get("adaptivity")
self.gfu = [] # store numerical solution at current timestep j [0], previous j-1 [1], j-2 [2] and so on depending on how many steps needed in bdf scheme
for i in range(self.time_init_params.attributes.get("BDF")+1):
self.gfu.append(ngs.GridFunction(self.fes))
self.gfu[i].Set(self.initial_cond)
#self.gfu[i].vec[0] = 0.0
#print("number of gfus = ",len(self.gfu), " for bdf ",self.time_init_params.attributes.get("BDF"))
self.gfu_extrap = ngs.GridFunction(self.fes)
tstep =self.time_init_params.attributes.get("tstep")
self.dt = ngs.Parameter(tstep)
self.bf = ngs.BilinearForm(self.fes)
self.lf = ngs.LinearForm(self.fes)
self.blowuprate = Blowuprate()
def ConstructForms(self,use_bdf=1):
if use_bdf == 1:
coeff_bdf = [1,0] # [1, -1] are the actual coeff, but we homogenize the dirichlet boundary condition by solving for dgfu := gfu[0] - gfu[1]
coeff_extrap = [1]
elif use_bdf == 2:
coeff_bdf =[1.5,-0.5,0.5] # [1.5, -2,0.5] same as bdf1
coeff_extrap = [2,-1]
else:
raise RuntimeError("Only BDF 1 and 2 supported!")
self.gfu_extrap.vec[:]=0
for i in range(len(coeff_extrap)):
self.gfu_extrap.vec.data +=coeff_extrap[i]*self.gfu[i+1].vec
u = self.fes.TrialFunction()
v = self.fes.TestFunction()
self.bf = ngs.BilinearForm(self.fes)
self.lf = ngs.LinearForm(self.fes)
if self.time_init_params.attributes.get("solver_type")=="convex_lin":
self.bf += coeff_bdf[0] * u * v *ngs.CoefficientFunction((ngs.x)) *ngs.dx
self.bf += self.dt* ngs.grad(u) * ngs.grad(v)*ngs.CoefficientFunction((ngs.x))* ngs.dx
self.bf += self.dt* u * v*ngs.CoefficientFunction((1/(ngs.x))) *ngs.dx
self.lf += - self.dt * ngs.grad(self.gfu[1])* ngs.grad(v)*ngs.CoefficientFunction((ngs.x)) * ngs.dx
self.lf += - self.dt* self.gfu[1] * v*ngs.CoefficientFunction((1/(ngs.x))) *ngs.dx
self.lf += self.dt* ngs.CoefficientFunction((1/(ngs.x))) *self.gfu_extrap* v *ngs.dx
self.lf += - self.dt * ngs.CoefficientFunction((1/(ngs.x))) * ngs.sin(2*self.gfu_extrap)/(2) * v * ngs.dx
for i in range(1,len(coeff_bdf)):
self.lf += -coeff_bdf[i]*self.gfu[i]*v*ngs.CoefficientFunction((ngs.x))*ngs.dx
elif self.time_init_params.attributes.get("solver_type")=="max_lin":
self.bf += coeff_bdf[0] * u * v * ngs.dx
self.bf += self.dt* ngs.grad(u) * ngs.grad(v)* ngs.dx
self.bf += -self.dt * ngs.CoefficientFunction((1/ngs.x)) * ngs.grad(u) * v* ngs.dx
self.bf += self.dt * ngs.CoefficientFunction((1/(ngs.x*ngs.x))) * ngs.sin(2*self.gfu_extrap)/(2*self.gfu_extrap)*u * v* ngs.dx
self.lf += - self.dt * ngs.grad(self.gfu[1])* ngs.grad(v) * ngs.dx
self.lf += self.dt * ngs.CoefficientFunction((1/ngs.x)) * ngs.grad(self.gfu[1]) * v *ngs.dx
self.lf += -self.dt * ngs.CoefficientFunction((1/(ngs.x*ngs.x))) * ngs.sin(2*self.gfu_extrap)/(2*self.gfu_extrap) * self.gfu[1] * v * ngs.dx
for i in range(1,len(coeff_bdf)):
self.lf += -coeff_bdf[i]*self.gfu[i]*v*ngs.dx
else:
raise RuntimeError("solver type ",self.time_init_params.attributes.get("solver_type"), " not supported. Use convex_lin or max_lin!")
def Solve(self):
self.Solve2(self.time_init_params.attributes.get("tstart"),self.time_init_params.attributes.get("tend"),self.time_init_params.attributes.get("BDF"))
def Solve2(self,tstart,tend,use_bdf,one_step=False):
if use_bdf == 1:
coeff_bdf = [1.0,0.0] # [1, -1] are the actual coeff, but we homogenize the dirichlet boundary condition by solving for dgfu := gfu[0] - gfu[1]
coeff_extrap = [1]
elif use_bdf == 2:
coeff_bdf =[1.5,-0.5,0.5] # [1.5, -2,0.5] same as bdf1
coeff_extrap = [2,-1]
else:
raise RuntimeError("Only BDF 1 and 2 supported!")
new_tstep = self.dt.Get()
number_elements = self.fes.ndof -1 # linear fes
new_mesh = self.mesh
new_fes = self.fes
new_gfu = self.gfu
new_gfu_extrap = self.gfu_extrap
t = tstart
## first steps needed to do multi-step bdf
bdf_counter = use_bdf
while bdf_counter>1 :
bdf_counter = bdf_counter -1
self.Solve2(t,tend,bdf_counter,True)
t = t+self.dt.Get()
self.ConstructForms(use_bdf)
debug_mode = self.output_params.attributes.get("debug")
use_adaptivity = self.adaptivity_params.attributes.get("use_adaptivity")
counter = 0
tsteps_arr = []
counter_tmin = 0
while abs(tend-t)>1e-12 and tend>t:
t = t+self.dt.Get()
counter = counter + 1
if debug_mode:
print("\rt=",t,end="")
for i in range(len(self.gfu)-1,0,-1):
self.gfu[i].vec.data = self.gfu[i-1].vec
self.gfu_extrap.vec[:]=0
for i in range(len(coeff_extrap)):
self.gfu_extrap.vec.data +=coeff_extrap[i]*self.gfu[i+1].vec
with ngs.TaskManager():
self.bf.Assemble()
self.lf.Assemble()
self.gfu[0].vec.data += self.bf.mat.Inverse(freedofs=self.fes.FreeDofs())*self.lf.vec
if use_adaptivity == 1:
self.blowuprate.value = self.adaptivity_params.attributes.get("blowup_rate")
if use_bdf > 1:
raise RuntimeError("Only BDF1 supported in adaptive mode!")
new_mesh,new_fes,new_gfu,new_gfu_extrap,number_elements,new_tstep,new_tcurrent,refined_something,counter_tmin = SpaceTimeAdaptivity(self.mesh,self.gfu,self.gfu_extrap,self.fes,self.dt.Get(),t,tend,len(self.gfu[0].vec)-1,self.adaptivity_params,counter_tmin)
mesh_size_el = ngs.Integrate(ngs.CoefficientFunction(1)*ngs.dx,self.mesh,element_wise=True)
min_mesh_size = min(mesh_size_el)
if min_mesh_size < self.adaptivity_params.attributes.get("eps_mesh_size"):
print("\nReach minimal mesh size of ", min_mesh_size)
break
self.dt.Set(new_tstep)
if refined_something:
self.mesh = new_mesh
self.fes = ngs.H1(self.mesh,order=1,dirichlet="left|right")
self.gfu = []
for i in range(0,len(new_gfu)):
self.gfu.append(ngs.GridFunction(self.fes))
self.gfu[i].Interpolate(new_gfu[i])
self.gfu_extrap = ngs.GridFunction(self.fes)
self.gfu_extrap.Interpolate(new_gfu_extrap)
t = new_tcurrent
self.ConstructForms(use_bdf)
counter = counter - 1
continue
if self.blowuprate.value>0:
self.blowuprate.dgfudx.append((self.gfu[0])(self.mesh(self.mesh.ngmesh.Points()[2][0],0.0,0.0))/self.mesh.ngmesh.Points()[2][0] )
self.blowuprate.time.append(t)
if debug_mode and use_adaptivity:
print("\n ndof = ", self.fes.ndof)
print("\n timestep = ",self.dt.Get())
if one_step:
break
if debug_mode and use_adaptivity:
print(" Number of time steps needed = ",counter)
print("\n ndof = ", self.fes.ndof)
def ComputeH1Error(gfu_1,gfu_2,mesh_ref):
with ngs.TaskManager():
l2error = ngs.sqrt(ngs.Integrate((gfu_1-gfu_2)**2*ngs.dx,mesh_ref))
graderror = ngs.sqrt(ngs.Integrate((ngs.grad(gfu_1)-ngs.grad(gfu_2))**2*ngs.dx,mesh_ref))
h1error =ngs.sqrt(graderror**2+l2error**2)
#print("L2 error = ",l2error)
#print("H1 error = ",h1error)
return l2error, h1error
def ComputeH1RAndEnergyError(gfu_1,gfu_2,mesh_ref):
with ngs.TaskManager():
l2rerror = ngs.sqrt(ngs.Integrate((gfu_1-gfu_2)**2*ngs.CoefficientFunction((ngs.x))*ngs.dx,mesh_ref))
h1rerror = ngs.sqrt(ngs.Integrate(((ngs.grad(gfu_1)-ngs.grad(gfu_2))**2*ngs.CoefficientFunction((ngs.x))+(gfu_1-gfu_2)**2*ngs.CoefficientFunction((1/ngs.x)))*ngs.dx,mesh_ref))
energyerror = ngs.sqrt(ngs.Integrate(((ngs.grad(gfu_1)-ngs.grad(gfu_2))**2*ngs.CoefficientFunction((ngs.x))+ngs.sin((gfu_1-gfu_2))**2*ngs.CoefficientFunction((1/ngs.x)))*ngs.dx,mesh_ref))
#print("L2R error = ",l2rerror)
#print("H1R error = ",h1rerror)
#print("Energy error = ",energyerror)
return l2rerror, h1rerror, energyerror
## Add project root to sys.path
import sys
import os
sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
from RSHMHF_FE import parameters_RSHMHF
from RSHMHF_FE.solving_RSHMHF import *
from HMHF_FE.Dim2d import parameters_HMHF2D
from HMHF_FE.Dim2d.solving_HMHF2D import *
from HMHF_FE.Dim2d.spherical_transf import *
import os
print("############### Compute reference solution from RSHMHF ####################")
## initial condition for RSHMHF
## requirement: u_0(x=0) = 0
u0 = ngs.CoefficientFunction((0.5*math.pi*ngs.x*ngs.x))
## parameter input file
this_dir = os.path.dirname(__file__)
input_file = os.path.join(this_dir, "input_files/RSHMHF_parameters.txt")
time_init_params_ref,space_init_params_ref,output_params_ref, adaptivity_params_ref =parameters_RSHMHF.ReadAndPrintAllParametersFromFile(input_file,0)
time_init_params_ref.attributes.get("solver_type")
number_elements_ref = space_init_params_ref.attributes.get("number_elements")
use_bdf_ref = time_init_params_ref.attributes.get("BDF")
dt_ref = time_init_params_ref.attributes.get("tstep")
fes_order = space_init_params_ref.attributes.get("fes_order")
print("use bdf ", use_bdf_ref)
print("time steps = ",dt_ref)
print("ref number elements = ",number_elements_ref)
print("fes_order = ",fes_order)
print("use adaptivity = ",adaptivity_params_ref.attributes.get("use_adaptivity"))
a = space_init_params_ref.attributes.get("interval_start")
b = space_init_params_ref.attributes.get("interval_end")
mesh_ref = Create1DMesh(number_elements_ref,a,b)
use_dirichlet = space_init_params_ref.attributes.get("Dirichlet")
if use_dirichlet:
## by default: "left" and "right"
dirichlet_str = mesh_ref.GetBoundaries()[0] + "|" + mesh_ref.GetBoundaries()[1]
#print(dirichlet_str)
V_ref = ngs.H1(mesh_ref,order=fes_order,dirichlet=dirichlet_str)
#print(V_ref.FreeDofs())
else:
raise RuntimeError("Only with Dirichlet boundary condition!")
print("ndof = ",V_ref.ndof)
params_ref = {"time":time_init_params_ref,"output":output_params_ref,"adaptivity":adaptivity_params_ref}
params_ref["time"].attributes["solver_type"] = "max_lin"
solver_RSHMHF_ref = SolvingRSHMHF(mesh_ref,V_ref,u0,params_ref)
solver_RSHMHF_ref.Solve()
print("\n############### Compute solution HMHF 2D ####################")
## input file
this_dir = os.path.dirname(__file__)
input_file = os.path.join(this_dir, "input_files/HMHF_2D_PPFEM_parameters.txt")
params = parameters_HMHF2D.ReadAndPrintAllParametersFromFile(input_file,0)
time_init_params_ref = params.get("time")
space_init_params_ref = params.get("space")
output_params_ref = params.get("output")
adaptivity_params_ref = params.get("adaptivity")
fes_order = space_init_params_ref.attributes.get("fes_order")
print("fes order = ",fes_order)
print("domain = ", space_init_params_ref.attributes.get("domain"))
print("use adaptivity = ",adaptivity_params_ref.attributes.get("use_adaptivity"))
u0_1 = ngs.CoefficientFunction((ngs.x/ngs.sqrt(ngs.x*ngs.x+ngs.y*ngs.y)*ngs.sin(0.5*math.pi*(ngs.x*ngs.x+ngs.y*ngs.y) )))
u0_2 = ngs.CoefficientFunction((ngs.y/ngs.sqrt(ngs.x*ngs.x+ngs.y*ngs.y)*ngs.sin(0.5*math.pi*(ngs.x*ngs.x+ngs.y*ngs.y) )))
u0_3 = ngs.CoefficientFunction((ngs.cos(0.5*math.pi*(ngs.x*ngs.x+ngs.y*ngs.y) )))
u0 = ngs.CoefficientFunction((u0_1,u0_2,u0_3))
h_mesh = space_init_params_ref.attributes.get("mesh_size")
print("")
print("###### tau convergence ######")
print("mesh size = ", space_init_params_ref.attributes.get("mesh_size"))
print("Spherical transformation...")
tau_arr = [1/20,1/40,1/80,1/160,1/320]
mesh_2D_ref, fes_2D_ref, gfu_2D_ref, gradgfu_2D_ref = SphericalTransform(solver_RSHMHF_ref.gfu[0],solver_RSHMHF_ref.mesh,h_mesh,u0,fes_order)
print("...done!")
print("use bdf ", time_init_params_ref.attributes.get("BDF"))
solver_type = ["PPFEM","TFEM"]
for s in solver_type:
print("")
print("\nuse solver ",s)
l2errors = []
h1errors = []
params["time"].attributes["solver_type"] = s
for tau in tau_arr:
solver_2D = SolvingHMHF2D(u0,params)
solver_2D.ChangeTimeStep(tau)
print("\ntime step = ",tau)
solver_2D.Solve()
if s == "TFEM":
gradgfu = ngs.CoefficientFunction((ngs.grad(solver_2D.gfu[0].components[0])[0],ngs.grad(solver_2D.gfu[0].components[1])[0],ngs.grad(solver_2D.gfu[0].components[2])[0],ngs.grad(solver_2D.gfu[0].components[0])[1],ngs.grad(solver_2D.gfu[0].components[1])[1],ngs.grad(solver_2D.gfu[0].components[2])[1]), dims=(2,3))
l2error,h1error = ComputeH1Error2DForCF(gfu_2D_ref,ngs.CoefficientFunction((solver_2D.gfu[0].components[0],solver_2D.gfu[0].components[1],solver_2D.gfu[0].components[2])),gradgfu_2D_ref,gradgfu,solver_2D.mesh)
else:
l2error,h1error = ComputeH1Error2DForCF(gfu_2D_ref,solver_2D.gfu[0],gradgfu_2D_ref,ngs.grad(solver_2D.gfu[0]),solver_2D.mesh)
l2errors.append(l2error)
h1errors.append(h1error)
print("\nSummary")
print("time steps = ",tau_arr)
eoc_l2 = [ngs.log(l2errors[i+1]/l2errors[i])/ngs.log(0.5) for i in range(len(l2errors)-1)]
eoc_h1 = [ngs.log(h1errors[i+1]/h1errors[i])/ngs.log(0.5) for i in range(len(h1errors)-1)]
print("L2-errors = ",l2errors)
print("L2 eoc = ",eoc_l2)
print("")
print("H1-errors = ",h1errors)
print("H1 eoc = ",eoc_h1)
print("")
print("###### mesh size convergence ######")
h_arr = [1/2,1/4,1/8,1/16,1/32,1/64]
print("time steps = ", time_init_params_ref.attributes.get("tstep"))
for s in solver_type:
print("")
print("\nuse solver ",s)
params["time"].attributes["solver_type"] = s
l2errors = []
h1errors = []
for h_loop in h_arr:
params["space"].attributes["mesh_size"] = h_loop
print("\nmesh size = ",h_loop)
print("Spherical transformation...")
mesh_2D_ref, fes_2D_ref, gfu_2D_ref, gradgfu_2D_ref = SphericalTransform(solver_RSHMHF_ref.gfu[0],solver_RSHMHF_ref.mesh,h_loop,u0,fes_order)
print("...done!")
solver_2D = SolvingHMHF2D(u0,params)
solver_2D.Solve()
if s == "TFEM":
gradgfu = ngs.CoefficientFunction((ngs.grad(solver_2D.gfu[0].components[0])[0],ngs.grad(solver_2D.gfu[0].components[1])[0],ngs.grad(solver_2D.gfu[0].components[2])[0],ngs.grad(solver_2D.gfu[0].components[0])[1],ngs.grad(solver_2D.gfu[0].components[1])[1],ngs.grad(solver_2D.gfu[0].components[2])[1]), dims=(2,3))
l2error,h1error = ComputeH1Error2DForCF(gfu_2D_ref,ngs.CoefficientFunction((solver_2D.gfu[0].components[0],solver_2D.gfu[0].components[1],solver_2D.gfu[0].components[2])),gradgfu_2D_ref,gradgfu,solver_2D.mesh)
else:
l2error,h1error = ComputeH1Error2DForCF(gfu_2D_ref,solver_2D.gfu[0],gradgfu_2D_ref,ngs.grad(solver_2D.gfu[0]),solver_2D.mesh)
l2errors.append(l2error)
h1errors.append(h1error)
print("\nSummary")
print("mesh sizes = ",h_arr)
eoc_l2 = [ngs.log(l2errors[i+1]/l2errors[i])/ngs.log(0.5) for i in range(len(l2errors)-1)]
eoc_h1 = [ngs.log(h1errors[i+1]/h1errors[i])/ngs.log(0.5) for i in range(len(h1errors)-1)]
print(l2errors)
print("L2-errors = ",l2errors)
print("L2 eoc = ",eoc_l2)
print("")
print("H1-errors = ",h1errors)
print("H1 eoc = ",eoc_h1)
## Add project root to sys.path
import sys
import os
sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
from RSHMHF_FE.solving_RSHMHF import *
import os
import numpy as np
## initial condition
## requirement: u_0(x=0) = 0
u0 = ngs.CoefficientFunction((math.pi*ngs.x*(1-ngs.x)))
print("############### Compute reference solution ####################")
## input file
this_dir = os.path.dirname(__file__)
input_path = os.path.join(this_dir, "input_files/RSHMHF_parameters.txt")
time_init_params_ref,space_init_params_ref,output_params_ref, adaptivity_params_ref = ReadAndPrintAllParametersFromFile(input_path,0)
number_elements_ref = space_init_params_ref.attributes.get("number_elements")
use_bdf_ref = time_init_params_ref.attributes.get("BDF")
dt_ref = time_init_params_ref.attributes.get("tstep")
fes_order_ref = space_init_params_ref.attributes.get("fes_order")
## ToDo: Put all printing into some outputfile with output options
print("use bdf ", use_bdf_ref)
print("time steps = ",dt_ref)
print("ref number elements = ",number_elements_ref)
print("fes order = ",fes_order_ref)
print("use linearization = ", time_init_params_ref.attributes.get("solver_type"))
print("use adaptivity = ",adaptivity_params_ref.attributes.get("use_adaptivity"))
a = space_init_params_ref.attributes.get("interval_start")
b = space_init_params_ref.attributes.get("interval_end")
mesh_ref = Create1DMesh(number_elements_ref,a,b)
use_dirichlet = space_init_params_ref.attributes.get("Dirichlet")
if use_dirichlet:
## by default: "left" and "right"
dirichlet_str = mesh_ref.GetBoundaries()[0] + "|" + mesh_ref.GetBoundaries()[1]
#print(dirichlet_str)
V_ref = ngs.H1(mesh_ref,order=fes_order_ref,dirichlet=dirichlet_str)
#print(V_ref.FreeDofs())
else:
raise RuntimeError("Only with Dirichlet boundary condition!")
print("ndof = ",V_ref.ndof)
params_ref = {"time":time_init_params_ref,"output":output_params_ref,"adaptivity":adaptivity_params_ref}
solver_RSHMHF_ref = SolvingRSHMHF(mesh_ref,V_ref,u0,params_ref)
solver_RSHMHF_ref.Solve()
print("\n############### Compute solution ####################")
print("")
time_init_params,space_init_params,output_params, adaptivity_params = ReadAndPrintAllParametersFromFile(input_path,0)
use_bdf = 1
time_init_params.attributes["BDF"] = use_bdf
print("###### Mesh size convergence ######")
print("use bdf ", use_bdf)
solver_type = ["convex_lin","max_lin"]
fes_order = 1
space_init_params.attributes["fes_order"] = fes_order
print("\nuse fes order ", fes_order)
nmb_elements = [2,4,8,16]
for s in solver_type:
l2rerrors = []
h1rerrors = []
energyerrors = []
time_init_params.attributes["solver_type"] = s
print("\nuse solver type ",s)
for i in nmb_elements:
space_init_params.attributes["number_elements"] = i
use_bdf = time_init_params.attributes.get("BDF")
print("\nmesh size = ",1/i)
mesh = Create1DMesh(i,0,1)
use_dirichlet = space_init_params.attributes.get("Dirichlet")
if use_dirichlet:
dirichlet_str = mesh.GetBoundaries()[0] + "|" + mesh.GetBoundaries()[1]
V = ngs.H1(mesh,order=fes_order,dirichlet=dirichlet_str)
#print(V_ref.FreeDofs())
else:
raise RuntimeError("Only with Dirichlet boundary condition for now!")
params = {"time":time_init_params,"output":output_params,"adaptivity":adaptivity_params}
solver_RSHMHF = SolvingRSHMHF(mesh,V,u0,params)
solver_RSHMHF.Solve()
l2rerror,h1rerror,energyerror = ComputeH1RAndEnergyError(solver_RSHMHF.gfu[0],solver_RSHMHF_ref.gfu[0],solver_RSHMHF.mesh)
l2rerrors.append(l2rerror)
h1rerrors.append(h1rerror)
energyerrors.append(energyerror)
eoc_l2r = [ngs.log(l2rerrors[i+1]/l2rerrors[i])/ngs.log(0.5) for i in range(len(l2rerrors)-1)]
eoc_h1r = [ngs.log(h1rerrors[i+1]/h1rerrors[i])/ngs.log(0.5) for i in range(len(h1rerrors)-1)]
eoc_energy = [ngs.log(energyerrors[i+1]/energyerrors[i])/ngs.log(0.5) for i in range(len(energyerrors)-1)]
print("\nSummary")
print("mesh sizes = ",[1/2,1/4,1/8,1/16])
print("L2-errors = ",l2rerrors)
print("L2 eoc = ",eoc_l2r)
print("")
print("H1-errors = ",h1rerrors)
print("H1 eoc",eoc_h1r)
print("")
print("Energy-errors = ",energyerrors)
print("Energy eoc = ",eoc_energy)
print("\n###### Time step convergence ######")
print("")
time_init_params.attributes["BDF"] = 1
space_init_params.attributes["number_elements"] =number_elements_ref
dt_arr = [1.25e-2,6.25e-3,3.125e-3,3.125e-3/2,3.125e-3/4]
space_init_params.attributes["fes_order"] = fes_order
print("use fes order = ",fes_order)
print("use bdf ", use_bdf)
print("use mesh size = ",1/number_elements_ref)
mesh = Create1DMesh(number_elements_ref,0,1)
use_dirichlet = space_init_params.attributes.get("Dirichlet")
if use_dirichlet:
dirichlet_str = mesh.GetBoundaries()[0] + "|" + mesh.GetBoundaries()[1]
V = ngs.H1(mesh,order=fes_order,dirichlet=dirichlet_str)
#print(V_ref.FreeDofs())
else:
raise RuntimeError("Only with Dirichlet boundary condition for now!")
for s in solver_type:
l2rerrors = []
h1rerrors = []
energyerrors = []
time_init_params.attributes["solver_type"] = s
print("\nuse solver type ",s)
for dt in dt_arr:
time_init_params.attributes["tstep"] = dt
time_step = time_init_params.attributes.get("tstep")
print("\ntime steps = ",time_step)
params = {"time":time_init_params,"output":output_params,"adaptivity":adaptivity_params}
solver_RSHMHF = SolvingRSHMHF(mesh,V,u0,params)
solver_RSHMHF.Solve()
l2rerror,h1rerror,energyerror = ComputeH1RAndEnergyError(solver_RSHMHF.gfu[0],solver_RSHMHF_ref.gfu[0],solver_RSHMHF.mesh)
l2rerrors.append(l2rerror)
h1rerrors.append(h1rerror)
energyerrors.append(energyerror)
eoc_l2r = [ngs.log(l2rerrors[i+1]/l2rerrors[i])/ngs.log(0.5) for i in range(len(l2rerrors)-1)]
eoc_h1r = [ngs.log(h1rerrors[i+1]/h1rerrors[i])/ngs.log(0.5) for i in range(len(h1rerrors)-1)]
eoc_energy = [ngs.log(energyerrors[i+1]/energyerrors[i])/ngs.log(0.5) for i in range(len(energyerrors)-1)]
print("\nSummary")
print("time steps = ",dt_arr)
print("L2-errors = ",l2rerrors)
print("L2 eoc = ",eoc_l2r)
print("")
print("H1-errors = ",h1rerrors)
print("H1 eoc = ",eoc_h1r)
print("")
print("Energy-errors = ",energyerrors)
print("Energy eoc = ",eoc_energy)
# for now only unit disk allowed
domain unit_disk
# number of elements in interval I
mesh_size 0.00390625
# time parameters
tstart 0.0
tstep 1e-6
tend 0.1
# time stepping scheme
BDF 1
# 1 = true, 0 = false
# if 0, then zero Neumann used
Dirichlet 1
Dirichlet_bnd_name unit_circle
fes_order 1
# debug mode: print to terminal
debug 1
# PPFEM, TFEM or CPFEM
# CPFEM only convergence in mesh size feasible
solver_type PPFEM
# define the intervall in 1D: I := [interval_start,interval_end]
interval_start 0
interval_end 1
# number of elements in interval I, h = 1/number_elements
number_elements 16384
#number_elements 10
# different linearizations for non-linearity in RSHMHF
solver_type convex_lin
# time parameters
tstart 0.0
tstep 1e-6
#tstep 1e-2
tend 0.1
# time stepping scheme
BDF 2
fes_order 2
# 1 = true, 0 = false
# if 0, then zero Neumann used
Dirichlet 1
# debug mode: print to terminal
debug 1
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment